[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)
|
---|
| 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_x-min_x)/(nx_patch-1):max_x];
|
---|
| 59 | array_realy=[min_y:(max_y-min_y)/(ny_patch-1):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 Not-a-Number
|
---|
| 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.*A22-A12.*A21;
|
---|
| 132 | C2=A13.*A22-A12.*A23;
|
---|
| 133 | C3=A13.*A21-A11.*A23;
|
---|
| 134 | D1=B11.*B22-B12.*B21;
|
---|
| 135 | D2=B13.*B22-B12.*B23;
|
---|
| 136 | D3=B13.*B21-B11.*B23;
|
---|
| 137 | A1=(A22.*D1.*(C1.*D3-C3.*D1)+A21.*D1.*(C2.*D1-C1.*D2));
|
---|
| 138 | A2=(A12.*D1.*(C3.*D1-C1.*D3)+A11.*D1.*(C1.*D2-C2.*D1));
|
---|
| 139 | B1=(B22.*C1.*(C3.*D1-C1.*D3)+B21.*C1.*(C1.*D2-C2.*D1));
|
---|
| 140 | B2=(B12.*C1.*(C1.*D3-C3.*D1)+B11.*C1.*(C2.*D1-C1.*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=dXa-Lambda.*A1;
|
---|
| 147 | Va=dYa-Lambda.*A2;
|
---|
| 148 | Ub=dXb-Lambda.*B1;
|
---|
| 149 | Vb=dYb-Lambda.*B2;
|
---|
| 150 |
|
---|
| 151 | %%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
---|
| 152 | %Calculations of displacements and error
|
---|
| 153 | %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
---|
| 154 | U=(A22.*D2.*Ua-A12.*D2.*Va-B22.*C2.*Ub+B12.*C2.*Vb)./(C1.*D2-C2.*D1);
|
---|
| 155 | V=(A21.*D3.*Ua-A11.*D3.*Va-B21.*C3.*Ub+B11.*C3.*Vb)./(C3.*D1-C1.*D3);
|
---|
| 156 | W=(A22.*D1.*Ua-A12.*D1.*Va-B22.*C1.*Ub+B12.*C1.*Vb)./(C2.*D1-C1.*D2);
|
---|
| 157 | W1=(-A21.*D1.*Ua+A11.*D1.*Va+B21.*C1.*Ub-B11.*C1.*Vb)./(C1.*D3-C3.*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 |
|
---|