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 


5  %=======================================================================


6  % Copyright 20082014, LEGI UMR 5519 / CNRS UJF GINP, Grenoble, France


7  % http://www.legi.grenobleinp.fr


8  % Joel.Sommeria  Joel.Sommeria (A) legi.cnrs.fr


9  %


10  % This file is part of the toolbox UVMAT.


11  %


12  % UVMAT is free software; you can redistribute it and/or modify


13  % it under the terms of the GNU General Public License as published


14  % by the Free Software Foundation; either version 2 of the license,


15  % or (at your option) any later version.


16  %


17  % UVMAT is distributed in the hope that it will be useful,


18  % but WITHOUT ANY WARRANTY; without even the implied warranty of


19  % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the


20  % GNU General Public License (see LICENSE.txt) for more details.


21  %=======================================================================


22 


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


24 


25  [XmlDataA,error]=imadoc2struct(fileAxml);%read xml file associated to image series A


26  [XmlDataB,error]=imadoc2struct(fileBxml);%read xml file associated to image series B


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


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


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


30  ImageSizeA=XmlDataA.Camera.ImageSize;


31  ImageSizeB=XmlDataB.Camera.ImageSize;


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


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


34  if length(xindex)>=2


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


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


37  end


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


39  if length(xindex)>=2


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


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


42  end


43  end


44  end


45  end


46  %%%%%%%% added for Duran


47  if isfield(XmlDataA,'Npx')&&isfield(XmlDataA,'Npy')&&isfield(XmlDataB,'Npx')&&isfield(XmlDataB,'Npy')


48  npxA=XmlDataA.Npx;


49  npyA=XmlDataA.Npy;


50  npxB=XmlDataB.Npx;


51  npyB=XmlDataB.Npy;


52  end


53  %%%%%%%% added for Duran


54  if isempty(npxA) isempty(npxB)


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


56  return


57  elseif isempty(npxB)  isempty(npyB)


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


59  return


60  end


61  if isfield(XmlDataA,'GeometryCalib')


62  tsaiA=XmlDataA.GeometryCalib;


63  else


64  msgbox_uvmat('ERROR','no geometric calibration available for image A')


65  return


66  end


67  if isfield(XmlDataB,'GeometryCalib')


68  tsaiB=XmlDataB.GeometryCalib;


69  else


70  msgbox_uvmat('ERROR','no geometric calibration available for image B')


71  return


72  end


73 


74  %corners of each image in px coordinates:


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


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


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


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


79  %corners of each image in phys coordinates:


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


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


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


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


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


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


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


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


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


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


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


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


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


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


94 


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


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


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


98  XimaA=grid_imaA(ind_good,1);


99  YimaA=grid_imaA(ind_good,2);


100  XimaB=grid_imaB(ind_good,1);


101  YimaB=grid_imaB(ind_good,2);


102  grid_real_x=grid_real(ind_good,1);


103  grid_real_y=grid_real(ind_good,2);


104 


105  % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


106  % %read the velocity fields


107  % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


108  %read field A


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


110  %removes false vectors


111  if isfield(Field,'FF')


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


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


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


115  Field.V=Field.V(find(Field.FF==0));


116  end


117  %interpolate on the grid common to both images in phys coordinates


118  dXa= griddata_uvmat(Field.X,Field.Y,Field.U,XimaA,YimaA);


119  dYa= griddata_uvmat(Field.X,Field.Y,Field.V,XimaA,YimaA);


120  dt=Field.dt;


121  time=Field.Time;


122 


123  %read field B


124  [Field,VelTypeOut]=read_civxdata(file_B,[],vel_type);


125  if ~isequal(Field.dt,dt)


126  msgbox_uvmat('ERROR','different time intervals for the two velocity fields ')


127  return


128  end


129  if ~isequal(Field.Time,time)


130  msgbox_uvmat('ERROR','different times for the two velocity fields ')


131  return


132  end


133  %removes false vectors


134  if isfield(Field,'FF')


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


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


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


138  Field.V=Field.V(find(Field.FF==0));


139  end


140  %interpolate on XimaB


141  dXb=griddata_uvmat(Field.X,Field.Y,Field.U,XimaB,YimaB);


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


143  %eliminate NotaNumber


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


145  dXa=dXa(ind_Nan);


146  dYa=dYa(ind_Nan);


147  dXb=dXb(ind_Nan);


148  dYb=dYb(ind_Nan);


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


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


151 


152  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


153  %compute the differential coefficients of the geometric calibration


154  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


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


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


157 


158  C1=A11.*A22A12.*A21;


159  C2=A13.*A22A12.*A23;


160  C3=A13.*A21A11.*A23;


161  D1=B11.*B22B12.*B21;


162  D2=B13.*B22B12.*B23;


163  D3=B13.*B21B11.*B23;


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


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


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


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


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


169 


170  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


171  %Projection for compatible displacements


172  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


173  Ua=dXaLambda.*A1;


174  Va=dYaLambda.*A2;


175  Ub=dXbLambda.*B1;


176  Vb=dYbLambda.*B2;


177 


178  %%%%%%%%%%%%%%%%%%%%%%%%%%%%


179  %Calculations of displacements and error


180  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


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


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


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


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


185 


186  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));


187 


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


189  U=U(ind_error);


190  V=V(ind_error);


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


192  error=error(ind_error);


193 


194  %create nc grid file


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


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


197  Result.nb_dim=2;


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


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


200  Result.hart=0;


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


202  % cte.title='grid';


203  Result.civ=0;%not a civ file (no direct correspondance with an image)


204  % Result.ListDimName={'nb_vectors'}


205  % Result.DimValue=length(U);


206  Result.ListVarName={'vec_X','vec_Y','vec_U','vec_V','vec_W','vec_E'};


207  Result.VarDimName={'nb_vectors','nb_vectors','nb_vectors','nb_vectors','nb_vectors','nb_vectors'}


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


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


210  Result.vec_U=U/dt;


211  Result.vec_V=V/dt;


212  Result.vec_W=W/dt;


213  Result.vec_E=error;


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


215  error=struct2nc(file_st,Result);


216  display([file_st ' written'])


217 


218 


219 


220  %'pxcm_tsai': find differentials of the Tsai calibration


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


222  R=(a.R)';


223 


224  x=var_phys(:,1);


225  y=var_phys(:,2);


226 


227  if isfield(a,'PlanePos')


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


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


230  Z1=str2double(Rep(1));


231  Z2=str2double(Rep(2));


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


233  else


234  z=0;


235  end


236 


237  %transform coeff for differentiels


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


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


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


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


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


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


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


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


246 


247  % %dependence in x,y


248  % denom=(R(7)*x+R(8)*y+R(9)*z+a.Tz).*(R(7)*x+R(8)*y+R(9)*z+a.Tz);


249  % A11=(a.f*a.sx*(a.C11*ya.C1x*z+R(1)*a.TzR(7)*a.Tx)./denom)/a.dpx;


250  % A12=(a.f*a.sx*(a.C12*xa.C1y*z+R(2)*a.TzR(8)*a.Tx)./denom)/a.dpx;


251  % A21=(a.f*a.sx*(a.C21*ya.C2x*z+R(4)*a.TzR(7)*a.Ty)./denom)/a.dpy;


252  % A22=(a.f*(a.C22*xa.C2y*z+R(5)*a.TzR(8)*a.Ty)./denom)/a.dpy;


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


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


255 


256  %dependence in x,y


257  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));


258  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);


259  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);


260  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);


261  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);


262  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);


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


264 


265 


266 


267 


268 


269 


270 


271 

