LSM 62 symmetric: estimate parameters function [est_x,d_x,est_sigma_0,NN,Red,Xm,v] = ... LSM_62_sym_estimate_parameters(fe,fes,g,h,yi,zi,w,xa,pt) fe = NfxNf estimated image fes = NfxNf smoothed image g,h = NgxNg, NhxNh observations of image yi,zi = Kg, Kh list of observed pixels w = Kg+Kh vector of weights xa = actual estimates for geometric and radiometric transformation pt boolean for plot/print est_x = updated estimate (.x = 8-vector, .C = 8x8 CovM d_x = last correction of est_x.x est_sigma_0 = sqrt ofsqrt of estimated variance factor NN = number of observing pixels Red = redundancy Xm = design matrix v = residuals wf 11/2018
0001 % LSM 62 symmetric: estimate parameters 0002 % 0003 % function [est_x,d_x,est_sigma_0,NN,Red,Xm,v] = ... 0004 % LSM_62_sym_estimate_parameters(fe,fes,g,h,yi,zi,w,xa,pt) 0005 % 0006 % fe = NfxNf estimated image 0007 % fes = NfxNf smoothed image 0008 % g,h = NgxNg, NhxNh observations of image 0009 % yi,zi = Kg, Kh list of observed pixels 0010 % w = Kg+Kh vector of weights 0011 % xa = actual estimates for geometric and radiometric transformation 0012 % pt boolean for plot/print 0013 % 0014 % est_x = updated estimate (.x = 8-vector, .C = 8x8 CovM 0015 % d_x = last correction of est_x.x 0016 % est_sigma_0 = sqrt ofsqrt of estimated variance factor 0017 % NN = number of observing pixels 0018 % Red = redundancy 0019 % Xm = design matrix 0020 % v = residuals 0021 % 0022 % wf 11/2018 0023 0024 function [est_x,d_x,est_sigma_0,NN,Red,Xm,v] = ... 0025 LSM_62_sym_estimate_parameters(fe,fes,g,h,yi,zi,w,xa,pt) 0026 0027 %% affinities 0028 B = [xa(1) xa(3) xa(5);xa(2) xa(4) xa(6);0 0 1]; 0029 Bi = inv(B); 0030 S = [xa(7);xa(8)]; 0031 fak = S(1); 0032 faki = 1/fak; 0033 0034 %% half size of images 0035 Nf = (size(fe,1)-1)/2; 0036 Ng = (size(g,1)-1)/2; 0037 Nh = (size(h,1)-1)/2; 0038 0039 %% number of observations 0040 Kg = size(yi,2); 0041 Kh = size(zi,2); 0042 0043 % 0044 NN = Kg+Kh; 0045 NU = 8; 0046 % N - (8 geometric parameters + interior of fe) 0047 % interior = round(sqrt(Kg*Kh)) 0048 Red = NN-(NU+round(sqrt(Kg*Kh))); 0049 0050 % 0051 dg = zeros(NN,1); 0052 Xm = zeros(NN,8); 0053 0054 %% first image 0055 for k = 1:Kg 0056 i=yi(1,k)-Ng; % coordinates in g (centred) 0057 j=yi(2,k)-Ng; 0058 xh = B*[i;j;1]; % real coordinates in f (centred) 0059 xc=xh(1); 0060 yc=xh(2); 0061 x=xc+Nf+1; % pixel coordinates in f 0062 y=yc+Nf+1; 0063 % 0064 fij = LSM_f_cubic_interpolation(x,y,fe); % estimated signal 0065 fxij = LSM_fx_cubic_interpolation(x,y,fes); % gradient at signal 0066 fyij = LSM_fy_cubic_interpolation(x,y,fes); % consistent with fij 0067 % fxij = LSM_fx_cubic_smoothing(x,y,fes); % gradient at smoothed signal 0068 % fyij = LSM_fy_cubic_smoothing(x,y,fes); % see parallel implementation 0069 dg(k) = g(i+Ng+1,j+Ng+1)/255 - (fij-S(2))/S(1); % g: 255->1 0070 % Jacobian 0071 Xm(k,:) = [fxij*i fyij*i fxij*j fyij*j fxij fyij -(fij-S(2))*faki -1]*faki; 0072 0073 % Check Jacobian numerically, are ok 0074 % [my,Syy,J]=var_prop_classical(@f,mx,Sxx,p) 0075 % if k < 10 0076 % [fij_c,~,Jg] = var_prop_classical(@gj,xa,0.0000000001*ones(8),[i;j;fe(:)]) 0077 % (fij-S(2))/S(1) 0078 % Xm(k,:) 0079 % [fij_c-(fij-S(2))/S(1),Jg-Xm(k,:)] 0080 % end 0081 end 0082 0083 %% second image 0084 for k=1:Kh; 0085 i=zi(1,k)-Nh; % pixel coordinates in h (centred) 0086 j=zi(2,k)-Nh; 0087 xh = Bi*[i;j;1]; % real coordinates in f (centred) 0088 xc=xh(1); 0089 yc=xh(2); 0090 x=xc+Nf+1; % pixel coordinates in f 0091 y=yc+Nf+1; 0092 % 0093 fij = LSM_f_cubic_interpolation(x,y,fe); % estimated signal 0094 fxij = LSM_fx_cubic_interpolation(x,y,fes); % gradient at signal 0095 fyij = LSM_fy_cubic_interpolation(x,y,fes); % consistent with fij 0096 % fxij = LSM_fx_cubic_smoothing(x,y,fes); % gradient at smoothed signal 0097 % fyij = LSM_fy_cubic_smoothing(x,y,fes); % see parallel implementation 0098 pij = Bi(1:2,1:2)'*[fxij;fyij]; % gradient in f 0099 pxij = pij(1); 0100 pyij = pij(2); 0101 % residuals 0102 dg(Kg+k) = h(i+Nh+1,j+Nh+1)/255 - (S(1)*fij+S(2)); % h:255->1 0103 % Jacobians 0104 Xm(Kg+k,:) = -[pxij*xc pyij*xc pxij*yc pyij*yc pxij pyij -fij*faki -faki]*fak; 0105 0106 % Check Jacobian numerically, are ok 0107 % [my,Syy,J]=var_prop_classical(@f,mx,Sxx,p) 0108 % if k < 10 0109 % [fij_c,~,Jh] = var_prop_classical(@hk,xa,0.0000000001*ones(8),[i;j;fe(:)]) 0110 % (S(1)*fij+S(2)) 0111 % Xm(Kg+k,:) 0112 % [fij_c-(S(1)*fij+S(2)),Jh-Xm(Kg+k,:)] 0113 % end 0114 end 0115 0116 0117 %% Normal equation system 0118 N = Xm' * (Xm .* repmat(w,1,NU)); 0119 n = Xm' * (dg .* w); 0120 0121 %% Inverse and solution 0122 C_xx = inv(N); 0123 d_x = C_xx * n; 0124 0125 est_x.x = xa + d_x; 0126 est_x.C = C_xx / 255^2; % since weights refer to [0..255] 0127 0128 %% residuals 0129 v = Xm * d_x - dg; 0130 Omega = v' * (v.*w) * 255^2; % since weights refer to [0..255] 0131 est_sigma_0 = sqrt(Omega/Red); 0132 0133 0134 %% plot/print 0135 0136 if pt > 1 0137 N = N; 0138 n = n'; 0139 dx = d_x' 0140 est_x_ = est_x.x' 0141 sigma_0_est = est_sigma_0 0142 Covariance_2_sigma_rho(est_x.C); 0143 sigmas = diag(sqrt(est_x.C))'; 0144 B_est = reshape(est_x.x(1:6),2,3) 0145 Sigmas_B = reshape(sigmas(1:6),2,3) 0146 S_est = est_x.x(7:8)' 0147 Sigmas_S = sigmas(7:8) 0148 0149 figure('Name','Histograms of residuals') 0150 subplot(1,2,1) 0151 hist(v(1:Kg),ceil(sqrt(Kg))) 0152 title({['v\_g: mean = ',num2str(mean(v(1:Kg)))],... 0153 ['v\_g: std = ',num2str(std(v(1:Kg)))]}) 0154 subplot(1,2,2) 0155 hist(v(Kg+1:Kg+Kh),ceil(sqrt(Kh))) 0156 title({['v\_h: mean = ',num2str(mean(v(Kg+1:Kg+Kh)))],... 0157 ['v\_h: std = ',num2str(std(v(Kg+1:Kg+Kh)))]}) 0158 end 0159 0160 return