X  = rand(4,13);

if 0
  P1 = [1 0 0 0;
	0 1 0 0;
	0 0 1 0];

  P2 = [1 0 0 0
	0 1 0 0
	0 0 1 1];
else
  P1 = rand(3,4);
  P2 = rand(3,4);
end

b  = P1(:,4)-P1(:,1:3)*inv(P2(:,1:3))*P2(:,4);

E  = [  0   -b(3)  b(2)
       b(3)   0   -b(1) 
      -b(2)  b(1)   0 ];

F  = E*P1(:,1:3)*inv(P2(:,1:3));

[U,D,V] = svd(F);

e1 = U(:,3);
e2 = V(:,3);

d    = 1./diag(sqrt(D));
d(3) = 1;

UU = [ 0   -d(2) 0
      d(1)   0   0
       0     0   1] * U';

VV = V * diag(d);

rP1 = [UU' [0;0;0]];
rP2 = [VV  e2];

u1 = P1*X;u1 = u1./(ones(3,1)*u1(3,:));
u2 = P2*X;u2 = u2./(ones(3,1)*u2(3,:));

ru1 = rP1*X;
ru2 = rP2*X;

for i=1:size(X,2)
  XX(:,i) = vgg_X_from_xP_lin([u1(1:2,i) u2(1:2,i)],{P1, P2});
end




