Home > Matching_SYM_LSM > src > 2D_LSM_62 > 2D_LSM_62_sym > LSM_62_sym_estimate_parameters.m

LSM_62_sym_estimate_parameters

PURPOSE ^

LSM 62 symmetric: estimate parameters

SYNOPSIS ^

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)

DESCRIPTION ^

 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

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

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

Generated on Sun 19-Jul-2020 23:00:25 by m2html © 2005