[8]  1  %'RUN_STLIN': combine velocity fields for stereo PIV


 2  % file_A,file_B: input velocity files


 3  %vel_type: string ='civ1' or 'civ2'


 4  function RUN_STLIN(file_A,file_B,vel_type,file_st,nx_patch,ny_patch,thresh_patch,fileAxml,fileBxml)


 5 


 6  [XmlDataA,error]=imadoc2struct(fileAxml);


 7  [XmlDataB,error]=imadoc2struct(fileBxml);


 8  npxA=[]; npyA=[]; pxB=[]; npyB=[];


 9  if isfield(XmlDataA,'Camera') && isfield(XmlDataB,'Camera')


 10  if isfield(XmlDataA.Camera,'ImageSize')&& isfield(XmlDataB.Camera,'ImageSize')


 11  ImageSizeA=XmlDataA.Camera.ImageSize;


 12  ImageSizeB=XmlDataB.Camera.ImageSize;


 13  if ~isempty(ImageSizeA)&& ~isempty(ImageSizeB)


 14  xindex=findstr(ImageSizeA,'x');


 15  if length(xindex)>=2


 16  npxA=str2num(ImageSizeA(1:xindex(1)1));


 17  npyA=str2num(ImageSizeA(xindex(1)+1:xindex(2)1));


 18  end


 19  xindex=findstr(ImageSizeB,'x');


 20  if length(xindex)>=2


 21  npxB=str2num(ImageSizeB(1:xindex(1)1));


 22  npyB=str2num(ImageSizeB(xindex(1)+1:xindex(2)1));


 23  end


 24  end


 25  end


 26  end


 27  if isempty(npxA) isempty(npxB)


[38]  28  msgbox_uvmat('ERROR','The size of image A needs to be defined in the xml file ImaDoc')


[8]  29  return


 30  elseif isempty(npxB)  isempty(npyB)


[38]  31  msgbox_uvmat('ERROR','The size of image B needs to be defined in the xml file ImaDoc')


[8]  32  return


 33  end


 34  if isfield(XmlDataA,'GeometryCalib')


 35  tsaiA=XmlDataA.GeometryCalib;


 36  else


[38]  37  msgbox_uvmat('ERROR','no geometric calibration available for image A')


[8]  38  return


 39  end


 40  if isfield(XmlDataB,'GeometryCalib')


 41  tsaiB=XmlDataB.GeometryCalib;


 42  else


[38]  43  msgbox_uvmat('ERROR','no geometric calibration available for image B')


[8]  44  return


 45  end


 46 


 47  %corners of each image in real coordinates:


 48  cornerA(:,1)=[0 0 npxA npxA]';%x positions


 49  cornerA(:,2)=[0 npyA 0 npyA]';%y positions


 50  cornerB(:,1)=[0 0 npxB npxB]';%x positions


 51  cornerB(:,2)=[0 npyB 0 npyB]';%y positions


 52  [xyA(:,1),xyA(:,2)]=phys_XYZ(tsaiA,cornerA(:,1),cornerA(:,2));


 53  [xyB(:,1),xyB(:,2)]=phys_XYZ(tsaiB,cornerB(:,1),cornerB(:,2));


 54  max_x=max(max(xyA(:,1)),max(xyB(:,1)));%maximum on the 4 corners of the the images


 55  min_x=min(min(xyA(:,1)),min(xyB(:,1)));%minimum on the 4 corners of the the images


 56  max_y=max(max(xyA(:,2)),max(xyB(:,2)));


 57  min_y=min(min(xyA(:,2)),min(xyB(:,2)));


 58  array_realx=[min_x:(max_xmin_x)/(nx_patch1):max_x];


 59  array_realy=[min_y:(max_ymin_y)/(ny_patch1):max_y];


 60  [grid_realx,grid_realy]=meshgrid(array_realx,array_realy);


 61  grid_real(:,1)=reshape(grid_realx,nx_patch*ny_patch,1);


 62  grid_real(:,2)=reshape(grid_realy,nx_patch*ny_patch,1);


 63  grid_real(:,3)=zeros(nx_patch*ny_patch,1);


 64  [grid_imaA(:,1),grid_imaA(:,2)]=px_XYZ(tsaiA,grid_real(:,1),grid_real(:,2));


 65  [grid_imaB(:,1),grid_imaB(:,2)]=px_XYZ(tsaiB,grid_real(:,1),grid_real(:,2));


 66 


 67  flagA=grid_imaA(:,1)>0 & grid_imaA(:,1)<npxA & grid_imaA(:,2)>0 & grid_imaA(:,2)<npyA;


 68  flagB=grid_imaB(:,1)>0 & grid_imaB(:,1)<npxB & grid_imaB(:,2)>0 & grid_imaB(:,2)<npyB;


 69  ind_good=find(flagA==1&flagB==1);


 70  XimaA=grid_imaA(ind_good,1);


 71  YimaA=grid_imaA(ind_good,2);


 72  XimaB=grid_imaB(ind_good,1);


 73  YimaB=grid_imaB(ind_good,2);


 74  grid_real_x=grid_real(ind_good,1);


 75  grid_real_y=grid_real(ind_good,2);


 76 


 77  % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


 78  % %read the velocity fields


 79  % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


 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  %read field A


 84  [Field,VelTypeOut]=read_civxdata(file_A,[],vel_type);


 85  %interpolate on XimaA


 86  Field.X=Field.X(find(Field.FF==0));


 87  Field.Y=Field.Y(find(Field.FF==0));


 88  Field.U=Field.U(find(Field.FF==0));


 89  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 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)


