Changeset 264
 Timestamp:
 Oct 21, 2011, 6:32:03 AM (9 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

trunk/src/RUN_STLIN.m
r38 r264 4 4 function RUN_STLIN(file_A,file_B,vel_type,file_st,nx_patch,ny_patch,thresh_patch,fileAxml,fileBxml) 5 5 6 [XmlDataA,error]=imadoc2struct(fileAxml); 7 [XmlDataB,error]=imadoc2struct(fileBxml); 6 [XmlDataA,error]=imadoc2struct(fileAxml);%read xml file associated to image series A 7 [XmlDataB,error]=imadoc2struct(fileBxml);%read xml file associated to image series B 8 8 npxA=[]; npyA=[]; pxB=[]; npyB=[]; 9 9 if isfield(XmlDataA,'Camera') && isfield(XmlDataB,'Camera') … … 25 25 end 26 26 end 27 %%%%%%%% added for Duran 28 if isfield(XmlDataA,'Npx')&&isfield(XmlDataA,'Npy')&&isfield(XmlDataB,'Npx')&&isfield(XmlDataB,'Npy') 29 npxA=XmlDataA.Npx; 30 npyA=XmlDataA.Npy; 31 npxB=XmlDataB.Npx; 32 npyB=XmlDataB.Npy; 33 end 34 %%%%%%%% added for Duran 27 35 if isempty(npxA) isempty(npxB) 28 36 msgbox_uvmat('ERROR','The size of image A needs to be defined in the xml file ImaDoc') … … 45 53 end 46 54 47 %corners of each image in realcoordinates:55 %corners of each image in px coordinates: 48 56 cornerA(:,1)=[0 0 npxA npxA]';%x positions 49 57 cornerA(:,2)=[0 npyA 0 npyA]';%y positions 50 58 cornerB(:,1)=[0 0 npxB npxB]';%x positions 51 59 cornerB(:,2)=[0 npyB 0 npyB]';%y positions 60 %corners of each image in phys coordinates: 52 61 [xyA(:,1),xyA(:,2)]=phys_XYZ(tsaiA,cornerA(:,1),cornerA(:,2)); 53 62 [xyB(:,1),xyB(:,2)]=phys_XYZ(tsaiB,cornerB(:,1),cornerB(:,2)); … … 78 87 % %read the velocity fields 79 88 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 80 %81 % [dt,time1,pixcmx,pixcmy,vec_X,vec_Y,vec_Z,vec_U,vec_V,vec_W,vec_C,vec_F,fixflag,vel_type_out,error,nb_coord,nb_dim]...82 % =read_vel({filecell_ncA},{vel_type});83 89 %read field A 84 90 [Field,VelTypeOut]=read_civxdata(file_A,[],vel_type); 85 %interpolate on XimaA 91 %removes false vectors 92 if isfield(Field,'FF') 93 Field.X=Field.X(find(Field.FF==0)); 94 Field.Y=Field.Y(find(Field.FF==0)); 95 Field.U=Field.U(find(Field.FF==0)); 96 Field.V=Field.V(find(Field.FF==0)); 97 end 98 %interpolate on the grid common to both images in phys coordinates 99 dXa= griddata_uvmat(Field.X,Field.Y,Field.U,XimaA,YimaA); 100 dYa= griddata_uvmat(Field.X,Field.Y,Field.V,XimaA,YimaA); 101 dt=Field.dt; 102 time=Field.Time; 103 104 %read field B 105 [Field,VelTypeOut]=read_civxdata(file_B,[],vel_type); 106 if ~isequal(Field.dt,dt) 107 msgbox_uvmat('ERROR','different time intervals for the two velocity fields ') 108 return 109 end 110 if ~isequal(Field.Time,time) 111 msgbox_uvmat('ERROR','different times for the two velocity fields ') 112 return 113 end 114 %removes false vectors 115 if isfield(Field,'FF') 86 116 Field.X=Field.X(find(Field.FF==0)); 87 117 Field.Y=Field.Y(find(Field.FF==0)); 88 118 Field.U=Field.U(find(Field.FF==0)); 89 119 Field.V=Field.V(find(Field.FF==0)); 90 dXa= griddata_uvmat(Field.X,Field.Y,Field.U,XimaA,YimaA);91 dYa= griddata_uvmat(Field.X,Field.Y,Field.V,XimaA,YimaA);92 dt=Field.dt;93 time=Field.Time;94 95 %read field B96 % [dt,time2,pixcmx,pixcmy,vec_X,vec_Y,vec_Z,vec_U,vec_V,vec_W,vec_C,vec_F,fixflag,vel_type_out,error,nb_coord,nb_dim]...97 % =read_vel({file_B},{vel_type});98 [Field,VelTypeOut]=read_civxdata(file_B,FieldNames,vel_type);99 if ~isequal(Field.dt,dt)100 msgbox_uvmat('ERROR','different time intervals for the two velocity fields ')101 return102 end103 if ~isequal(Field.Time,time)104 msgbox_uvmat('ERROR','different times for the two velocity fields ')105 return106 120 end 107 121 %interpolate on XimaB 108 Field.X=Field.X(find(Field.FF==0));109 Field.Y=Field.Y(find(Field.FF==0));110 Field.U=Field.U(find(Field.FF==0));111 Field.V=Field.V(find(Field.FF==0));112 122 dXb=griddata_uvmat(Field.X,Field.Y,Field.U,XimaB,YimaB); 113 123 dYb=griddata_uvmat(Field.X,Field.Y,Field.V,XimaB,YimaB); … … 122 132 123 133 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 124 %compute the coefficients134 %compute the differential coefficients of the geometric calibration 125 135 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 126 127 128 136 [A11,A12,A13,A21,A22,A23]=pxcm_tsai(tsaiA,grid_phys1); 129 137 [B11,B12,B13,B21,B22,B23]=pxcm_tsai(tsaiB,grid_phys1); … … 175 183 % cte.title='grid'; 176 184 Result.civ=0;%not a civ file (no direct correspondance with an image) 177 Result.ListDimName={'nb_vectors'}178 Result.DimValue=length(U);179 Result.ListVarName={'vec_X' ;'vec_Y';'vec_U';'vec_V';'vec_W';'vec_E'};180 Result.VarDim Index: {[1] [1] [1] [1] [1] [1]}185 % Result.ListDimName={'nb_vectors'} 186 % Result.DimValue=length(U); 187 Result.ListVarName={'vec_X','vec_Y','vec_U','vec_V','vec_W','vec_E'}; 188 Result.VarDimName={'nb_vectors','nb_vectors','nb_vectors','nb_vectors','nb_vectors','nb_vectors'} 181 189 Result.vec_X= grid_phys1(ind_error,1); 182 190 Result.vec_Y= grid_phys1(ind_error,2); … … 192 200 193 201 %'pxcm_tsai': find differentials of the Tsai calibration 194 %195 202 function [A11,A12,A13,A21,A22,A23]=pxcm_tsai(a,var_phys) 196 a_read=a;197 198 203 R=(a.R)'; 199 204 … … 221 226 a.C2y=R(6)*R(8)R(9)*R(5); 222 227 223 224 %dependence in x,y225 denom=(R(7)*x+R(8)*y+R(9)*z+a.Tz).*(R(7)*x+R(8)*y+R(9)*z+a.Tz);226 A11=(a.f*a.sx*(a.C11*ya.C1x*z+R(1)*a.TzR(7)*a.Tx)./denom)/a.dpx;227 A12=(a.f*a.sx*(a.C12*xa.C1y*z+R(2)*a.TzR(8)*a.Tx)./denom)/a.dpx;228 A21=(a.f*a.sx*(a.C21*ya.C2x*z+R(4)*a.TzR(7)*a.Ty)./denom)/a.dpy;229 A22=(a.f*(a.C22*xa.C2y*z+R(5)*a.TzR(8)*a.Ty)./denom)/a.dpy;230 A13=(a.f*(a.C1x*x+a.C1y*y+R(3)*a.TzR(9)*a.Tx)./denom)/a.dpx;231 A23=(a.f*(a.C2x*x+a.C2y*y+R(6)*a.TzR(9)*a.Ty)./denom)/a.dpy;232 233 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%234 %Old Version for z=0235 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%236 % %'camera' coordinates237 % xc=R(1)*x+R(2)*y+a.Tx;238 % yc=R(4)*x+R(5)*y+a.Ty;239 % zc=R(7)*x+R(8)*y+a.Tz;240 % %undistorted image coordinates241 % Xu=a.f*xc./zc;242 % Yu=a.f*yc./zc;243 % %distorted image coordinates244 % distortion=(a.kappa1)*(Xu.*Xu+Yu.*Yu)+1; %!! intégrer derivation kappa245 % % distortion=1;246 % Xd=Xu./distortion;247 % Yd=Yu./distortion;248 % %pixel coordinates249 % X=Xd*a.sx/a.dpx+a.Cx;250 % Y=Yd/a.dpy+a.Cy;251 %252 % %transform coeff for differentiels253 % a.C11=R(1)*R(8)R(2)*R(7);254 % a.C12=R(2)*R(7)R(1)*R(8);255 % a.C21=R(4)*R(8)R(5)*R(7);256 % a.C22=R(5)*R(7)R(4)*R(8);257 % a.C1x=R(3)*R(7)R(9)*R(1);258 % a.C1y=R(3)*R(8)R(9)*R(2);259 % a.C2x=R(6)*R(7)R(9)*R(4);260 % a.C2y=R(6)*R(8)R(9)*R(5);261 %262 %263 228 % %dependence in x,y 264 % denom=(R(7)*x+R(8)*y+ a.Tz).*(R(7)*x+R(8)*y+a.Tz);265 % A11=(a.f*a.sx*(a.C11*y +R(1)*a.TzR(7)*a.Tx)./denom)/a.dpx;266 % A12=(a.f*a.sx*(a.C12*x +R(2)*a.TzR(8)*a.Tx)./denom)/a.dpx;267 % A21=(a.f*a.sx*(a.C21*y +R(4)*a.TzR(7)*a.Ty)./denom)/a.dpy;268 % A22=(a.f*(a.C22*x +R(5)*a.TzR(8)*a.Ty)./denom)/a.dpy;229 % denom=(R(7)*x+R(8)*y+R(9)*z+a.Tz).*(R(7)*x+R(8)*y+R(9)*z+a.Tz); 230 % A11=(a.f*a.sx*(a.C11*ya.C1x*z+R(1)*a.TzR(7)*a.Tx)./denom)/a.dpx; 231 % A12=(a.f*a.sx*(a.C12*xa.C1y*z+R(2)*a.TzR(8)*a.Tx)./denom)/a.dpx; 232 % A21=(a.f*a.sx*(a.C21*ya.C2x*z+R(4)*a.TzR(7)*a.Ty)./denom)/a.dpy; 233 % A22=(a.f*(a.C22*xa.C2y*z+R(5)*a.TzR(8)*a.Ty)./denom)/a.dpy; 269 234 % A13=(a.f*(a.C1x*x+a.C1y*y+R(3)*a.TzR(9)*a.Tx)./denom)/a.dpx; 270 235 % A23=(a.f*(a.C2x*x+a.C2y*y+R(6)*a.TzR(9)*a.Ty)./denom)/a.dpy; 271 % 272 273 274 275 276 277 278 279 280 281 236 237 %dependence in x,y 238 denom=(R(7)*x+R(8)*y+R(9)*z+a.Tx_Ty_Tz(3)).*(R(7)*x+R(8)*y+R(9)*z+a.Tx_Ty_Tz(3)); 239 A11=(a.fx_fy(1)*(a.C11*ya.C1x*z+R(1)*a.Tx_Ty_Tz(3)R(7)*a.Tx_Ty_Tz(1))./denom); 240 A12=(a.fx_fy(1)*(a.C12*xa.C1y*z+R(2)*a.Tx_Ty_Tz(3)R(8)*a.Tx_Ty_Tz(1))./denom); 241 A21=(a.fx_fy(1)*(a.C21*ya.C2x*z+R(4)*a.Tx_Ty_Tz(3)R(7)*a.Tx_Ty_Tz(2))./denom); 242 A22=(a.fx_fy(2)*(a.C22*xa.C2y*z+R(5)*a.Tx_Ty_Tz(3)R(8)*a.Tx_Ty_Tz(2))./denom); 243 A13=(a.fx_fy(2)*(a.C1x*x+a.C1y*y+R(3)*a.Tx_Ty_Tz(3)R(9)*a.Tx_Ty_Tz(1))./denom); 244 A23=(a.fx_fy(2)*(a.C2x*x+a.C2y*y+R(6)*a.Tx_Ty_Tz(3)R(9)*a.Tx_Ty_Tz(2))./denom); 245 246 247 248 249 250 251 252
Note: See TracChangeset
for help on using the changeset viewer.