Triangolazione e trasformazione lineare diretta

Question

Segue la geometria multivista di Hartley/Zisserman, Algoritmo 12: il metodo di triangolazione ottimale (p318), ho ottenuto i punti immagine corrispondenti xhat1 e xhat2 (passaggio 10). Nel passaggio 11, è necessario calcolare il punto 3D Xhat. Uno di questi metodi è Direct Linear Transform (DLT), menzionato in 12.2 (p312) e 4.1 (p88).

Il metodo omogeneo (DLT), p312-313, afferma che trova una soluzione come unità singolare vettore corrispondente al più piccolo valore singolare di A, così,

A = [xhat1(1) * P1(3,:)' - P1(1,:)' ; 
     xhat1(2) * P1(3,:)' - P1(2,:)' ; 
     xhat2(1) * P2(3,:)' - P2(1,:)' ; 
     xhat2(2) * P2(3,:)' - P2(2,:)' ]; 

[Ua Ea Va] = svd(A); 
Xhat = Va(:,end); 

plot3(Xhat(1),Xhat(2),Xhat(3), 'r.'); 

Tuttavia, A è un 16x1 matrice, risultante in un Va che è 1x1.

Cosa sto facendo male (e una correzione) nell'ottenere il punto 3D?

Per quello che i suoi dati di esempio vale la pena:

xhat1 = 

    1.0e+009 * 

    4.9973 
    -0.2024 
    0.0027 


xhat2 = 

    1.0e+011 * 

    2.0729 
    2.6624 
    0.0098 


P1 = 

    699.6674   0 392.1170   0 
     0 701.6136 304.0275   0 
     0   0 1.0000   0 


P2 = 

    1.0e+003 * 

    -0.7845 0.0508 -0.1592 1.8619 
    -0.1379 0.7338 0.1649 0.6825 
    -0.0006 0.0001 0.0008 0.0010 


A = <- my computation 

    1.0e+011 * 

    -0.0000 
     0 
    0.0500 
     0 
     0 
    -0.0000 
    -0.0020 
     0 
    -1.3369 
    0.2563 
    1.5634 
    2.0729 
    -1.7170 
    0.3292 
    2.0079 
    2.6624 

Aggiornamento codice di lavoro per il settore XI in algoritmo

% xi 
A = [xhat1(1) * P1(3,:) - P1(1,:) ; 
    xhat1(2) * P1(3,:) - P1(2,:) ; 
    xhat2(1) * P2(3,:) - P2(1,:) ; 
    xhat2(2) * P2(3,:) - P2(2,:) ]; 

A(1,:) = A(1,:)/norm(A(1,:)); 
A(2,:) = A(2,:)/norm(A(2,:)); 
A(3,:) = A(3,:)/norm(A(3,:)); 
A(4,:) = A(4,:)/norm(A(4,:)); 

[Ua Ea Va] = svd(A); 
X = Va(:,end); 
X = X/X(4); % 3D Point 

Answer 1

Come è menzionato nel libro (sec 12.2), p^{i T} sono le righe di P. Pertanto, non è necessario trasporre lo P1(k,:) (ad esempio, la formulazione corretta è A = [xhat1(1) * P1(3,:) - P1(1,:) ; ...).

Spero sia stato solo un errore di battitura.

Inoltre, si consiglia di normalizzare ogni fila di A con la sua norma L2, vale a dire per tutti i

A(i,:) = A(i,:)/norm(A(i,:));

E se si vuole tracciare i punti 3D triangolari, bisogna normalizzare Xhat prima tracciando (il suo senza senso contrario), vale a dire

Xhat = Xhat/Xhat(4);