demo_check_of_signs_single_and_two_view_geometry

% demo: check of signs single and two view geometry

 1. Epipoles as projections and vectors of determinants 
 2. Signs of epipolar lines 
 3. Positivity of projection centre Z_h > 0
 4. l' e'+l'' e'' = 0, PCV (12.286)
 5. Sign of reconstructed 3D lines, PCV (13.288)

 Wolfgang Förstner
 wfoerstn@uni-bonn.de

 last changes: Susanne Wenzel 06/18
 wenzel@igg.uni-bonn.de

0001 %% demo: check of signs single and two view geometry
0002 %
0003 % 1. Epipoles as projections and vectors of determinants
0004 % 2. Signs of epipolar lines
0005 % 3. Positivity of projection centre Z_h > 0
0006 % 4. l' e'+l'' e'' = 0, PCV (12.286)
0007 % 5. Sign of reconstructed 3D lines, PCV (13.288)
0008 %
0009 % Wolfgang Förstner
0010 % wfoerstn@uni-bonn.de
0011 %
0012 % last changes: Susanne Wenzel 06/18
0013 % wenzel@igg.uni-bonn.de
0014 
0015 clc
0016 
0017 addpath(genpath('../General-Functions'))
0018 
0019 disp('----------------------------------------------------------------')
0020 disp('----- check Werner/Pajdla_s ls*es+lss*ess = 0 PCV (12.286) -----')
0021 disp('----------------------------------------------------------------')
0022 
0023 disp('Create Fundamental matrix F, assuming randomly given positive projection centres Z1 and Z2 ')
0024 % random Projection matrices with positive projection centres
0025 P1  = randn(3,4);                    % random P
0026 DP1 = det(P1(1:3,1:3));              % determinant of A
0027 P1  = P1/sign(DP1);                  % proper P, see PCV p. 474
0028                                      % --> positive projection centre
0029 % Projection centre as intersection of principal planes, see PCV (12.45)
0030 Z1  = -calc_Gammadual(calc_Pi(P1(1,:)')*P1(2,:)')'*P1(3,:)';
0031 
0032 % same for P2
0033 P2  = randn(3,4);
0034 DP2 = det(P2(1:3,1:3));
0035 P2  = P2/sign(DP2);
0036 Z2  = -calc_Gammadual(calc_Pi(P2(1,:)')*P2(2,:)')'*P2(3,:)';
0037 
0038 % corresponding Q's
0039 Q1 = calc_Q_from_P(P1);             
0040 Q2 = calc_Q_from_P(P2);
0041 
0042 % Fundamental matrix, see PCV (13.19)
0043 F = Q1*calc_Dual*Q2'                                                       %#ok<NOPTS>
0044 
0045 % epipoles, see PCV (13.71)
0046 e1 = P1*Z2;
0047 e2 = P2*Z1;
0048 
0049 % epipoles, see PCV (13.72)
0050 e1_det = [...
0051     det([P2(1,:)',P2(2,:)',P2(3,:)',P1(1,:)']);...
0052     det([P2(1,:)',P2(2,:)',P2(3,:)',P1(2,:)']);...
0053     det([P2(1,:)',P2(2,:)',P2(3,:)',P1(3,:)']);...
0054       ];
0055 
0056 disp(' ')
0057 disp('----- Check epipoles -----')  
0058 check_e1_minus_e1det_zero = (e1 - e1_det)';
0059 disp('Epipole e1 = P1 times Z2 (13.71) minus e1 from camera planes in the ')
0060 disp(['projection matrices (13.72) should be zero vector:   [', num2str(check_e1_minus_e1det_zero),']'])
0061 
0062 check_e1F_zero=(e1'*F)';
0063 check_Fe2_zero=(F*e2)';
0064 disp(' ')
0065 disp(['e1^T times F should be zero vector:   [', num2str(check_e1F_zero'),']^T'])
0066 disp(['F times e2 should be zero vector:     [', num2str(check_Fe2_zero),']'])
0067 
0068 disp(' ')
0069 disp('----- Create random 3D Line and according projections into images ...')
0070 % random 3D Line
0071 X1 = randn(4,1);            % two random 3D points
0072 X2 = randn(4,1);
0073 L  = calc_Pi(X1)*X2;        % resulting 3D line
0074 l1 = Q1*L;                  % lines in images
0075 l2 = Q2*L;
0076 
0077 disp('Object point X1 is projected into image 1 as x11, into image 2 as x12')
0078 disp('Object point X2 is projected into image 2 as x21, into image 2 as x22')
0079 % projections
0080 x11 = P1*X1;                % first point in images
0081 x12 = P2*X1;
0082 x21 = P1*X2;                % second point in images
0083 x22 = P2*X2;
0084 
0085 disp(' ')
0086 disp('----- Calculate epipolar lines l1 and l2 of X1 ...')
0087 % epipolar lines l1 and l2 of X1
0088 l1x = cross(e1,x11)                                                        %#ok<NOPTS>
0089 l1xF = F*x12;
0090 l2x = cross(e2,x12)                                                        %#ok<NOPTS>
0091 l2xF = F'*x11;
0092 
0093 disp(' ')
0094 disp('----- Check epipolar lines')
0095 check_l1x_minus_l1xF_zero = (l1x-l1xF)';
0096 check_l2x_minus_l2xF_zero = (l2x-l2xF)';
0097 disp(['epipolar line l1 = e1 cross x11 minus l1 = F times x12 should be zero vector:     [', ...
0098     num2str(check_l1x_minus_l1xF_zero),']'])
0099 disp(['epipolar line l2 = e2 cross x12 minus l2 = F^T times x11 should be zero vector:   [', ...
0100     num2str(check_l2x_minus_l2xF_zero),']'])
0101 
0102 check_l_on_x_1_zero = (l1x'*x11)';
0103 disp(['image point x11 should be on l1, thus l1^T times x11 should be zero:              ', num2str(check_l_on_x_1_zero)])
0104 check_e_on_x_1_zero = (l1x'*e1)';
0105 disp(['epipole e1 should be on l1, thus l1^T times e1 should be zero:                    ', num2str(check_e_on_x_1_zero)])
0106 
0107 check_l_on_x_2_zero = (l2x'*x12)';
0108 disp(['image point x12 should be on l2, thus l2^T times x12 should be zero:              ', num2str(check_l_on_x_2_zero)])
0109 check_e_on_x_2_zero = (l2x'*e2)';
0110 disp(['epipole e2 should be on l2, thus l2^T times e2 should be zero:                    ', num2str(check_e_on_x_2_zero)])
0111 
0112 %% start of check
0113 disp(' ')
0114 disp('---- Check positive sign and value of Z0 vs abs|det(A)| , see PCV p. 474 -----')
0115 Z_s_positive_Dets       = [Z1',abs(DP1);Z2',abs(DP2)]                      %#ok<NOPTS>
0116 
0117 
0118 disp('----- Check sign of l1*e1 and l2*e2 -----')
0119 lses   = l1'*e1;
0120 lssess = l2'*e2;
0121 check_lses_plus_lssess_zero = lses+lssess;
0122 disp(['PCV (13.287)   l1*e1 + l2*e2 should be zero:        ', num2str(check_lses_plus_lssess_zero)])
0123 check_QP_minus_PiZ_zero = calc_Dual*Q1'*P1-calc_Pi(Z1);
0124 disp(['PCV (13.288)   Q_dual(P1) - Pi(Z1) should be zero:  ', num2str(check_lses_plus_lssess_zero)])
0125 disp(' ')
0126 %%
0127 disp('Check sign of 3D line reconstruction PCV (13.288)')
0128 A1 = P1'*l1;
0129 A2 = P2'*l2;
0130 Lc = sign(l1'*e1)*calc_Pidual(A1)*A2;   
0131 disp( ['L = A2 cap A1 =                                    [',num2str(Lc'),']'])
0132 % Lc/Lc(1)*sign(Lc(1));
0133 % L/L(1)*sign(L(1));
0134 factor_12 = sign(l1'*e1);
0135 disp(['Sign of 3D line: ', num2str(factor_12)])
0136 check_L_and_Lreconstructed_zero = (Lc/Lc(1)*sign(Lc(1))-L/L(1)*sign(L(1)))';
0137 disp(['Diff L and reconstructed 3D line should be zero:   [', num2str(check_L_and_Lreconstructed_zero),']'])
0138 
0139 disp(' ')
0140 Lc = sign(l2'*e2)*calc_Pidual(A2)*A1;
0141 disp( ['L = A1 cap A2 =                                    [',num2str(Lc'),']'])
0142 % Lc/Lc(1)*sign(Lc(1));
0143 % L/L(1)*sign(L(1));
0144 factor_21 = sign(l2'*e2);
0145 disp(['Sign of 3D line: ', num2str(factor_21)])
0146 check_zero = (Lc/Lc(1)*sign(Lc(1))-L/L(1)*sign(L(1)))';
0147 disp(['Diff L and reconstructed 3D line should be zero:   [', num2str(check_zero),']'])
0148 
0149 
0150 
0151 
0152