Description of demo_estimate_2D_homography_from_point

0001 %% demo_estimate_2D_homography_from_point_pairs
0002 %
0003 % generate 3x3 grid of points
0004 % transform them with a homography
0005 % add noise, possibly correlated between correspoinding points
0006 %
0007 % estimate homography (algebraical and ML-estimation)
0008 %
0009 % check estimation result (algebraic and ML-estimation)
0010 %
0011 % show scatter plots for transformed points
0012 %      with theoretical covariance matrices
0013 %
0014 % see PCV Sect. 10.6.3, Fig. 10.28
0015 % choose correlation in lines 36 ff.
0016 %
0017 % Wolfgang Förstner 2015
0018 % wfoerstn@uni-bonn.de
0019 
0020 
0021 close all
0022 clearvars
0023 %clear all
0024 
0025 addpath(genpath('../../General-Functions'))
0026 
0027 global plot_option
0028 global print_option
0029 global print_option_estimation
0030 
0031 disp('=========================================')
0032 disp('      Estimate 2D homography             ')
0033 disp('-----------------------------------------')
0034 sugr_INIT
0035 ss = plot_init;
0036 
0037 %% control parameters
0038 
0039 % random seed
0040 init_rand      = 1;
0041 %rinit_rand      = 0;
0042 
0043 % number of cases
0044 %N_samples     = 300;
0045 %N_samples     = 100;
0046 N_samples     = 50;
0047 %N_samples     = 1;
0048 disp(['Number of samples                      : ', num2str(N_samples)])
0049 
0050 sigma_x       = 0.001;
0051 disp(['Standard deviation of directions       : ', num2str(sigma_x)])
0052 % correlation between points
0053 %rho           = -0.98;
0054 %rho           =     0;
0055 rho           = +0.98;
0056 
0057 disp(['Correlation of homologeous points      : ', num2str(rho)])
0058 
0059 % Threshold for iteration
0060 T       = 10^(-8);
0061 maxiter = 10;
0062 
0063 %significance level
0064 S   = 0.999;
0065 Slo = (1-S)/2;
0066 Sup = 1-Slo;
0067 
0068 %% random number intialization
0069 init_rand_seed(init_rand);
0070 
0071 %% Generate points in unit square
0072 sigma_x       = 0.001;                  % std. dev.
0073 sigma_y       = sigma_x;                  % std. dev.
0074 N             = 9;                      % number of lines
0075 H_tilde_e     = [1   -0.2  2.6;...
0076     0.1  1.2 -0.3; ...
0077     0.25 0.2    1];          % true homography matrix
0078 H_tilde       = sugr_Homography_2D(H_tilde_e);   % true homography
0079 factor_p      = 100;
0080 b_random      = 0;                    % random points in [-1,+1]^2
0081 
0082 
0083 print_option_estimation=1;
0084 if N_samples >1
0085     print_option_estimation=0;
0086 end
0087 % do not plot samples
0088 plot_option = 2;
0089 print_option = 2;
0090 if N_samples > 1
0091     plot_option = 0;
0092     print_option = 0;
0093 end
0094 
0095 
0096 % generate true points
0097 PPt = sugr_generate_true_2D_point_pairs_homography(H_tilde.H,N,b_random);
0098 
0099 
0100 est_H_samples_a    = zeros(N_samples,9);
0101 est_H_samples_ml   = zeros(N_samples,9);
0102 sigma_0s           = zeros(N_samples,1);
0103 figure('Color','w','Position',[20 ss(2)/3 ss(1)/3 ss(2)/2])
0104 hold on
0105 
0106 start=cputime;
0107 for i=1:N_samples
0108     if N_samples < 100
0109         fprintf(num2str(i)),fprintf(',')
0110         if mod(i,10)==0
0111             disp(' ')
0112         end
0113     else
0114         if mod(i,10)==0
0115             fprintf(num2str(i)),fprintf(',')
0116         end
0117         if mod(i,100)==0
0118             disp(' ')
0119         end
0120     end
0121     %disp(['sample: ',num2str(i)])
0122     % perturb point pairs
0123     PP = sugr_perturb_2D_point_pairs(PPt,sigma_x,sigma_y,rho,factor_p);
0124     
0125     %% Estimate point algebraically
0126     est_H_a = sugr_estimation_algebraic_Homography_2D_from_point_pairs(PP);
0127     
0128     %% Estimate point maximum likelihood
0129     [est_H_ml,sigma_0,R] = sugr_estimation_ml_Homography_2D_from_point_pairs(PP,est_H_a.H,T,maxiter);
0130     
0131     est_H_samples_a(i,:)    = est_H_a.H(:);
0132     est_H_samples_ml(i,:)   = est_H_ml.H(:);
0133     sigma_0s(i)             = sigma_0;
0134 end
0135 disp(['CPU time for ', num2str(N_samples),' samples               : ',num2str(cputime-start)]);
0136 
0137 
0138 %% evaluate result
0139 
0140 
0141 
0142 if N_samples > 2
0143     
0144     %% plot transferred points in 3x3 grid
0145     
0146     k=0;
0147     xp=zeros(N*N_samples,1);
0148     yp=zeros(N*N_samples,1);
0149     for i=1:N_samples
0150         for n=1:N
0151             %       given noisy points
0152             %[xn,yn,Cxy]   = sugr_select_Point_Pair_2D(PP,n);
0153             xn = PP.h(n,1:3)';
0154             yn = PP.h(n,4:6)';
0155             xne = xn(1:2)/xn(3);
0156             % true points as reference
0157             %[xt,yt,Cxy] = sugr_select_Point_Pair_2D(PPt,n);
0158             xt = PPt.h(n,1:3)';
0159             yt = PPt.h(n,4:6)';
0160             xte = xt(1:2)/xt(3);
0161             yte = yt(1:2)/yt(3);
0162             % true point transformed with estimated homography algebraic
0163             %y_a = sugr_transform_with_Homography_2D(est_H,xt);
0164             est_H = reshape(est_H_samples_a(i,:),3,3);
0165             y_a = est_H * xt;
0166             %ye_a = y_a.h(1:2)/y_a.h(3);
0167             ye_a = y_a(1:2)/y_a(3);
0168             %
0169             k=k+1;
0170             xp(k) = yte(1)+factor_p*(ye_a(1)-yte(1));
0171             yp(k) = yte(2)+factor_p*(ye_a(2)-yte(2));
0172             %plot(yte(1)+factor_p*(ye_a(1)-yte(1)),yte(2)+factor_p*(ye_a(2)-yte(2)),'.b');
0173         end
0174     end
0175     figure(1)
0176     hold on
0177     plot(xp,yp,'.b');
0178     
0179     
0180     k=0;
0181     xp=zeros(N*N_samples,1);
0182     yp=zeros(N*N_samples,1);
0183     for i=1:N_samples
0184         for n=1:N
0185             %       given noisy points
0186             %[xn,yn,Cxy]   = sugr_select_Point_Pair_2D(PP,n);
0187             xn = PP.h(n,1:3)';
0188             yn = PP.h(n,4:6)';
0189             xne = xn(1:2)/xn(3);
0190             % true points as reference
0191             %[xt,yt,Cxy] = sugr_select_Point_Pair_2D(PPt,n);
0192             xt = PPt.h(n,1:3)';
0193             yt = PPt.h(n,4:6)';
0194             xte = xt(1:2)/xt(3);
0195             yte = yt(1:2)/yt(3);
0196             % true point transformed with estimated homography algebraic
0197             %y_a = sugr_transform_with_Homography_2D(est_H,xt);
0198             est_H = reshape(est_H_samples_ml(i,:),3,3);
0199             y_a = est_H * xt;
0200             %ye_a = y_a.h(1:2)/y_a.h(3);
0201             ye_a = y_a(1:2)/y_a(3);
0202             %
0203             k=k+1;
0204             xp(k) = yte(1)+factor_p*(ye_a(1)-yte(1));
0205             yp(k) = yte(2)+factor_p*(ye_a(2)-yte(2));
0206             %plot(yte(1)+factor_p*(ye_a(1)-yte(1)),yte(2)+factor_p*(ye_a(2)-yte(2)),'.b');
0207         end
0208     end
0209     figure('Color','w','Position',[100+ss(1)/3  ss(2)/3 ss(1)/3 ss(2)/2])
0210     hold on
0211     plot(xp,yp,'.b');
0212     
0213     
0214     %
0215     est_H_mean_ml = mean(est_H_samples_ml);
0216     est_H_cov_ml  = cov(est_H_samples_ml);
0217     est_H_mean_a = mean(est_H_samples_a);
0218     est_H_cov_a  = cov(est_H_samples_a);
0219     
0220     est_H_ml_emp   = sugr_Homography_2D(reshape(est_H_mean_ml,3,3),est_H_cov_ml);
0221     est_H_a_emp   = sugr_Homography_2D(reshape(est_H_mean_a,3,3),est_H_cov_a);
0222     
0223     
0224     
0225     figure(1)
0226     %sugr_plot_Homography_2D(est_H_a_emp,'ok','-k',1,3*factor_p,'H');
0227     sugr_plot_Homography_2D(est_H_a,'ok','-k',1,3*factor_p,'H');
0228     title('ALG: samples with theoretical CovM')
0229     axis equal
0230     figure(2)
0231     sugr_plot_Homography_2D(est_H_ml,'ok','-k',1,3*factor_p,'H');
0232     title('ML: samples with theoretical CovM')
0233     axis equal
0234     
0235     R = 2*N-8;
0236     
0237     % for reducing the homography parameters (spectral normalization)
0238     Jhr = sugr_get_Jacobian_Jrh_Homography_2D(H_tilde.H);
0239     % Algebraic
0240     check_estimation_result(R,zeros(8,1),est_H_a.Crr,ones(N_samples,1),...
0241         est_H_samples_a*Jhr',S,' ALG')
0242     
0243     
0244     % Maximum likelihood
0245     check_estimation_result(R,zeros(8,1),est_H_ml.Crr,sigma_0s.^2,...
0246         est_H_samples_ml*Jhr',S,' ML');
0247     
0248 end
0249 
0250 % save(strcat('test_sugr_homography_estimation_a_ml_',num2str(rho*1000,'%4d')))
demo_estimate_2D_homography_from_point_pairs

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SOURCE CODE
