demo_epipolar_geometry

% Epipolar geometry

 The script allows to determine the relative orientation (F-matrix) 
 of two images from 8 corresponding points and to show the 
 epipolar lines of points in the right image 
 togehter with their uncertainty.

 The uncertainty of the measured points can be chosen to be 
  1. the rounding error(appr. 1/3 pixel) or 
  2. derived from the image content. 
     Then the script can be used to explore 
     the uncertainty of matching based on the sstructure tensor of
     the image, see PCV (13.93)
 This is controlled by the variable 'cov_im'

0001 %% Epipolar geometry
0002 %
0003 % The script allows to determine the relative orientation (F-matrix)
0004 % of two images from 8 corresponding points and to show the
0005 % epipolar lines of points in the right image
0006 % togehter with their uncertainty.
0007 %
0008 % The uncertainty of the measured points can be chosen to be
0009 %  1. the rounding error(appr. 1/3 pixel) or
0010 %  2. derived from the image content.
0011 %     Then the script can be used to explore
0012 %     the uncertainty of matching based on the sstructure tensor of
0013 %     the image, see PCV (13.93)
0014 % This is controlled by the variable 'cov_im'
0015 
0016 % The uncertainty of the epipoles (if visible) and
0017 % the epipolar lines is visualized.
0018 %
0019 % The 8 corresponding points and the fundamental matrix may be stored,
0020 % and read in later. As soon as 8 correspoinding points are measured,
0021 % the old file will be owerwritten.
0022 %
0023 % 1. Choose two corresponding images
0024 %    (They are assumed to be in the folder ../Data/images)
0025 %    Six image pairs are prepared, so only the number is to be chosen
0026 %    (see the parameter 'image_pair')
0027 % 2. Decide on whether you want to measure 8 corresponidng points or
0028 %    read data from file (see the parameter 'readX')
0029 %    (they are lying in the folder /Data/ObservedImageCoordinates
0030 % 3. If chosen: measure the corresponding points
0031 %    one after the other in the first and the second image
0032 % 4. Measure individual points in the first image (now number 3)
0033 %    to get the epipolar line in the second image (number 4) visualized.
0034 %    In addition the expected precision of the measured
0035 %    parallaxes/disparities px' and py'is shown in image 3,
0036 %    and the standard deviations of px' and py' are give on the console.
0037 % 5. Step 4 can be repeated until right mouse button is pressed.
0038 %
0039 % The folder '../Data/images' contains the original images.
0040 % The folder '../Data/ObservedImageCoordinates' contains the measured points
0041 % for each image pair twice
0042 %    - as 'XXX.mat'-file which can be overwritten by measuring the
0043 %      8 corresponding points
0044 %    - as 'XXX-original.mat'-file, which contains already measured points,
0045 %      which can be used to recover from wrongly measured points.
0046 %
0047 % Wolfgang Förstner 4/2018
0048 % wfoerstn@uni-bonn.de
0049 %
0050 % last changes: Susanne Wenzel 05/18
0051 % wenzel@igg.uni-bonn.de
0052 
0053 function demo_epipolar_geometry()
0054 
0055 addpath(genpath('../../General-Functions'))
0056 addpath('../Functions')
0057 addpath(genpath('../Data'))
0058 
0059 ss = plot_init();
0060 close all
0061 
0062 % magnification = magnification factor for plotting
0063 plot_params.magnification_m = 300; % factor for measured point
0064 plot_params.magnification   =  3; % factor for epipolar geometry
0065 
0066 %% Control parameters
0067 % choose number of image pair
0068 image_pair = 1;
0069 
0070 % choose whether to measure 8 image pairs or read data from file
0071 % if readX = 0 then measure image correspondencies
0072 %          = 1 read them from file
0073 readX = 1;
0074 
0075 % choose what uncertainty is shown for the expected parallaxes/disparities
0076 % cov_im = 0 the uncertainty is the rounding error
0077 %        ~=0 the uncertainty is derived from the local structure tensor
0078 cov_im = 1;
0079 
0080 %%
0081 switch image_pair
0082     case 1
0083         %% Cathedral baseline in viewing direction
0084         Image_name_l = 'IMG_9473.JPG';
0085         Image_name_r = 'IMG_9474.JPG';
0086         plot_params.magnification_m = 200; % factor for measured point
0087         plot_params.magnification = 4;   % factor for epipolar geometry
0088     case 2
0089         %% Church
0090         Image_name_l = '008.JPG';
0091         Image_name_r = '009.JPG';
0092         plot_params.magnification_m = 100; % factor for measured point
0093         plot_params.magnification = 1;  % factor for epipolar geometry
0094     case 3
0095         %% Upside down building
0096         Image_name_l = 'IMG_1936.JPG';
0097         Image_name_r = 'IMG_1937.JPG';
0098         plot_params.magnification_m = 300; % factor for measured point
0099         plot_params.magnification = 6;   % factor for epipolar geometry
0100     case 4
0101         %% Cathedral b/w
0102         Image_name_l = 'IMG_9473-bw.JPG';
0103         Image_name_r = 'IMG_9474-bw.JPG';
0104         plot_params.magnification_m = 200; % factor for measured point
0105         plot_params.magnification = 4;   % factor for epipolar geometry
0106     case 5
0107         %% corn rigth left
0108         Image_name_l = '1336-34.jpg';
0109         Image_name_r = '1336-95.jpg';
0110         plot_params.magnification_m = 200; % factor for measured point
0111         plot_params.magnification = 2;   % factor for epipolar geometry
0112     case 6
0113         %% corn up/down
0114         Image_name_l=  '1337-1.jpg';
0115         Image_name_r = '1337-2.jpg';
0116         plot_params.magnification_m = 200; % factor for measured point
0117         plot_params.magnification = 2;   % factor for epipolar geometry
0118     otherwise
0119         disp('wrong number of image pair')          
0120 end
0121 
0122 %% default parameters
0123 
0124 % linewidths for plotting rectangles
0125 plot_params.f1 = 9;    % black (background)
0126 plot_params.f2 = 4;    % yellow (foreground)
0127 
0128 im_path = '../Data/images/';
0129 pt_path = '../Data/ObservedImageCoordinates/';
0130 
0131 % control parameters for structure tensor
0132 structure_tensor_params.tau = 0.7;   % differentiation scale
0133 structure_tensor_params.sigma = 3;   % integration scale for structure tensor
0134 structure_tensor_params.integration_width = ceil(1.5*structure_tensor_params.sigma)*2+1;
0135 structure_tensor_params.sigma_n = 0.03; % noise standard deviation
0136 
0137 %% initiate SUGR
0138 sugr_INIT
0139 
0140 %% read images
0141 Image_l = imread(fullfile(im_path,Image_name_l));
0142 Image_r = imread(fullfile(im_path,Image_name_r));
0143 coordinate_name = fullfile(pt_path,Image_name_l(1:end-4));
0144 
0145 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0146 %% show images
0147 f1 = figure('name','Left image');
0148 imshow(Image_l);
0149 set(f1,'Position',[0 ss(2)/2.4 ss(1)/2 ss(2)/2]);
0150 
0151 f2 = figure('name','Right image');
0152 imshow(Image_r);
0153 set(f2,'Position',[ss(1)/2 ss(2)/2.4 ss(1)/2 ss(2)/2]);
0154 
0155 if readX == 0   % measure image correspondencies
0156     
0157     [x1, x2] = measure_homologeous_Points(f1, f2, Image_l, Image_r, 8);
0158     
0159     if x1 == -1
0160         error('Not enough points. We need at least 8 homologeous points!')
0161         return                                                             %#ok<UNRCH>
0162     end
0163 
0164     X = [x1', x2'];
0165     
0166     save(coordinate_name,'X');
0167     
0168 else  % read measured image correspondencies from file
0169     
0170     load(coordinate_name,'X');
0171     
0172     N = size(X,1);                                                          %#ok<NODEF>
0173     
0174     figure(f1)
0175     hold on    
0176     for n = 1:N
0177         plot_square_with_background(X(n,1),X(n,2),100,plot_params.f1,plot_params.f2);
0178     end
0179     
0180     figure(f2)
0181     hold on    
0182     for n = 1:N
0183         plot_square_with_background(X(n,3),X(n,4),100,plot_params.f1,plot_params.f2);
0184     end
0185 end
0186 
0187 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0188 %% estimate F
0189 
0190 % Stochastic modell of observed image points
0191 Cxx = [ ones(8,1),zeros(8,4),ones(8,1),...
0192         zeros(8,4),ones(8,1),zeros(8,4),ones(8,1) ]...
0193         /12;
0194 
0195 % estimate F from measured images coordinates
0196 [F,Cff,el,Cll,er,Crr] = sugr_estimate_F_and_epipoles_algebraically(X,Cxx);
0197 %[F,Cff,el,er] = sugr_estimate_F_and_epipoles_algebraically_n(X,Cxx);
0198 
0199 %% start showing epipolar lines in Figure 4 for point in Figure 3
0200 
0201 %% plot epipole in left image
0202 f3 = figure('name','Left image');
0203 imshow(Image_l);
0204 set(f3,'Position',[0 0 ss(1)/2 ss(2)/2]);
0205 hold on
0206 plot_epipole(f3,el,Cll,plot_params)
0207 %% plot epipole in right image
0208 f4 = figure('name','Right image');
0209 imshow(Image_r);
0210 set(f4,'Position',[ss(1)/2 0 ss(1)/2 ss(2)/2]);
0211 hold on
0212 plot_epipole(f4,er,Crr,plot_params)
0213 
0214 
0215 while true 
0216     %% wait until a button is pressed
0217     disp('Measure point in Figure 3 to see epipolar line in Figure 4.')
0218     disp('Press right mouse button to exit')
0219     
0220 
0221     %% measure poin in left image
0222     x = measure_Point(f3, Image_l);
0223     
0224     if x == -1
0225         title('Interactively stopped')
0226         disp('Interactively stopped')
0227         break
0228     end
0229     
0230     [im_window, success] = get_image_window(Image_l, x, structure_tensor_params.integration_width);
0231     
0232     if success==1
0233         
0234         if cov_im == 0
0235             Cov = eye(2)/12;
0236             std_xy = sqrt(diag(Cov)');
0237         else
0238             T   = structure_tensor(im_window,...
0239                 structure_tensor_params.tau,...
0240                 structure_tensor_params.sigma);
0241             
0242             Cov = structure_tensor_params.sigma_n^2 * ...
0243                 inv( T * structure_tensor_params.integration_width^2 ...
0244                 + eye(2)*10^(-6) );                                        %#ok<MINV>
0245             std_xy = sqrt(diag(Cov)');
0246         end
0247         disp(['standard  deviations xy in pixel: ', num2str(std_xy)])        
0248         
0249         %plot epipole into left image again
0250         plot_epipole(f3,el,Cll,plot_params)        
0251                 
0252         %% plot ellipse in left image
0253         plot_ellipse_mc(x, Cov*plot_params.magnification_m^2,'-k',plot_params.f1);
0254         plot_ellipse_mc(x, Cov*plot_params.magnification_m^2,'-y',plot_params.f2);
0255         %hold off
0256 
0257         %% refresh image right
0258         figure(f4)
0259         imshow(Image_r);
0260         hold on
0261         
0262         
0263         % epipolar line
0264         lss =  F'*[x;1];
0265         
0266         % covariance matrix of lss due to uncertainty of F and point
0267         J = kron(eye(3),[x;1]');
0268         Cllh = J*Cff*J' + F'*[Cov,zeros(2,1);zeros(1,3)]*F;
0269         
0270         % plot epipolar line into right image
0271         epi = sugr_Line_2D(lss,Cllh);
0272         sugr_plot_Line_2D(epi,'-k','-k',plot_params.f1,plot_params.magnification);
0273         sugr_plot_Line_2D(epi,'-y','-y',plot_params.f2,plot_params.magnification);
0274         
0275         % plot epipole into right image again
0276         plot_epipole(f4,er,Crr,plot_params)
0277         %sugr_plot_Point_2D(er,'.k','-m',3,plot_params.magnification);
0278         
0279         hold off     
0280 
0281     end
0282 end
0283 
0284 return
0285 
0286 
0287 function plot_epipole(f,e,Cee,plot_params)
0288 
0289 figure(f)
0290 
0291 if abs(e(3)) > 0.0001
0292     
0293     % Plot standard ellipse of epipole
0294     ers=sugr_Point_2D(e,Cee);
0295     
0296     sugr_plot_Point_2D(ers,'.k','-w',plot_params.f1,plot_params.magnification)
0297     sugr_plot_Point_2D(ers,'.k','-m',plot_params.f2,plot_params.magnification)
0298     
0299     % plot text
0300     [e_E,C_EE]=sugr_get_Euclidean_Point_2D(ers);   
0301     text(e_E(1)+ 2*sqrt(C_EE(1,1))*plot_params.magnification,...
0302          e_E(2)+ 2*sqrt(C_EE(2,2))*plot_params.magnification,...
0303         'Epipole','FontSize',14,'FontWeight','bold','Color','m')
0304     
0305 end