0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
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
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
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
0042 Kg = size(yi,2);
0043 Kh = size(zi,2);
0044
0045
0046 NN = Kg+Kh;
0047 NU = 8;
0048
0049
0050 Red = NN-(NU+round(sqrt(Kg*Kh)));
0051
0052
0053 dg = zeros(NN,1);
0054 Xm = zeros(NN,8);
0055
0056 Derx = [3 10 3; 0 0 0; -3 -10 -3]/ 32;
0057 fesx = conv2(fes,Derx ,'same');
0058 fesy = conv2(fes,Derx','same');
0059
0060 fsxy(:,:,1)=fe';
0061 fsxy(:,:,2)=fesx';
0062 fsxy(:,:,3)=fesy';
0063
0064
0065
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
0068 BCi = inv(BC);
0069
0070 TBC = maketform('affine',BCi');
0071
0072 xdata = [1 size(g,1)]; ydata = [1 size(g,2)];
0073
0074 fsxygo = imtransform(fsxy,TBC,'bicubic','xdata',xdata,'ydata',ydata);
0075
0076
0077 if size(fsxygo,1) < size(g,1)
0078
0079 Mfo = size(fsxygo,1);
0080 Nfo = (Mfo-1)/2;
0081
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
0091 feh = fsxyg(:,:,1)';
0092 fesxg = fsxyg(:,:,2)';
0093 fesyg = fsxyg(:,:,3)';
0094
0095 feij_v = feh(yi(1,:)'+Mg*(yi(2,:)'-1));
0096
0097 fxij_v = fesxg(yi(1,:)'+Mg*(yi(2,:)'-1));
0098
0099 fyij_v = fesyg(yi(1,:)'+Mg*(yi(2,:)'-1));
0100
0101 dg_v = g(yi(1,:)'+Mg*(yi(2,:)'-1))/255 - (feij_v-S(2))/S(1);
0102
0103 ij_v = [yi(1,:)',yi(2,:)',ones(Kg,1)] * ([1 0 -Ng-1;0 1 -Ng-1;0 0 1])';
0104
0105 xc_v = [yi(1,:)',yi(2,:)',ones(Kg,1)] * (B * [1 0 -Ng-1;0 1 -Ng-1;0 0 1])';
0106
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
0113
0114
0115
0116
0117
0118
0119
0120
0121
0122
0123
0124
0125
0126
0127
0128
0129
0130
0131
0132
0133
0134
0135
0136
0137
0138
0139
0140
0141
0142
0143
0144
0145
0146
0147
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
0150 BCh = inv(BChi);
0151
0152 TBCh = maketform('affine',BCh');
0153
0154 xdata = [1 size(h,1)]; ydata = [1 size(h,2)];
0155
0156 fsxyho = imtransform(fsxy,TBCh,'bicubic','xdata',xdata,'ydata',ydata);
0157
0158
0159 if size(fsxyho,1) < size(h,1)
0160
0161 Mfo = size(fsxyho,1);
0162 Nfo = (Mfo-1)/2;
0163
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
0173 feh = fsxyh(:,:,1)';
0174 fesxh = fsxyh(:,:,2)';
0175 fesyh = fsxyh(:,:,3)';
0176
0177
0178 feij_v = feh(zi(1,:)'+Mh*(zi(2,:)'-1));
0179
0180 fxij_v = fesxh(zi(1,:)'+Mh*(zi(2,:)'-1));
0181
0182 fyij_v = fesyh(zi(1,:)'+Mh*(zi(2,:)'-1));
0183
0184 dg_v = [dg_v;...
0185 h(zi(1,:)'+Mh*(zi(2,:)'-1))/255 - (S(1)*feij_v+S(2))];
0186
0187 xc_v = [zi(1,:)',zi(2,:)',ones(Kh,1)] * (Bi * [1 0 -Nh-1;0 1 -Nh-1;0 0 1])';
0188
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
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
0201
0202
0203
0204
0205
0206
0207
0208
0209
0210
0211
0212
0213
0214
0215
0216
0217
0218
0219
0220
0221
0222
0223
0224
0225
0226
0227
0228
0229
0230
0231
0232
0233
0234
0235
0236
0237
0238
0239 Xm = Xm_v;
0240 dg = dg_v;
0241
0242
0243 N = Xm' * (Xm .* repmat(w,1,NU));
0244 n = Xm' * (dg .* w);
0245
0246
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;
0252
0253
0254 v = Xm * d_x - dg;
0255 Omega = v' * (v.*w) * 255^2;
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