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