For future reference using Matlab / Octave syntax!
The join of two points in Plücker coordinates can be expressed as follows
% line = point join point
function L=join(A, B)
L=[
A(1)*B(2)-A(2)*B(1);
A(1)*B(3)-A(3)*B(1);
A(1)*B(4)-A(4)*B(1);
A(2)*B(3)-A(3)*B(2);
A(2)*B(4)-A(4)*B(2);
A(3)*B(4)-A(4)*B(3)
];
end % function
These are the 6 distinct values from the anti-symmetric matrix
Lx=B*A'-A*B'
A point on the backprojection ray can be found
X=pinv(P)*x
where
x=[u v 1]'
is the image point at pixel location (u,v) and
pinv(P)
the pseudoinverse of the projection matrix.
The camera center can be found as the null space of the projection matrix
C=null(P);
C=C/C(4)
The Plücker coordinates of the backprojection ray are thus
L=join(X,C)
For those interested in oriented projective geometry: If you normalize the projection matrix as follows
% Get length of principal ray
m3n=norm(P(3,1:3));
% Enforce positivity of determinant
if (det(P(:,1:3))<0)
m3n=-m3n;
end % if
% Normalize
P=P/m3n;
Then the determinant of the left 3x3 matrix is positive (i.e. right handed systems) and L will point from C to X.
Update: The other day, a co-worker asked me for a matrix-form solution. If you multiply the following matrix with an image point [x0 x1 x2]' you get the Plücker coordinates directly:
A=[
- p12*p23 + p13*p22, + p02*p23 - p03*p22, + p03*p12 - p02*p13
+ p11*p23 - p13*p21, - p01*p23 + p03*p21, + p01*p13 - p03*p11
- p11*p22 + p12*p21, + p01*p22 - p02*p21, + p02*p11 - p01*p12
- p10*p23 + p13*p20, + p00*p23 - p03*p20, + p03*p10 - p00*p13
+ p10*p22 - p12*p20, - p00*p22 + p02*p20, + p00*p12 - p02*p10
- p10*p21 + p11*p20, + p00*p21 - p01*p20, + p01*p10 - p00*p11
]
The line direction is therefore
d =
p01*p12*x2 - p02*p11*x2 - p01*p22*x1 + p02*p21*x1 + p11*p22*x0 - p12*p21*x0
p02*p10*x2 - p00*p12*x2 + p00*p22*x1 - p02*p20*x1 - p10*p22*x0 + p12*p20*x0
p00*p11*x2 - p01*p10*x2 - p00*p21*x1 + p01*p20*x1 + p10*p21*x0 - p11*p20*x0
and the line moment is:
m =
p03*p10*x2 - p00*p13*x2 + p00*p23*x1 - p03*p20*x1 - p10*p23*x0 + p13*p20*x0
p03*p11*x2 - p01*p13*x2 + p01*p23*x1 - p03*p21*x1 - p11*p23*x0 + p13*p21*x0
p03*p12*x2 - p02*p13*x2 + p02*p23*x1 - p03*p22*x1 - p12*p23*x0 + p13*p22*x0
Obviously, one can pull out the [x2 x1 x0] from the expression to get matrix form. The following MATLAB symbolic code was used to generate this solution:
syms p00 p01 p02 p03 p10 p11 p12 p13 p20 p21 p22 p23 real
P=[ p00 p01 p02 p03
p10 p11 p12 p13
p20 p21 p22 p23
]
% Some Image Point
syms x0 x1 x2 real
x=[x0 x1 x2]'
% Source Position
C=null(P)
% Backprojection of x
Xtilde=pinv(P)*x
% Backprojection line (Plücker Matrix)
Lx=simplify(C*Xtilde' - Xtilde*C')
% It's homogeneous...
arbitrary_scale=(p00*p11*p22 - p00*p12*p21 - p01*p10*p22 + p01*p12*p20 + p02*p10*p21 - p02*p11*p20)
Lx=Lx*arbitrary_scale
% Plücker Coordinates:
L=[Lx(1,2) Lx(1,3) Lx(1,4) Lx(2,3) Lx(2,4) Lx(3,4)];
% Direction
d=[-L(3) -L(5) -L(6)]';
% Moment
m=[L(4) -L(2) L(1)]';
% In matrix form
A=[
- p12*p23 + p13*p22, + p02*p23 - p03*p22, + p03*p12 - p02*p13
+ p11*p23 - p13*p21, - p01*p23 + p03*p21, + p01*p13 - p03*p11
- p11*p22 + p12*p21, + p01*p22 - p02*p21, + p02*p11 - p01*p12
- p10*p23 + p13*p20, + p00*p23 - p03*p20, + p03*p10 - p00*p13
+ p10*p22 - p12*p20, - p00*p22 + p02*p20, + p00*p12 - p02*p10
- p10*p21 + p11*p20, + p00*p21 - p01*p20, + p01*p10 - p00*p11
]
% Verification: (should be zero)
simplify(A*x - L')