Changeset 264


Ignore:
Timestamp:
Oct 21, 2011, 6:32:03 AM (9 years ago)
Author:
sommeria
Message:

stereo option repaired

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/src/RUN_STLIN.m

    r38 r264  
    44function RUN_STLIN(file_A,file_B,vel_type,file_st,nx_patch,ny_patch,thresh_patch,fileAxml,fileBxml)
    55               
    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
    88 npxA=[]; npyA=[]; pxB=[]; npyB=[];
    99 if isfield(XmlDataA,'Camera') && isfield(XmlDataB,'Camera')
     
    2525     end
    2626 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
    2735 if isempty(npxA) ||isempty(npxB)
    2836     msgbox_uvmat('ERROR','The size of image A needs to be defined in the xml file ImaDoc')
     
    4553 end
    4654 
    47  %corners of each image in real coordinates:
     55 %corners of each image in px coordinates:
    4856 cornerA(:,1)=[0 0 npxA npxA]';%x positions
    4957 cornerA(:,2)=[0 npyA 0 npyA]';%y positions
    5058 cornerB(:,1)=[0 0 npxB npxB]';%x positions
    5159 cornerB(:,2)=[0 npyB 0 npyB]';%y positions
     60  %corners of each image in phys coordinates:
    5261[xyA(:,1),xyA(:,2)]=phys_XYZ(tsaiA,cornerA(:,1),cornerA(:,2));
    5362[xyB(:,1),xyB(:,2)]=phys_XYZ(tsaiB,cornerB(:,1),cornerB(:,2));
     
    7887% %read the velocity fields
    7988% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    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});
    8389%read field A
    8490[Field,VelTypeOut]=read_civxdata(file_A,[],vel_type);
    85 %interpolate on XimaA
     91%removes false vectors
     92if 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));
     97end
     98%interpolate on the grid common to both images in phys coordinates
     99dXa= griddata_uvmat(Field.X,Field.Y,Field.U,XimaA,YimaA);
     100dYa= griddata_uvmat(Field.X,Field.Y,Field.V,XimaA,YimaA);
     101dt=Field.dt;
     102time=Field.Time;
     103
     104%read field B
     105[Field,VelTypeOut]=read_civxdata(file_B,[],vel_type);
     106if ~isequal(Field.dt,dt)
     107    msgbox_uvmat('ERROR','different time intervals for the two velocity fields ')
     108     return
     109end
     110if ~isequal(Field.Time,time)
     111    msgbox_uvmat('ERROR','different times for the two velocity fields ')
     112     return
     113end
     114%removes false vectors
     115if isfield(Field,'FF')
    86116Field.X=Field.X(find(Field.FF==0));
    87117Field.Y=Field.Y(find(Field.FF==0));
    88118Field.U=Field.U(find(Field.FF==0));
    89119Field.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 B
    96 % [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      return
    102 end
    103 if ~isequal(Field.Time,time)
    104     msgbox_uvmat('ERROR','different times for the two velocity fields ')
    105      return
    106120end
    107121%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));
    112122dXb=griddata_uvmat(Field.X,Field.Y,Field.U,XimaB,YimaB);
    113123dYb=griddata_uvmat(Field.X,Field.Y,Field.V,XimaB,YimaB);
     
    122132 
    123133%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    124 %compute the coefficients
     134%compute the differential coefficients of the geometric calibration
    125135%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    126 
    127 
    128136[A11,A12,A13,A21,A22,A23]=pxcm_tsai(tsaiA,grid_phys1);
    129137[B11,B12,B13,B21,B22,B23]=pxcm_tsai(tsaiB,grid_phys1);
     
    175183% cte.title='grid';
    176184Result.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.VarDimIndex: {[1]  [1]  [1]  [1]  [1]  [1]}
     185% Result.ListDimName={'nb_vectors'}
     186% Result.DimValue=length(U);
     187Result.ListVarName={'vec_X','vec_Y','vec_U','vec_V','vec_W','vec_E'};
     188Result.VarDimName={'nb_vectors','nb_vectors','nb_vectors','nb_vectors','nb_vectors','nb_vectors'}
    181189Result.vec_X= grid_phys1(ind_error,1);
    182190Result.vec_Y= grid_phys1(ind_error,2);
     
    192200
    193201%'pxcm_tsai': find differentials of the Tsai calibration
    194 %
    195202function [A11,A12,A13,A21,A22,A23]=pxcm_tsai(a,var_phys)
    196 a_read=a;
    197 
    198203R=(a.R)';
    199204
     
    221226a.C2y=R(6)*R(8)-R(9)*R(5);
    222227
    223 
    224 %dependence in x,y
    225 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*y-a.C1x*z+R(1)*a.Tz-R(7)*a.Tx)./denom)/a.dpx;
    227 A12=(a.f*a.sx*(a.C12*x-a.C1y*z+R(2)*a.Tz-R(8)*a.Tx)./denom)/a.dpx;
    228 A21=(a.f*a.sx*(a.C21*y-a.C2x*z+R(4)*a.Tz-R(7)*a.Ty)./denom)/a.dpy;
    229 A22=(a.f*(a.C22*x-a.C2y*z+R(5)*a.Tz-R(8)*a.Ty)./denom)/a.dpy;
    230 A13=(a.f*(a.C1x*x+a.C1y*y+R(3)*a.Tz-R(9)*a.Tx)./denom)/a.dpx;
    231 A23=(a.f*(a.C2x*x+a.C2y*y+R(6)*a.Tz-R(9)*a.Ty)./denom)/a.dpy;
    232 
    233 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    234 %Old Version for z=0
    235 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    236 % %'camera' coordinates
    237 % 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 coordinates
    241 % Xu=a.f*xc./zc;
    242 % Yu=a.f*yc./zc;
    243 % %distorted image coordinates
    244 % distortion=(a.kappa1)*(Xu.*Xu+Yu.*Yu)+1; %!! intégrer derivation kappa
    245 % % distortion=1;
    246 % Xd=Xu./distortion;
    247 % Yd=Yu./distortion;
    248 % %pixel coordinates
    249 % X=Xd*a.sx/a.dpx+a.Cx;
    250 % Y=Yd/a.dpy+a.Cy;
    251 %
    252 % %transform coeff for differentiels
    253 % 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 %
    263228% %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.Tz-R(7)*a.Tx)./denom)/a.dpx;
    266 % A12=(a.f*a.sx*(a.C12*x+R(2)*a.Tz-R(8)*a.Tx)./denom)/a.dpx;
    267 % A21=(a.f*a.sx*(a.C21*y+R(4)*a.Tz-R(7)*a.Ty)./denom)/a.dpy;
    268 % A22=(a.f*(a.C22*x+R(5)*a.Tz-R(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*y-a.C1x*z+R(1)*a.Tz-R(7)*a.Tx)./denom)/a.dpx;
     231% A12=(a.f*a.sx*(a.C12*x-a.C1y*z+R(2)*a.Tz-R(8)*a.Tx)./denom)/a.dpx;
     232% A21=(a.f*a.sx*(a.C21*y-a.C2x*z+R(4)*a.Tz-R(7)*a.Ty)./denom)/a.dpy;
     233% A22=(a.f*(a.C22*x-a.C2y*z+R(5)*a.Tz-R(8)*a.Ty)./denom)/a.dpy;
    269234% A13=(a.f*(a.C1x*x+a.C1y*y+R(3)*a.Tz-R(9)*a.Tx)./denom)/a.dpx;
    270235% A23=(a.f*(a.C2x*x+a.C2y*y+R(6)*a.Tz-R(9)*a.Ty)./denom)/a.dpy;
    271 %
    272 
    273 
    274 
    275 
    276 
    277 
    278 
    279 
    280 
    281 
     236
     237%dependence in x,y
     238denom=(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));
     239A11=(a.fx_fy(1)*(a.C11*y-a.C1x*z+R(1)*a.Tx_Ty_Tz(3)-R(7)*a.Tx_Ty_Tz(1))./denom);
     240A12=(a.fx_fy(1)*(a.C12*x-a.C1y*z+R(2)*a.Tx_Ty_Tz(3)-R(8)*a.Tx_Ty_Tz(1))./denom);
     241A21=(a.fx_fy(1)*(a.C21*y-a.C2x*z+R(4)*a.Tx_Ty_Tz(3)-R(7)*a.Tx_Ty_Tz(2))./denom);
     242A22=(a.fx_fy(2)*(a.C22*x-a.C2y*z+R(5)*a.Tx_Ty_Tz(3)-R(8)*a.Tx_Ty_Tz(2))./denom);
     243A13=(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);
     244A23=(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.