% Check correctness of Sym-LSM The demo routine demo_LSM_simulated.m is meant to check the correctness of the implementation based on simulated data. It allows to monitor the individual iterations for a single case, or to statistically test, whether the resultant parameters and their covariance matrix are coherent with the theoretical values, the true values of the parameters and the theore2018tical covariance matrix derived by the estimation procedure (Cramer-Rao bound). symmetric model =============== g(z) = 1/a_7 * (f(B * z + a ) - a_8) + m_h(z) h(z) = a_7 * f(Bi * z + ai) + a_8 + n_h(z) B=[a1 a3 a5] S = [a7 a8] [a2 a4 a6] [0 1 ] [0 0 1 ] A = B^2 R = S^2 main parameters =============== img_filename = name of image file for generating true image 'corner': corner image, simulated 'filename': read from file s,t = parameters for noise sigma_n^2=s+t*f Tg = geometric parameters 2x3 matrix [a1 a3 a5] [a2 a4 a6] Tr = radiometric parameters [a_7,a_8]' Ng = size of first image Nh = size of second image sigma = smoothing kernel size, if = 0: no smoothing init_rand = number of random generator, 0 --> arbitrary N_samples = number of samples generated (for checking covariance matrix) N_iter = number of iterations wf 11/2018, 01/2020, 7/2020
0001 %% Check correctness of Sym-LSM 0002 % 0003 % The demo routine demo_LSM_simulated.m is meant to check the 0004 % correctness of the implementation based on simulated data. It 0005 % allows to monitor the individual iterations for a single case, 0006 % or to statistically test, whether the resultant parameters and 0007 % their covariance matrix are coherent with the theoretical values, 0008 % the true values of the parameters and the theore2018tical covariance 0009 % matrix derived by the estimation procedure (Cramer-Rao bound). 0010 % 0011 % symmetric model 0012 % =============== 0013 % g(z) = 1/a_7 * (f(B * z + a ) - a_8) + m_h(z) 0014 % h(z) = a_7 * f(Bi * z + ai) + a_8 + n_h(z) 0015 % B=[a1 a3 a5] S = [a7 a8] 0016 % [a2 a4 a6] [0 1 ] 0017 % [0 0 1 ] 0018 % A = B^2 R = S^2 0019 % 0020 % main parameters 0021 % =============== 0022 % img_filename = name of image file for generating true image 0023 % 'corner': corner image, simulated 0024 % 'filename': read from file 0025 % s,t = parameters for noise sigma_n^2=s+t*f 0026 % Tg = geometric parameters 2x3 matrix 0027 % [a1 a3 a5] 0028 % [a2 a4 a6] 0029 % Tr = radiometric parameters [a_7,a_8]' 0030 % Ng = size of first image 0031 % Nh = size of second image 0032 % sigma = smoothing kernel size, if = 0: no smoothing 0033 % init_rand = number of random generator, 0 --> arbitrary 0034 % N_samples = number of samples generated 0035 % (for checking covariance matrix) 0036 % N_iter = number of iterations 0037 % 0038 % wf 11/2018, 01/2020, 7/2020 0039 0040 clearvars 0041 close all 0042 %clc 0043 0044 disp(' --------------------------------------------------------------------- ') 0045 disp(' ---- symmetric LSM with 6 geometric and 2 radiometric parameters ---- ') 0046 disp(' --------------------------------------------------------------------- ') 0047 0048 addpath(genpath('src/')) 0049 0050 0051 %% === Choose your own parameters ======================================= 0052 % --- initial random number (default = 0) ------------------------------- 0053 par.init_rand = 2; 0054 % init_rand = 0; 0055 0056 % --- Number of samples (default = 50) ---------------------------------- 0057 par.N_samples = 50; 0058 % par.N_samples = 1; 0059 0060 % --- Type of data, in [0..4], default = 4 ------------------------------ 0061 % par.type_data = 1; % 1 square 0062 % par.type_data = 2; % 2 waves 0063 % par.type_data = 3; % 3 MIT-image, subwindow 0064 par.type_data = 4; % 4 smoothed noise 0065 0066 % --- test symmetry, only has an effect if N_samples=1 ----------------- 0067 par.test_symmetry=0; 0068 %par.test_symmetry=1; 0069 0070 % --- use warping routine --- 0071 par.warping_routine =1 ; 0072 % par.warping_routine =0 ; 0073 0074 %% === Set and display main parameters ================================== 0075 par = simulated_set_parameters(par); 0076 0077 0078 simulated_display_parameters(par); 0079 0080 par = simulated_prepare_true_values(par); 0081 0082 %% generate reference image ============================================ 0083 0084 [I_true,f_true,g_true,h_true,var] = ... 0085 generate_reference_image(par); 0086 0087 %% start samples, may be only one ====================================== 0088 0089 % set statistics 0090 sigma_0s = zeros(par.N_samples,1); 0091 est_xs = zeros(par.N_samples,8); 0092 N_iters = zeros(par.N_samples,1); 0093 NNs = zeros(par.N_samples,1); 0094 init_randoms = zeros(par.N_samples,1); 0095 0096 % for measuring CPU time 0097 lsm_time=0; 0098 r_sample = 0; 0099 for sample = 1:par.N_samples 0100 0101 monitor_samples(sample,par.N_samples); 0102 0103 % --- Generate f_true, g, and h from I_true ------------------------ 0104 st = randn('state'); 0105 par.init_randoms(sample)= st(2); 0106 0107 [g,h,std_i_g,std_i_h] = LSM_62_sym_generate_noisy_observing_images... 0108 (g_true,h_true,par); 0109 0110 if sample == 1 0111 figure('name','first sample image pair') 0112 subplot(1,2,1) 0113 imshow(g/255); 0114 title('noisy left image') 0115 subplot(1,2,2) 0116 imshow(h/255); 0117 title('noisy right image') 0118 end 0119 0120 % --- adapt noise: take smoothing of signal into account ----------- 0121 vgm=par.vg + 1/sqrt(par.R_a(1)) * std_i_g^2*0.38; 0122 vhm=par.vh + sqrt(par.R_a(1)) * std_i_h^2*0.38; 0123 0124 0125 % ---- estimate transformation (LSM_62_main) ----------------------- 0126 % set parameters 0127 A_a = par.A_a; % 3x3 matrix, approx. geom. transf. 0128 R_a = par.R_a; % 1x2 vector, approx. rad. transf. 0129 sigma_smooth = par.sigma_smooth; % smoothing parameter 0130 max_iter = par.max_iter; % maximal number of iterations 0131 Nf_MIN = par.Nf_MIN; % minimum size of overlap 0132 pt = par.plot_type; % output level 0133 0134 tic; 0135 0136 if par.warping_routine 0137 [est_x,est_sigma_0,NN,Nf,Red,N_iter] = ... 0138 LSM_62_sym_warp_main(g,h,vgm,vhm,A_a,R_a,sigma_smooth,max_iter,Nf_MIN,pt); 0139 else 0140 [est_x,est_sigma_0,NN,Nf,Red,N_iter] = ... 0141 LSM_62_sym_main(g,h,vgm,vhm,A_a,R_a,sigma_smooth,max_iter,Nf_MIN,pt); 0142 end 0143 0144 0145 est_xi = est_x; 0146 if Red > 0 && N_iter < par.max_iter 0147 % successful matching 0148 r_sample = r_sample+1; 0149 0150 lsm_time=lsm_time+toc; 0151 0152 sigma_0s(r_sample) = est_sigma_0; 0153 est_xs(r_sample,:) = est_x.x'; 0154 N_iters(r_sample) = N_iter; 0155 NNs(r_sample) = NN; 0156 0157 % If N_samples = 1 and test_symmetry = 1 0158 if par.N_samples == 1 && par.test_symmetry 0159 disp('--- start symmetric estimation h->g ---') 0160 % set parameters 0161 Ai_a = inv(par.A_a); 0162 Ri_a = [1/par.R_a(1),-par.R_a(2)/par.R_a(1)]; 0163 % estimate 0164 [est_xi,est_sigma_0i,NN,Nf,Redi,N_iter] = ... 0165 LSM_62_sym_warp_main... 0166 (h,g,vhm,vgm,Ai_a,Ri_a,sigma_smooth,max_iter,Nf_MIN,pt); 0167 0168 if par.plot_type > 1 0169 Ai_est = [reshape(est_xi.x(1:6),2,3);0 0 1] 0170 Ri_est = [est_xi.x(7:8)';0 1] 0171 end 0172 end 0173 end 0174 0175 end 0176 0177 0178 %% analyse 0179 0180 par.plot_type = 1; 0181 0182 sym_lsm_analysis(... 0183 r_sample,... 0184 est_sigma_0,... 0185 est_x,... 0186 est_xi,... 0187 est_xs, sigma_0s, N_iters,... 0188 Red,... 0189 lsm_time,... 0190 par); 0191 0192 0193 disp(' --------------------------------- ') 0194 disp(' ---- end demo symmetric LSM ---- ') 0195 disp(' --------------------------------- ') 0196 0197 0198 % file_name = ['LSM_sym_par_pyr_result_',num2str(init_rand_prolog),'_',num2str(init_rand_samples),'_',num2str(type_data,'%02.f'),'.mat'] 0199 % save(file_name) 0200 0201 0202 return