principal_distance_from_F

% principal_distance_from_F: determina single c from fundamental matrix

 c = principal_distance_from_F(F)

 F = fundamental matrix, rank 2-matrix
     assumed for conditioned cooridates
 
 c = principal distance

 Reference:
 P. Sturm, Z.L. Cheng, P.C.Y. Chen, A.N. Poo 
 Focal length calibration from two views:
 method and analysis of singular cases
 Computer Vision and Image Understanding 99 (2005) 58–95

 Wolfgang Förstner 3/2013
 wfoerstn@uni-bonn.de

0001 %% principal_distance_from_F: determina single c from fundamental matrix
0002 %
0003 % c = principal_distance_from_F(F)
0004 %
0005 % F = fundamental matrix, rank 2-matrix
0006 %     assumed for conditioned cooridates
0007 %
0008 % c = principal distance
0009 %
0010 % Reference:
0011 % P. Sturm, Z.L. Cheng, P.C.Y. Chen, A.N. Poo
0012 % Focal length calibration from two views:
0013 % method and analysis of singular cases
0014 % Computer Vision and Image Understanding 99 (2005) 58–95
0015 %
0016 % Wolfgang Förstner 3/2013
0017 % wfoerstn@uni-bonn.de
0018 
0019 
0020 function c = principal_distance_from_F(F)
0021 
0022 % enforce long principal distance as prior
0023 factor = 1;
0024 K0 = diag([factor,factor,1]);
0025 F  = K0*F*K0;
0026 
0027 % svd
0028 [U,D,V] = svd(F);
0029 a = D(1,1);
0030 b = D(2,2);
0031 
0032 %% build coefficients
0033 
0034 % 1. linear equation
0035 a1 = a*U(3,1)*U(3,2)*(1-V(3,1)^2)+b*V(3,1)*V(3,2)*(1-U(3,2)^2);
0036 a0 = U(3,2)*V(3,1)*(a*U(3,1)*V(3,1)+b*U(3,2)*V(3,2));
0037 
0038 % 2. linear equation
0039 b1 = a*V(3,1)*V(3,2)*(1-U(3,1)^2)+b*U(3,1)*U(3,2)*(1-V(3,2)^2);
0040 b0 = U(3,1)*V(3,2)*(a*U(3,1)*V(3,1)+b*U(3,2)*V(3,2));
0041 
0042 % quadratic equation
0043 c2 = a^2*(1-U(3,1)^2)*(1-V(3,1)^2)-b^2*(1-U(3,2)^2)*(1-V(3,2)^2);
0044 c1 = a^2*(U(3,1)^2+V(3,1)^2-2*U(3,1)^2*V(3,1)^2)-...
0045      b^2*(U(3,2)^2+V(3,2)^2-2*U(3,2)^2*V(3,2)^2);
0046 c0 = a^2*U(3,1)^2*V(3,1)^2-b^2*U(3,2)^2*V(3,2)^2;
0047 
0048 f_squares = roots([c2,c1,c0]);
0049 % fs = sqrt(f_squares);
0050 
0051 % flin1 = -a0/a1;
0052 % flin2 = -b0/b1;
0053 
0054 c = 0;
0055 for i = [1,2]
0056     f2 = f_squares(i);
0057     if isreal(f2) && f2 > 0
0058         if abs(a1*f2+a0) < 10^(-10) && abs(b1*f2+b0) < 10^(-10)
0059            c = sqrt(f2);
0060         end
0061     end
0062 end
0063 
0064 c = factor*c;
0065 
0066