% LS-matching 6+2 of two images g and h, symmetric model: observed images g and h referring to unknown reference f g(y) = a_7^-1 * (f(B * y + a) - a_8) + n_g(y) h(z) = a_7 * (f(B^-1 * (z - a) + a_8) + n_h(z) B_a = [u1 u3 u5; S = [u7 u8 u2 u4 u6; 0 1] 0 0 1] geometric transformation A = B^2 radiometric transformation R = S^2 all images are square: size (M x M), M odd radius N=(M-1)/2, centre at (N+1,N+1) call: [est_x,est_sigma_0,Red,N_iter] = LSM_62_sym_main(g,h,vg,vh,A_a,R_a,sigma_smooth,max_iter,pt) g = Ng x Ng matrix of first image h = Nh x Nh matrix of second image vg = 256-vector of noise variances [0:1:255] of g vh = 256-vector of noise variances [0:1:255] of h A_a = 3x3 matrix, approximate transformation geometric parameters A_a = [a b c; d e f; 0 0 1] zh = A_a * yh R_a = 1x2 vector of approximate contrast and intensity deviation h = R_a(1)* g + R_a(2) sigma_smooth = smoothing parameter for true image at first iteration max_iter = maximum number of iterations Nf_MIN = minimum size of overlap (radius of window f) pt in (0,1,2) plot output control est_x = estimated transformation parameters est_x.x = parameters est_x.C = covariance matrix 1..6: geometric parameters 7..8: intenisty parameters est_sigma_0 estimated factor sigma_0 (should be close to 1) NN number of observing pixels Mf size of estimated image f Red redundancy N_iter number of iterations needed max_ratio maximum ratio of corrections of parameters/sigma wf 7/2020
0001 %% LS-matching 6+2 of two images g and h, symmetric 0002 % 0003 % model: observed images g and h referring to unknown reference f 0004 % g(y) = a_7^-1 * (f(B * y + a) - a_8) + n_g(y) 0005 % h(z) = a_7 * (f(B^-1 * (z - a) + a_8) + n_h(z) 0006 % 0007 % B_a = [u1 u3 u5; S = [u7 u8 0008 % u2 u4 u6; 0 1] 0009 % 0 0 1] 0010 % geometric transformation A = B^2 0011 % radiometric transformation R = S^2 0012 % 0013 % all images are square: 0014 % size (M x M), M odd 0015 % radius N=(M-1)/2, 0016 % centre at (N+1,N+1) 0017 % 0018 % call: 0019 % [est_x,est_sigma_0,Red,N_iter] = 0020 % LSM_62_sym_main(g,h,vg,vh,A_a,R_a,sigma_smooth,max_iter,pt) 0021 % 0022 % g = Ng x Ng matrix of first image 0023 % h = Nh x Nh matrix of second image 0024 % vg = 256-vector of noise variances [0:1:255] of g 0025 % vh = 256-vector of noise variances [0:1:255] of h 0026 % A_a = 3x3 matrix, approximate transformation geometric parameters 0027 % A_a = [a b c; 0028 % d e f; 0029 % 0 0 1] 0030 % zh = A_a * yh 0031 % R_a = 1x2 vector of approximate contrast and intensity deviation 0032 % h = R_a(1)* g + R_a(2) 0033 % sigma_smooth = smoothing parameter for true image at first iteration 0034 % max_iter = maximum number of iterations 0035 % Nf_MIN = minimum size of overlap (radius of window f) 0036 % pt in (0,1,2) plot output control 0037 % 0038 % est_x = estimated transformation parameters 0039 % est_x.x = parameters 0040 % est_x.C = covariance matrix 0041 % 1..6: geometric parameters 0042 % 7..8: intenisty parameters 0043 % est_sigma_0 estimated factor sigma_0 (should be close to 1) 0044 % NN number of observing pixels 0045 % Mf size of estimated image f 0046 % Red redundancy 0047 % N_iter number of iterations needed 0048 % max_ratio maximum ratio of corrections of parameters/sigma 0049 % 0050 % wf 7/2020 0051 0052 function [est_x,est_sigma_0,NN,Mf,Red,N_iter,max_ratio] =... 0053 LSM_62_sym_main(g,h,vg,vh,A_a,R_a,sigma_smooth,max_iter,Nf_MIN,pt) 0054 0055 0056 % check for postiveness 0057 if det(A_a) < 0 0058 disp('Geometric affinity is not positive definite') 0059 est_x = 0; 0060 est_sigma_0 = 0; 0061 Red = -1; 0062 NN = 0; 0063 Mf = 0; 0064 N_iter = iter; 0065 return 0066 end 0067 0068 %% test CPU-time 0069 NO_TEST_CPU = 1; 0070 sum_cpu=0; 0071 0072 %% determine approximate values: take squareroot of given affinities 0073 0074 % take squareroot from given affinities 0075 B_a = sqrtm(A_a); 0076 a7 = sqrt(R_a(1)); 0077 a8 = R_a(2)/(1+a7); 0078 S_a = [a7,a8]'; 0079 0080 % set approximate values 0081 Be = B_a(1:2,:); 0082 xa = [Be(:);S_a]; 0083 0084 if pt > 0 0085 [im_lr,weite,border] = fill_left_right_images_0(g,h,max_iter); 0086 end 0087 0088 if pt > 1 0089 fig_conv = show_main_convergence_0(xa); 0090 end 0091 0092 %% start iterations ============================================= 0093 est_sigma_0 = 1; 0094 max_ratio = 10; 0095 0096 for iter=1:max_iter 0097 if pt > 0 0098 disp(['Iteration No.: ',num2str(iter), ' -----------------------------']) 0099 end 0100 0101 %% find observation positions 0102 if pt > 1 0103 tic; 0104 end 0105 if max_ratio > 0.5 0106 [Nf,yi,zi,w] = LSM_62_sym_find_observation_positions(... 0107 g,h,vg,vh,xa,Nf_MIN,NO_TEST_CPU*pt); 0108 end 0109 if Nf < Nf_MIN 0110 disp('Overlap is too small') 0111 est_x = 0; 0112 est_sigma_0 = 0; 0113 Red = -1; 0114 NN = 0; 0115 Mf = Nf; 0116 N_iter = iter; 0117 max_ratio = 100; 0118 return 0119 end 0120 0121 %% determine estimated fe and smoothed fes 0122 0123 0124 [fe,fes,fg,fh] = LSM_62_sym_estimate_f(... 0125 g,h,vg,vh,Nf,xa,sigma_smooth,NO_TEST_CPU*pt); 0126 if pt > 1 0127 cpu=toc; 0128 sum_cpu=sum_cpu+cpu; 0129 end 0130 % adapt smoothing parameter 0131 sigma_smooth=sigma_smooth*0.5; 0132 0133 if pt > 0 0134 im_lr = fill_left_right_images(... 0135 im_lr,fg,fh,weite,iter); 0136 end 0137 0138 %% estimate parameters 0139 old_est_sigma_0=est_sigma_0; 0140 if pt > 1 0141 old_cpu=cpu; 0142 tic; 0143 end 0144 0145 % ------------- main routine -------------------- 0146 [est_xa,d_xa,est_sigma_0,NN,Red,Xm,dg] = ... 0147 LSM_62_sym_estimate_parameters(... 0148 fe,fes,g,h,yi,zi,w,xa,NO_TEST_CPU*pt); 0149 % ------------- main routine -------------------- 0150 0151 if pt > 1 0152 cpu=cpu+toc; 0153 old_cpu=cpu; 0154 sum_cpu=sum_cpu+cpu; 0155 end 0156 if pt > 1 0157 show_main_convergence(... 0158 fig_conv,xa,est_xa,old_est_sigma_0,est_sigma_0,iter,old_cpu,cpu,sum_cpu); 0159 end 0160 xa = est_xa.x; 0161 0162 % 0163 %% close iteration 0164 0165 if pt > 1 0166 check=(Xm'*diag(w)*dg); 0167 check_AtWv=check' 0168 end 0169 max_ratio = max(abs(d_xa./sqrt(diag(est_xa.C)))); 0170 %% plot/print 0171 if pt > 0 0172 display(['N = ',num2str(NN),', sigma_0 = ', num2str(est_sigma_0),', ratio = ',num2str(max_ratio)]) 0173 end 0174 0175 if max_ratio < 0.1 0176 break 0177 end 0178 0179 0180 end 0181 0182 N_iter=iter ; 0183 Mf = size(fe,1); 0184 0185 0186 if pt > 0 0187 figure('name','undistorted images') 0188 imshow(im_lr) 0189 0190 for iter = 1:max_iter 0191 text(border+(iter+0.5)*weite,1.75*border+weite,num2str(iter)) 0192 end 0193 end 0194 %% Finale 0195 0196 est_x = finalize_sym_lsm(... 0197 iter,... 0198 xa,... 0199 est_xa,... 0200 pt); 0201 0202 0203 return