[38]  100  msgbox_uvmat('ERROR','different time intervals for the two velocity fields ')


[8]  101  return


 102  end


 103  if ~isequal(Field.Time,time)


[38]  104  msgbox_uvmat('ERROR','different times for the two velocity fields ')


[8]  105  return


 106  end


 107  %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  dXb=griddata_uvmat(Field.X,Field.Y,Field.U,XimaB,YimaB);


 113  dYb=griddata_uvmat(Field.X,Field.Y,Field.V,XimaB,YimaB);


 114  %eliminate NotaNumber


 115  ind_Nan=find(and(~isnan(dXa),~isnan(dXb)));


 116  dXa=dXa(ind_Nan);


 117  dYa=dYa(ind_Nan);


 118  dXb=dXb(ind_Nan);


 119  dYb=dYb(ind_Nan);


 120  grid_phys1(:,1)=grid_real_x(ind_Nan);


 121  grid_phys1(:,2)=grid_real_y(ind_Nan);


 122 


 123  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


 124  %compute the coefficients


 125  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


 126 


 127 


 128  [A11,A12,A13,A21,A22,A23]=pxcm_tsai(tsaiA,grid_phys1);


 129  [B11,B12,B13,B21,B22,B23]=pxcm_tsai(tsaiB,grid_phys1);


 130 


 131  C1=A11.*A22A12.*A21;


 132  C2=A13.*A22A12.*A23;


 133  C3=A13.*A21A11.*A23;


 134  D1=B11.*B22B12.*B21;


 135  D2=B13.*B22B12.*B23;


 136  D3=B13.*B21B11.*B23;


 137  A1=(A22.*D1.*(C1.*D3C3.*D1)+A21.*D1.*(C2.*D1C1.*D2));


 138  A2=(A12.*D1.*(C3.*D1C1.*D3)+A11.*D1.*(C1.*D2C2.*D1));


 139  B1=(B22.*C1.*(C3.*D1C1.*D3)+B21.*C1.*(C1.*D2C2.*D1));


 140  B2=(B12.*C1.*(C1.*D3C3.*D1)+B11.*C1.*(C2.*D1C1.*D2));


 141  Lambda=(A1.*dXa+A2.*dYa+B1.*dXb+B2.*dYb)./(A1.*A1+A2.*A2+B1.*B1+B2.*B2);


 142 


 143  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


 144  %Projection for compatible displacements


 145  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


 146  Ua=dXaLambda.*A1;


 147  Va=dYaLambda.*A2;


 148  Ub=dXbLambda.*B1;


 149  Vb=dYbLambda.*B2;


 150 


 151  %%%%%%%%%%%%%%%%%%%%%%%%%%%%


 152  %Calculations of displacements and error


 153  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


 154  U=(A22.*D2.*UaA12.*D2.*VaB22.*C2.*Ub+B12.*C2.*Vb)./(C1.*D2C2.*D1);


 155  V=(A21.*D3.*UaA11.*D3.*VaB21.*C3.*Ub+B11.*C3.*Vb)./(C3.*D1C1.*D3);


 156  W=(A22.*D1.*UaA12.*D1.*VaB22.*C1.*Ub+B12.*C1.*Vb)./(C2.*D1C1.*D2);


 157  W1=(A21.*D1.*Ua+A11.*D1.*Va+B21.*C1.*UbB11.*C1.*Vb)./(C1.*D3C3.*D1);


 158 


 159  error=sqrt((A1.*dXa+A2.*dYa+B1.*dXb+B2.*dYb).*(A1.*dXa+A2.*dYa+B1.*dXb+B2.*dYb)./(A1.*A1+A2.*A2+B1.*B1+B2.*B2));


 160 


 161  ind_error=(find(error<thresh_patch));


 162  U=U(ind_error);


 163  V=V(ind_error);


 164  W=W(ind_error);%correction for water interface


 165  error=error(ind_error);


 166 


 167  %create nc grid file


 168  Result.ListGlobalAttribute={'nb_coord','nb_dim','constant_pixcm','absolut_time_T0','hart','dt','civ'};


 169  Result.nb_coord=3;%grid file, no velocity


 170  Result.nb_dim=2;


 171  Result.constant_pixcm=0;%no linear correspondance with images


 172  Result.absolut_time_T0=time;%absolute time of the field


 173  Result.hart=0;


 174  Result.dt=dt;%time interval for image correlation (put by default)


 175  % cte.title='grid';


 176  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.VarDimIndex: {[1] [1] [1] [1] [1] [1]}


 181  Result.vec_X= grid_phys1(ind_error,1);


 182  Result.vec_Y= grid_phys1(ind_error,2);


 183  Result.vec_U=U/dt;


 184  Result.vec_V=V/dt;


 185  Result.vec_W=W/dt;


 186  Result.vec_E=error;


 187  % error=write_netcdf(file_st,cte,fieldlabels,grid_phys);


 188  error=struct2nc(file_st,Result);


 189  display([file_st ' written'])


 190 


 191 


 192 


