Triangulación y transformación lineal directa

Question

Siguiendo la Geometría Multiview de Hartley/Zisserman, Algoritmo 12: El método de triangulación óptimo (p318), obtuve los puntos de imagen correspondientes xhat1 y xhat2 (paso 10). En el paso 11, uno necesita calcular el punto 3D Xhat. Uno de estos métodos es la Transformada Lineal Directa (DLT), mencionada en 12.2 (p312) y 4.1 (p88).

El método homogéneo (DLT), p312-313, afirma que encuentra una solución como el vector singular unidad correspondiente al valor singular más pequeño de A, por lo tanto,

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.'); 

Sin embargo, A es un 16x1 matriz, lo que resulta en una Va que es 1x1.

¿Qué estoy haciendo mal (y una solución) para obtener el punto 3D?

Por lo que su datos de ejemplo vale:

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 

actualización código de trabajo para la sección XI en el 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

Como se menciona en el libro (s 12.2), p^{i T} son las filas de P. Por lo tanto, no es necesario transponer P1(k,:) (es decir, la formulación correcta es A = [xhat1(1) * P1(3,:) - P1(1,:) ; ...).

Espero que haya sido solo un error tipográfico.

Además, se recomienda para normalizar cada fila de A con su norma L2, es decir, para todos los i

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

Y si desea representar los puntos 3D triangulados, hay que normalizar Xhat antes el trazado (su sentido de otra manera), es decir,

Xhat = Xhat/Xhat(4);