Home > . > RUN_STLIN.m

# RUN_STLIN

## PURPOSE

'RUN_STLIN': combine velocity fields for stereo PIV

## SYNOPSIS

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

## DESCRIPTION

    'RUN_STLIN': combine velocity fields for stereo PIV
     file_A,file_B: input velocity files
    vel_type: string ='civ1' or 'civ2'

## CROSS-REFERENCE INFORMATION

This function calls:

- griddata_uvmat 'griddata_uvmat': function griddata_uvmat(vec2_X,vec2_Y,vec2_U,vec_X,vec_Y,'linear')
- imadoc2struct 'imadoc2struct': reads the xml file for image documentation
- phys_XYZ 'phys_XYZ':transforms image (px) to real world (phys) coordinates using geometric calibration parameters
- px_XYZ 'px_XYZ': transform phys coordinates to image coordinates (px)
- pxcm_tsai 'pxcm_tsai': find differentials of the Tsai calibration
- read_civxdata 'read_civxdata': reads civx data from netcdf files
- struct2nc 'struct2nc': create a netcdf file from a Matlab structure
- warndlg_uvmat 'warndlg_uvmat': display warning message (error, warning, confirmation) in a given figure

This function is called by:

- civ 'civ': function associated with the interface 'civ.fig' for PIV, spline interpolation and stereo PIV (patch)