[38]  193  %'pxcm_tsai': find differentials of the Tsai calibration


 194  %


 195  function [A11,A12,A13,A21,A22,A23]=pxcm_tsai(a,var_phys)


 196  a_read=a;


[8]  197 


[38]  198  R=(a.R)';


[8]  199 


[38]  200  x=var_phys(:,1);


 201  y=var_phys(:,2);


[8]  202 


[38]  203  if isfield(a,'PlanePos')


 204  prompt={'Plane 1 Index','Plane 2 Index'};


 205  Rep=inputdlg(prompt,'Target displacement test');


 206  Z1=str2double(Rep(1));


 207  Z2=str2double(Rep(2));


 208  z=(a.PlanePos(Z2,3)+a.PlanePos(Z1,3))/2


 209  else


 210  z=0;


 211  end


[8]  212 


[38]  213  %transform coeff for differentiels


 214  a.C11=R(1)*R(8)R(2)*R(7);


 215  a.C12=R(2)*R(7)R(1)*R(8);


 216  a.C21=R(4)*R(8)R(5)*R(7);


 217  a.C22=R(5)*R(7)R(4)*R(8);


 218  a.C1x=R(3)*R(7)R(9)*R(1);


 219  a.C1y=R(3)*R(8)R(9)*R(2);


 220  a.C2x=R(6)*R(7)R(9)*R(4);


 221  a.C2y=R(6)*R(8)R(9)*R(5);


[8]  222 


 223 


[38]  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*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;


[8]  232 


[38]  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  %


 263  % %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;


 269  % A13=(a.f*(a.C1x*x+a.C1y*y+R(3)*a.TzR(9)*a.Tx)./denom)/a.dpx;


 270  % A23=(a.f*(a.C2x*x+a.C2y*y+R(6)*a.TzR(9)*a.Ty)./denom)/a.dpy;


 271  %


[8]  272 


 273 


 274 


[38]  275 


 276 


 277 


 278 


 279 


 280 


 281 

