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)


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


29  return


30  elseif isempty(npxB)  isempty(npyB)


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


32  return


33  end


34  if isfield(XmlDataA,'GeometryCalib')


35  tsaiA=XmlDataA.GeometryCalib;


36  else


37  warndlg_uvmat('no geometric calibration available for image A','ERROR')


38  return


39  end


40  if isfield(XmlDataB,'GeometryCalib')


41  tsaiB=XmlDataB.GeometryCalib;


42  else


43  warndlg_uvmat('no geometric calibration available for image B','ERROR')


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)


100  warndlg_uvmat('different time intervals for the two velocity fields ','ERROR')


101  return


102  end


103  if ~isequal(Field.Time,time)


104  warndlg_uvmat('different times for the two velocity fields ','ERROR')


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 


193 


194 


195 


196 


197 


198 


199 


200 


201 


202 

