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 2008-2014, LEGI UMR 5519 / CNRS UJF G-INP, Grenoble, France
|
---|
7 | % http://www.legi.grenoble-inp.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_x-min_x)/(nx_patch-1):max_x];
|
---|
87 | array_realy=[min_y:(max_y-min_y)/(ny_patch-1):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 Not-a-Number
|
---|
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.*A22-A12.*A21;
|
---|
159 | C2=A13.*A22-A12.*A23;
|
---|
160 | C3=A13.*A21-A11.*A23;
|
---|
161 | D1=B11.*B22-B12.*B21;
|
---|
162 | D2=B13.*B22-B12.*B23;
|
---|
163 | D3=B13.*B21-B11.*B23;
|
---|
164 | A1=(A22.*D1.*(C1.*D3-C3.*D1)+A21.*D1.*(C2.*D1-C1.*D2));
|
---|
165 | A2=(A12.*D1.*(C3.*D1-C1.*D3)+A11.*D1.*(C1.*D2-C2.*D1));
|
---|
166 | B1=(B22.*C1.*(C3.*D1-C1.*D3)+B21.*C1.*(C1.*D2-C2.*D1));
|
---|
167 | B2=(B12.*C1.*(C1.*D3-C3.*D1)+B11.*C1.*(C2.*D1-C1.*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=dXa-Lambda.*A1;
|
---|
174 | Va=dYa-Lambda.*A2;
|
---|
175 | Ub=dXb-Lambda.*B1;
|
---|
176 | Vb=dYb-Lambda.*B2;
|
---|
177 |
|
---|
178 | %%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
---|
179 | %Calculations of displacements and error
|
---|
180 | %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
---|
181 | U=(A22.*D2.*Ua-A12.*D2.*Va-B22.*C2.*Ub+B12.*C2.*Vb)./(C1.*D2-C2.*D1);
|
---|
182 | V=(A21.*D3.*Ua-A11.*D3.*Va-B21.*C3.*Ub+B11.*C3.*Vb)./(C3.*D1-C1.*D3);
|
---|
183 | W=(A22.*D1.*Ua-A12.*D1.*Va-B22.*C1.*Ub+B12.*C1.*Vb)./(C2.*D1-C1.*D2);
|
---|
184 | W1=(-A21.*D1.*Ua+A11.*D1.*Va+B21.*C1.*Ub-B11.*C1.*Vb)./(C1.*D3-C3.*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*y-a.C1x*z+R(1)*a.Tz-R(7)*a.Tx)./denom)/a.dpx;
|
---|
250 | % A12=(a.f*a.sx*(a.C12*x-a.C1y*z+R(2)*a.Tz-R(8)*a.Tx)./denom)/a.dpx;
|
---|
251 | % A21=(a.f*a.sx*(a.C21*y-a.C2x*z+R(4)*a.Tz-R(7)*a.Ty)./denom)/a.dpy;
|
---|
252 | % A22=(a.f*(a.C22*x-a.C2y*z+R(5)*a.Tz-R(8)*a.Ty)./denom)/a.dpy;
|
---|
253 | % A13=(a.f*(a.C1x*x+a.C1y*y+R(3)*a.Tz-R(9)*a.Tx)./denom)/a.dpx;
|
---|
254 | % A23=(a.f*(a.C2x*x+a.C2y*y+R(6)*a.Tz-R(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*y-a.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*x-a.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*y-a.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*x-a.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 |
|
---|