Home > Matching_SYM_LSM > src > 2D_LSM_62 > 2D_LSM_62_sym_warp > LSM_62_sym_warp_estimate_parameters.m

LSM_62_sym_warp_estimate_parameters

PURPOSE ^

% LSM 62 symmetric: estimate parameters

SYNOPSIS ^

function [est_x,d_x,est_sigma_0,NN,Red,Xm,v] =LSM_62_sym_warp_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_warp_estimate_parameters(fe,fes,g,h,yi,zi,w,xa,pt)

 fe          = MfxMf estimated image
 fes         = MfxMf smoothed image
 g,h         = MgxMg, MhxMh 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 of sqrt of estimated variance factor
 NN          = number of observing pixels
 Red         = redundancy
 Xm          = design matrix
 v           = residuals

 wf 11/2018, 1/2020

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_warp_estimate_parameters(fe,fes,g,h,yi,zi,w,xa,pt)
0005 %
0006 % fe          = MfxMf estimated image
0007 % fes         = MfxMf smoothed image
0008 % g,h         = MgxMg, MhxMh 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 of sqrt 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, 1/2020
0023 
0024 function [est_x,d_x,est_sigma_0,NN,Red,Xm,v]  = ...
0025     LSM_62_sym_warp_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 Mg = size(g,1);
0037 Mh = size(h,1);
0038 Ng = (Mg-1)/2;
0039 Nh = (Mh-1)/2;
0040 
0041 %% number of observations
0042 Kg   = size(yi,2);
0043 Kh   = size(zi,2);
0044 
0045 %% keynumbers of estimation
0046 NN  = Kg+Kh;
0047 NU  = 8;
0048 % N - (8 geometric parameters + interior of fe)
0049 % interior = round(sqrt(Kg*Kh))
0050 Red = NN-(NU+round(sqrt(Kg*Kh)));
0051 
0052 %
0053 dg = zeros(NN,1);
0054 Xm = zeros(NN,8);
0055 %% Prepare derivatives of f
0056 Derx = [3 10 3; 0 0 0; -3 -10 -3]/ 32; % Scharr's filter
0057 fesx = conv2(fes,Derx ,'same');
0058 fesy = conv2(fes,Derx','same');
0059 % concatenate for simultaneous warping
0060 fsxy(:,:,1)=fe';
0061 fsxy(:,:,2)=fesx';
0062 fsxy(:,:,3)=fesy';
0063 
0064 %% first image
0065 % backward  complete transform
0066 BC = [1 0 Nf+1;0 1 Nf+1;0 0 1] * B * [1 0 -Ng-1;0 1 -Ng-1;0 0 1];
0067 % forward complete transform (for maketform.m)
0068 BCi  = inv(BC);
0069 % provide transformation (transposed inverse)
0070 TBC = maketform('affine',BCi');
0071 % specify target window = g
0072 xdata  = [1 size(g,1)];  ydata = [1 size(g,2)];
0073 % simultanous warping of f, fx, fy
0074 fsxygo  = imtransform(fsxy,TBC,'bicubic','xdata',xdata,'ydata',ydata);
0075 
0076 % enforce same size if smaller than g
0077 if size(fsxygo,1) < size(g,1)
0078     % size of fsxygo
0079     Mfo = size(fsxygo,1);
0080     Nfo = (Mfo-1)/2;
0081     % possible shift
0082     shift = max(0,Ng-Nfo);
0083     T = [1 0 0; 0 1 0; shift shift 1];
0084     TT = maketform('affine',T);
0085     fsxyg = imtransform(fsxygo,TT,'XData',xdata,'YData',ydata);
0086 else
0087     fsxyg = fsxygo;
0088 end
0089 
0090 % extract
0091 feh   = fsxyg(:,:,1)';
0092 fesxg = fsxyg(:,:,2)';
0093 fesyg = fsxyg(:,:,3)';
0094 % estimated signal
0095 feij_v  = feh(yi(1,:)'+Mg*(yi(2,:)'-1));
0096 % fx of estimated signal
0097 fxij_v  = fesxg(yi(1,:)'+Mg*(yi(2,:)'-1));
0098 % fy of estimated signal
0099 fyij_v  = fesyg(yi(1,:)'+Mg*(yi(2,:)'-1));
0100 % residuals
0101 dg_v    = g(yi(1,:)'+Mg*(yi(2,:)'-1))/255 - (feij_v-S(2))/S(1);
0102 % transformed coordinates (Kg*3 matrix)
0103 ij_v    = [yi(1,:)',yi(2,:)',ones(Kg,1)] * ([1 0 -Ng-1;0 1 -Ng-1;0 0 1])';
0104 % transformed coordinates (Kg*3 matrix)
0105 xc_v    = [yi(1,:)',yi(2,:)',ones(Kg,1)] * (B * [1 0 -Ng-1;0 1 -Ng-1;0 0 1])';
0106 % Design matrix for additive correction of affinity
0107 Xm_v    = [ fxij_v.*ij_v(:,1)   fyij_v.*ij_v(:,1) ...
0108     fxij_v.*ij_v(:,2)   fyij_v.*ij_v(:,2) ...
0109     fxij_v              fyij_v ...
0110     -(feij_v-S(2))*faki -ones(Kg,1)]*faki;
0111 
0112 % % % sequential (old code)
0113 % profile=zeros(Kg,5);
0114 % for k = 1:Kg
0115 %     i=yi(1,k)-Ng-1;                 % coordinates in g (centred)
0116 %     j=yi(2,k)-Ng-1;
0117 %     xh = B*[i;j;1];                 % real coordinates in f (centred)
0118 %     xc=xh(1);
0119 %     yc=xh(2);
0120 %     x=xc+Nf+1;                      % pixel coordinates in f
0121 %     y=yc+Nf+1;
0122 %     fij  = LSM_f_cubic_interpolation(x,y,fe);    % estimated signal
0123 %     fxij = LSM_fx_cubic_interpolation(x,y,fes);  % gradient at smoothed signal
0124 %     fyij = LSM_fy_cubic_interpolation(x,y,fes);
0125 %
0126 %     dg(k)  = g(i+Ng+1,j+Ng+1)/255 - (fij-S(2))/S(1);    % g: 255->1
0127 %     % Jacobian (see report)
0128 %     Xm(k,:) = [fxij*i fyij*i fxij*j fyij*j fxij fyij -(fij-S(2))*faki -1]*faki;
0129 %
0130 %
0131 %     profile(k,:)=[xc,yc,fij,fxij,fyij];
0132 %     %     w(k) = 1/((1/w(k) - 1/12)*...
0133 %     %         factors_V_cubic_interpolation(x)*...
0134 %     %         factors_V_cubic_interpolation(y) +1/12);
0135 % %     % Check Jacobian numerically, are ok
0136 % %     % [my,Syy,J]=var_prop_classical(@f,mx,Sxx,p)
0137 % %         if k < 10
0138 % %             [fij_c,~,Jg] = var_prop_classical(@gj,xa,0.0000000001*ones(8),[i;j;fe(:)])
0139 % %             (fij-S(2))/S(1)
0140 % %             Xm(k,:)
0141 % %             [fij_c-(fij-S(2))/S(1),Jg-Xm(k,:)]
0142 % %         end
0143 % end
0144 
0145 
0146 %% second image
0147 % backward  complete transform
0148 BChi = [1 0 Nf+1;0 1 Nf+1;0 0 1] * Bi * [1 0 -Nh-1;0 1 -Nh-1;0 0 1];
0149 % forward complete transform (for maketform.m)
0150 BCh  = inv(BChi);
0151 % provide transformation (transposed)
0152 TBCh = maketform('affine',BCh');
0153 % specify target window = h
0154 xdata  = [1 size(h,1)];  ydata = [1 size(h,2)];
0155 % simultanous warping of f, fx, fy
0156 fsxyho  = imtransform(fsxy,TBCh,'bicubic','xdata',xdata,'ydata',ydata);
0157 
0158 % enforce same size
0159 if size(fsxyho,1) < size(h,1)
0160     % size of fsxyho
0161     Mfo = size(fsxyho,1);
0162     Nfo = (Mfo-1)/2;
0163     % possible shift
0164     shift = max(0,Nh-Nfo);
0165     T = [1 0 0; 0 1 0; shift shift 1];
0166     TT = maketform('affine',T);
0167     fsxyh = imtransform(fsxyho,TT,'XData',xdata,'YData',ydata);
0168 else
0169     fsxyh = fsxyho;
0170 end
0171 
0172 % extract the three arrays f,fx,fy
0173 feh   = fsxyh(:,:,1)';
0174 fesxh = fsxyh(:,:,2)';
0175 fesyh = fsxyh(:,:,3)';
0176 
0177 % estimated signal
0178 feij_v  = feh(zi(1,:)'+Mh*(zi(2,:)'-1));
0179 % fx of estimated signal
0180 fxij_v  = fesxh(zi(1,:)'+Mh*(zi(2,:)'-1));
0181 % fy of estimated signal
0182 fyij_v  = fesyh(zi(1,:)'+Mh*(zi(2,:)'-1));
0183 % residuals
0184 dg_v    = [dg_v;...
0185     h(zi(1,:)'+Mh*(zi(2,:)'-1))/255 - (S(1)*feij_v+S(2))];
0186 % transformed coordinates (Kh*3 matrix)
0187 xc_v    = [zi(1,:)',zi(2,:)',ones(Kh,1)] * (Bi * [1 0 -Nh-1;0 1 -Nh-1;0 0 1])';
0188 % transform gradients
0189 pij_v   = [fxij_v,fyij_v]*Bi(1:2,1:2);
0190 pxij_v  = pij_v(:,1);
0191 pyij_v  = pij_v(:,2);
0192 % Design matrix
0193 Xm_v    = [Xm_v;
0194     -[pxij_v.*xc_v(:,1) pyij_v.*xc_v(:,1) ...
0195     pxij_v.*xc_v(:,2) pyij_v.*xc_v(:,2) ...
0196     pxij_v            pyij_v ...
0197     -feij_v*faki      -ones(Kh,1)*faki]*fak];
0198 
0199 
0200 % sequential (old code)
0201 % profileh=zeros(Kh,7);
0202 % fijh = zeros(Mh);
0203 % for k=1:Kh;
0204 %     i=zi(1,k)-Nh-1;                   % pixel coordinates in h (centred)
0205 %     j=zi(2,k)-Nh-1;
0206 %     xh = Bi*[i;j;1];                % real coordinates in f (centred)
0207 %     xc=xh(1);
0208 %     yc=xh(2);
0209 %     x=xc+Nf+1;                      % pixel coordinates in f
0210 %     y=yc+Nf+1;
0211 %
0212 %     fij  = LSM_f_cubic_interpolation(x,y,fe);    % estimated signal
0213 %     fxij = LSM_fx_cubic_smoothing(x,y,fes);  % gradient at smoothed signal
0214 %     fyij = LSM_fy_cubic_smoothing(x,y,fes);
0215 %
0216 %     pij  = Bi(1:2,1:2)'*[fxij;fyij]; % gradient in f
0217 %     pxij = pij(1);
0218 %     pyij = pij(2);
0219 %     % residuals
0220 %     dg(Kg+k)   = h(i+Nh+1,j+Nh+1)/255 - (S(1)*fij+S(2));     % h:255->1
0221 %     % Jacobians  (see report)
0222 %     Xm(Kg+k,:) = -[pxij*xc pyij*xc pxij*yc pyij*yc pxij pyij -fij*faki -faki]*fak;
0223 %
0224 %     profileh(k,:)=[xc,yc,fij,fxij,fyij,pxij,pyij];
0225 %     fijh(zi(1,k),zi(2,k))= dg(Kg+k);
0226 %
0227 %
0228 %     % Check Jacobian numerically, are ok
0229 %     % [my,Syy,J]=var_prop_classical(@f,mx,Sxx,p)
0230 %     %     if k < 10
0231 %     %         [fij_c,~,Jh] = var_prop_classical(@hk,xa,0.0000000001*ones(8),[i;j;fe(:)])
0232 %     %         (S(1)*fij+S(2))
0233 %     %         Xm(Kg+k,:)
0234 %     %         [fij_c-(S(1)*fij+S(2)),Jh-Xm(Kg+k,:)]
0235 %     %     end
0236 % end
0237 
0238 %% take parallel results
0239 Xm = Xm_v;
0240 dg = dg_v;
0241 
0242 %% Normal equation system
0243 N  = Xm' * (Xm   .* repmat(w,1,NU));
0244 n  = Xm' * (dg .* w);
0245 
0246 %% Inverse and solution
0247 C_xx     = inv(N);
0248 d_x      = C_xx * n;
0249 
0250 est_x.x  = xa + d_x;
0251 est_x.C  = C_xx / 255^2;       % since weights refer to [0..255]
0252 
0253 %% residuals
0254 v     = Xm * d_x - dg;
0255 Omega = v' * (v.*w) * 255^2;   % since weights refer to [0..255]
0256 est_sigma_0 = sqrt(Omega/Red);
0257 
0258 
0259 if pt > 1
0260     N           = N;
0261     n           = n';
0262     dx          = d_x'
0263     est_x_      = est_x.x'
0264     sigma_0_est = est_sigma_0
0265     covariance_2_sigma_rho(est_x.C);
0266     sigmas      = diag(sqrt(est_x.C))';
0267     B_est       = reshape(est_x.x(1:6),2,3)
0268     Sigmas_B    = reshape(sigmas(1:6),2,3)
0269     S_est       = est_x.x(7:8)'
0270     Sigmas_S    = sigmas(7:8)
0271     
0272     figure(13)
0273     subplot(1,2,1)
0274     hist(v(1:Kg),ceil(sqrt(Kg)))
0275     title({['v\_g: mean = ',num2str(mean(v(1:Kg)))],...
0276         ['v\_g: std  = ',num2str(std(v(1:Kg)))]})
0277     subplot(1,2,2)
0278     hist(v(Kg+1:Kg+Kh),ceil(sqrt(Kh)))
0279     title({['v\_h: mean = ',num2str(mean(v(Kg+1:Kg+Kh)))],...
0280         ['v\_h: std  = ',num2str(std(v(Kg+1:Kg+Kh)))]})
0281     
0282 end
0283 
0284 return

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