## SOURCE CODE

    0001 %'RUN_STLIN': combine velocity fields for stereo PIV
    0002 % file_A,file_B: input velocity files
    0003 %vel_type: string ='civ1' or 'civ2'
    0004 function RUN_STLIN(file_A,file_B,vel_type,file_st,nx_patch,ny_patch,thresh_patch,fileAxml,fileBxml)
    0005
    0006  [XmlDataA,error]=imadoc2struct(fileAxml);
    0007  [XmlDataB,error]=imadoc2struct(fileBxml);
    0008  npxA=[]; npyA=[]; pxB=[]; npyB=[];
    0009  if isfield(XmlDataA,'Camera') && isfield(XmlDataB,'Camera')
    0010       if isfield(XmlDataA.Camera,'ImageSize')&& isfield(XmlDataB.Camera,'ImageSize')
    0011           ImageSizeA=XmlDataA.Camera.ImageSize;
    0012           ImageSizeB=XmlDataB.Camera.ImageSize;
    0013           if ~isempty(ImageSizeA)&& ~isempty(ImageSizeB)
    0014                xindex=findstr(ImageSizeA,'x');
    0015                if length(xindex)>=2
    0016                     npxA=str2num(ImageSizeA(1:xindex(1)-1));
    0017                     npyA=str2num(ImageSizeA(xindex(1)+1:xindex(2)-1));
    0018                end
    0019                xindex=findstr(ImageSizeB,'x');
    0020                if length(xindex)>=2
    0021                     npxB=str2num(ImageSizeB(1:xindex(1)-1));
    0022                     npyB=str2num(ImageSizeB(xindex(1)+1:xindex(2)-1));
    0023                end
    0024           end
    0025      end
    0026  end
    0027  if isempty(npxA) ||isempty(npxB)
    0028      warndlg_uvmat('The size of image A needs to be defined in the xml file ImaDoc','ERROR')
    0029      return
    0030  elseif isempty(npxB) || isempty(npyB)
    0031       warndlg_uvmat('The size of image B needs to be defined in the xml file ImaDoc','ERROR')
    0032      return
    0033  end
    0034  if isfield(XmlDataA,'GeometryCalib')
    0035      tsaiA=XmlDataA.GeometryCalib;
    0036  else
    0037      warndlg_uvmat('no geometric calibration available for image A','ERROR')
    0038      return
    0039  end
    0040  if isfield(XmlDataB,'GeometryCalib')
    0041      tsaiB=XmlDataB.GeometryCalib;
    0042  else
    0043      warndlg_uvmat('no geometric calibration available for image B','ERROR')
    0044      return
    0045  end
    0046
    0047  %corners of each image in real coordinates:
    0048  cornerA(:,1)=[0 0 npxA npxA]';%x positions
    0049  cornerA(:,2)=[0 npyA 0 npyA]';%y positions
    0050  cornerB(:,1)=[0 0 npxB npxB]';%x positions
    0051  cornerB(:,2)=[0 npyB 0 npyB]';%y positions
    0052 [xyA(:,1),xyA(:,2)]=phys_XYZ(tsaiA,cornerA(:,1),cornerA(:,2));
    0053 [xyB(:,1),xyB(:,2)]=phys_XYZ(tsaiB,cornerB(:,1),cornerB(:,2));
    0054  max_x=max(max(xyA(:,1)),max(xyB(:,1)));%maximum on the 4 corners of the the images
    0055  min_x=min(min(xyA(:,1)),min(xyB(:,1)));%minimum on the 4 corners of the the images
    0056  max_y=max(max(xyA(:,2)),max(xyB(:,2)));
    0057  min_y=min(min(xyA(:,2)),min(xyB(:,2)));
    0058  array_realx=[min_x:(max_x-min_x)/(nx_patch-1):max_x];
    0059  array_realy=[min_y:(max_y-min_y)/(ny_patch-1):max_y];
    0060  [grid_realx,grid_realy]=meshgrid(array_realx,array_realy);
    0061  grid_real(:,1)=reshape(grid_realx,nx_patch*ny_patch,1);
    0062  grid_real(:,2)=reshape(grid_realy,nx_patch*ny_patch,1);
    0063  grid_real(:,3)=zeros(nx_patch*ny_patch,1);
    0064 [grid_imaA(:,1),grid_imaA(:,2)]=px_XYZ(tsaiA,grid_real(:,1),grid_real(:,2));
    0065 [grid_imaB(:,1),grid_imaB(:,2)]=px_XYZ(tsaiB,grid_real(:,1),grid_real(:,2));
    0066
    0067  flagA=grid_imaA(:,1)>0 & grid_imaA(:,1)<npxA & grid_imaA(:,2)>0 & grid_imaA(:,2)<npyA;
    0068  flagB=grid_imaB(:,1)>0 & grid_imaB(:,1)<npxB & grid_imaB(:,2)>0 & grid_imaB(:,2)<npyB;
    0069  ind_good=find(flagA==1&flagB==1);
    0070  XimaA=grid_imaA(ind_good,1);
    0071  YimaA=grid_imaA(ind_good,2);
    0072  XimaB=grid_imaB(ind_good,1);
    0073  YimaB=grid_imaB(ind_good,2);
    0074  grid_real_x=grid_real(ind_good,1);
    0075  grid_real_y=grid_real(ind_good,2);
    0076
    0077 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    0078 % %read the velocity fields
    0079 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    0080 %
    0081 % [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]...
    0082 %     =read_vel({filecell_ncA},{vel_type});
    0083 %read field A
    0084 [Field,VelTypeOut]=read_civxdata(file_A,[],vel_type);
    0085 %interpolate on XimaA
    0086 Field.X=Field.X(find(Field.FF==0));
    0087 Field.Y=Field.Y(find(Field.FF==0));
    0088 Field.U=Field.U(find(Field.FF==0));
    0089 Field.V=Field.V(find(Field.FF==0));
    0090 dXa= griddata_uvmat(Field.X,Field.Y,Field.U,XimaA,YimaA);
    0091 dYa= griddata_uvmat(Field.X,Field.Y,Field.V,XimaA,YimaA);
    0092 dt=Field.dt;
    0093 time=Field.Time;
    0094
    0095 %read field B
    0096 % [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]...
    0097 %     =read_vel({file_B},{vel_type});
    0098 [Field,VelTypeOut]=read_civxdata(file_B,FieldNames,vel_type);
    0099 if ~isequal(Field.dt,dt)
    0100     warndlg_uvmat('different time intervals for the two velocity fields ','ERROR')
    0101      return
    0102 end
    0103 if ~isequal(Field.Time,time)
    0104     warndlg_uvmat('different times for the two velocity fields ','ERROR')
    0105      return
    0106 end
    0107 %interpolate on XimaB
    0108 Field.X=Field.X(find(Field.FF==0));
    0109 Field.Y=Field.Y(find(Field.FF==0));
    0110 Field.U=Field.U(find(Field.FF==0));
    0111 Field.V=Field.V(find(Field.FF==0));
    0112 dXb=griddata_uvmat(Field.X,Field.Y,Field.U,XimaB,YimaB);
    0113 dYb=griddata_uvmat(Field.X,Field.Y,Field.V,XimaB,YimaB);
    0114 %eliminate Not-a-Number
    0115 ind_Nan=find(and(~isnan(dXa),~isnan(dXb)));
    0116 dXa=dXa(ind_Nan);
    0117 dYa=dYa(ind_Nan);
    0118 dXb=dXb(ind_Nan);
    0119 dYb=dYb(ind_Nan);
    0120 grid_phys1(:,1)=grid_real_x(ind_Nan);
    0121 grid_phys1(:,2)=grid_real_y(ind_Nan);
    0122
    0123 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    0124 %compute the coefficients
    0125 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    0126
    0127
    0128 [A11,A12,A13,A21,A22,A23]=pxcm_tsai(tsaiA,grid_phys1);
    0129 [B11,B12,B13,B21,B22,B23]=pxcm_tsai(tsaiB,grid_phys1);
    0130
    0131 C1=A11.*A22-A12.*A21;
    0132 C2=A13.*A22-A12.*A23;
    0133 C3=A13.*A21-A11.*A23;
    0134 D1=B11.*B22-B12.*B21;
    0135 D2=B13.*B22-B12.*B23;
    0136 D3=B13.*B21-B11.*B23;
    0137 A1=(A22.*D1.*(C1.*D3-C3.*D1)+A21.*D1.*(C2.*D1-C1.*D2));
    0138 A2=(A12.*D1.*(C3.*D1-C1.*D3)+A11.*D1.*(C1.*D2-C2.*D1));
    0139 B1=(B22.*C1.*(C3.*D1-C1.*D3)+B21.*C1.*(C1.*D2-C2.*D1));
    0140 B2=(B12.*C1.*(C1.*D3-C3.*D1)+B11.*C1.*(C2.*D1-C1.*D2));
    0141 Lambda=(A1.*dXa+A2.*dYa+B1.*dXb+B2.*dYb)./(A1.*A1+A2.*A2+B1.*B1+B2.*B2);
    0142
    0143 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    0144 %Projection for compatible displacements
    0145 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    0146 Ua=dXa-Lambda.*A1;
    0147 Va=dYa-Lambda.*A2;
    0148 Ub=dXb-Lambda.*B1;
    0149 Vb=dYb-Lambda.*B2;
    0150
    0151 %%%%%%%%%%%%%%%%%%%%%%%%%%%%
    0152 %Calculations of displacements and error
    0153 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    0154 U=(A22.*D2.*Ua-A12.*D2.*Va-B22.*C2.*Ub+B12.*C2.*Vb)./(C1.*D2-C2.*D1);
    0155 V=(A21.*D3.*Ua-A11.*D3.*Va-B21.*C3.*Ub+B11.*C3.*Vb)./(C3.*D1-C1.*D3);
    0156 W=(A22.*D1.*Ua-A12.*D1.*Va-B22.*C1.*Ub+B12.*C1.*Vb)./(C2.*D1-C1.*D2);
    0157 W1=(-A21.*D1.*Ua+A11.*D1.*Va+B21.*C1.*Ub-B11.*C1.*Vb)./(C1.*D3-C3.*D1);
    0158
    0159 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));
    0160
    0161 ind_error=(find(error<thresh_patch));
    0162 U=U(ind_error);
    0163 V=V(ind_error);
    0164 W=W(ind_error);%correction for water interface
    0165 error=error(ind_error);
    0166
    0167 %create nc grid file
    0168 Result.ListGlobalAttribute={'nb_coord','nb_dim','constant_pixcm','absolut_time_T0','hart','dt','civ'};
    0169 Result.nb_coord=3;%grid file, no velocity
    0170 Result.nb_dim=2;
    0171 Result.constant_pixcm=0;%no linear correspondance with images
    0172 Result.absolut_time_T0=time;%absolute time of the field
    0173 Result.hart=0;
    0174 Result.dt=dt;%time interval for image correlation (put  by default)
    0175 % cte.title='grid';
    0176 Result.civ=0;%not a civ file (no direct correspondance with an image)
    0177 Result.ListDimName={'nb_vectors'}
    0178 Result.DimValue=length(U);
    0179 Result.ListVarName={'vec_X';'vec_Y';'vec_U';'vec_V';'vec_W';'vec_E'};
    0180 Result.VarDimIndex: {[1]  [1]  [1]  [1]  [1]  [1]}
    0181 Result.vec_X= grid_phys1(ind_error,1);
    0182 Result.vec_Y= grid_phys1(ind_error,2);
    0183 Result.vec_U=U/dt;
    0184 Result.vec_V=V/dt;
    0185 Result.vec_W=W/dt;
    0186 Result.vec_E=error;
    0187 % error=write_netcdf(file_st,cte,fieldlabels,grid_phys);
    0188 error=struct2nc(file_st,Result);
    0189 display([file_st ' written'])
    0190
    0191
    0192
    0193
    0194
    0195
    0196
    0197
    0198
    0199
    0200
    0201
    0202