source: trunk/src/RUN_STLIN.m @ 38

Last change on this file since 38 was 38, checked in by sommeria, 11 years ago

field transforms put in subdir transform_field. cleaning of obsolete functions

File size: 9.7 KB
Line 
1%'RUN_STLIN': combine velocity fields for stereo PIV
2% file_A,file_B: input velocity files
3%vel_type: string ='civ1' or 'civ2'
4function 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     msgbox_uvmat('ERROR','The size of image A needs to be defined in the xml file ImaDoc')
29     return
30 elseif isempty(npxB) || isempty(npyB)
31      msgbox_uvmat('ERROR','The size of image B needs to be defined in the xml file ImaDoc')
32     return
33 end
34 if isfield(XmlDataA,'GeometryCalib')
35     tsaiA=XmlDataA.GeometryCalib;
36 else
37     msgbox_uvmat('ERROR','no geometric calibration available for image A')
38     return
39 end
40 if isfield(XmlDataB,'GeometryCalib')
41     tsaiB=XmlDataB.GeometryCalib;
42 else
43     msgbox_uvmat('ERROR','no geometric calibration available for image B')
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
86Field.X=Field.X(find(Field.FF==0));
87Field.Y=Field.Y(find(Field.FF==0));
88Field.U=Field.U(find(Field.FF==0));
89Field.V=Field.V(find(Field.FF==0));
90dXa= griddata_uvmat(Field.X,Field.Y,Field.U,XimaA,YimaA);
91dYa= griddata_uvmat(Field.X,Field.Y,Field.V,XimaA,YimaA);
92dt=Field.dt;
93time=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);
99if ~isequal(Field.dt,dt)
100    msgbox_uvmat('ERROR','different time intervals for the two velocity fields ')
101     return
102end
103if ~isequal(Field.Time,time)
104    msgbox_uvmat('ERROR','different times for the two velocity fields ')
105     return
106end
107%interpolate on XimaB
108Field.X=Field.X(find(Field.FF==0));
109Field.Y=Field.Y(find(Field.FF==0));
110Field.U=Field.U(find(Field.FF==0));
111Field.V=Field.V(find(Field.FF==0));
112dXb=griddata_uvmat(Field.X,Field.Y,Field.U,XimaB,YimaB);
113dYb=griddata_uvmat(Field.X,Field.Y,Field.V,XimaB,YimaB);
114%eliminate Not-a-Number
115ind_Nan=find(and(~isnan(dXa),~isnan(dXb)));
116dXa=dXa(ind_Nan);
117dYa=dYa(ind_Nan);
118dXb=dXb(ind_Nan);
119dYb=dYb(ind_Nan);
120grid_phys1(:,1)=grid_real_x(ind_Nan);
121grid_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
131C1=A11.*A22-A12.*A21;
132C2=A13.*A22-A12.*A23;
133C3=A13.*A21-A11.*A23;
134D1=B11.*B22-B12.*B21;
135D2=B13.*B22-B12.*B23;
136D3=B13.*B21-B11.*B23;
137A1=(A22.*D1.*(C1.*D3-C3.*D1)+A21.*D1.*(C2.*D1-C1.*D2));
138A2=(A12.*D1.*(C3.*D1-C1.*D3)+A11.*D1.*(C1.*D2-C2.*D1));
139B1=(B22.*C1.*(C3.*D1-C1.*D3)+B21.*C1.*(C1.*D2-C2.*D1));
140B2=(B12.*C1.*(C1.*D3-C3.*D1)+B11.*C1.*(C2.*D1-C1.*D2));
141Lambda=(A1.*dXa+A2.*dYa+B1.*dXb+B2.*dYb)./(A1.*A1+A2.*A2+B1.*B1+B2.*B2);
142
143%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
144%Projection for compatible displacements
145%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
146Ua=dXa-Lambda.*A1;
147Va=dYa-Lambda.*A2;
148Ub=dXb-Lambda.*B1;
149Vb=dYb-Lambda.*B2;
150
151%%%%%%%%%%%%%%%%%%%%%%%%%%%%
152%Calculations of displacements and error
153%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
154U=(A22.*D2.*Ua-A12.*D2.*Va-B22.*C2.*Ub+B12.*C2.*Vb)./(C1.*D2-C2.*D1);
155V=(A21.*D3.*Ua-A11.*D3.*Va-B21.*C3.*Ub+B11.*C3.*Vb)./(C3.*D1-C1.*D3);
156W=(A22.*D1.*Ua-A12.*D1.*Va-B22.*C1.*Ub+B12.*C1.*Vb)./(C2.*D1-C1.*D2);
157W1=(-A21.*D1.*Ua+A11.*D1.*Va+B21.*C1.*Ub-B11.*C1.*Vb)./(C1.*D3-C3.*D1);
158
159error=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
161ind_error=(find(error<thresh_patch));
162U=U(ind_error);
163V=V(ind_error);
164W=W(ind_error);%correction for water interface
165error=error(ind_error);
166
167%create nc grid file
168Result.ListGlobalAttribute={'nb_coord','nb_dim','constant_pixcm','absolut_time_T0','hart','dt','civ'};
169Result.nb_coord=3;%grid file, no velocity
170Result.nb_dim=2;
171Result.constant_pixcm=0;%no linear correspondance with images
172Result.absolut_time_T0=time;%absolute time of the field
173Result.hart=0;
174Result.dt=dt;%time interval for image correlation (put  by default)
175% cte.title='grid';
176Result.civ=0;%not a civ file (no direct correspondance with an image)
177Result.ListDimName={'nb_vectors'}
178Result.DimValue=length(U);
179Result.ListVarName={'vec_X';'vec_Y';'vec_U';'vec_V';'vec_W';'vec_E'};
180Result.VarDimIndex: {[1]  [1]  [1]  [1]  [1]  [1]}
181Result.vec_X= grid_phys1(ind_error,1);
182Result.vec_Y= grid_phys1(ind_error,2);
183Result.vec_U=U/dt;
184Result.vec_V=V/dt;
185Result.vec_W=W/dt;
186Result.vec_E=error;
187% error=write_netcdf(file_st,cte,fieldlabels,grid_phys);
188error=struct2nc(file_st,Result);
189display([file_st ' written'])
190
191
192
193%'pxcm_tsai': find differentials of the Tsai calibration
194%
195function [A11,A12,A13,A21,A22,A23]=pxcm_tsai(a,var_phys)
196a_read=a;
197
198R=(a.R)';
199
200x=var_phys(:,1);
201y=var_phys(:,2);
202
203if isfield(a,'PlanePos')
204    prompt={'Plane 1 Index','Plane 2 Index'};
205    Rep=inputdlg(prompt,'Target displacement test');
206    Z1=str2double(Rep(1));
207    Z2=str2double(Rep(2));
208    z=(a.PlanePos(Z2,3)+a.PlanePos(Z1,3))/2
209else
210    z=0;
211end
212
213%transform coeff for differentiels
214a.C11=R(1)*R(8)-R(2)*R(7);
215a.C12=R(2)*R(7)-R(1)*R(8);
216a.C21=R(4)*R(8)-R(5)*R(7);
217a.C22=R(5)*R(7)-R(4)*R(8);
218a.C1x=R(3)*R(7)-R(9)*R(1);
219a.C1y=R(3)*R(8)-R(9)*R(2);
220a.C2x=R(6)*R(7)-R(9)*R(4);
221a.C2y=R(6)*R(8)-R(9)*R(5);
222
223
224%dependence in x,y
225denom=(R(7)*x+R(8)*y+R(9)*z+a.Tz).*(R(7)*x+R(8)*y+R(9)*z+a.Tz);
226A11=(a.f*a.sx*(a.C11*y-a.C1x*z+R(1)*a.Tz-R(7)*a.Tx)./denom)/a.dpx;
227A12=(a.f*a.sx*(a.C12*x-a.C1y*z+R(2)*a.Tz-R(8)*a.Tx)./denom)/a.dpx;
228A21=(a.f*a.sx*(a.C21*y-a.C2x*z+R(4)*a.Tz-R(7)*a.Ty)./denom)/a.dpy;
229A22=(a.f*(a.C22*x-a.C2y*z+R(5)*a.Tz-R(8)*a.Ty)./denom)/a.dpy;
230A13=(a.f*(a.C1x*x+a.C1y*y+R(3)*a.Tz-R(9)*a.Tx)./denom)/a.dpx;
231A23=(a.f*(a.C2x*x+a.C2y*y+R(6)*a.Tz-R(9)*a.Ty)./denom)/a.dpy;
232
233%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
234%Old Version for z=0
235%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
236% %'camera' coordinates
237% xc=R(1)*x+R(2)*y+a.Tx;
238% yc=R(4)*x+R(5)*y+a.Ty;
239% zc=R(7)*x+R(8)*y+a.Tz;
240% %undistorted image coordinates
241% Xu=a.f*xc./zc;
242% Yu=a.f*yc./zc;
243% %distorted image coordinates
244% distortion=(a.kappa1)*(Xu.*Xu+Yu.*Yu)+1; %!! intégrer derivation kappa
245% % distortion=1;
246% Xd=Xu./distortion;
247% Yd=Yu./distortion;
248% %pixel coordinates
249% X=Xd*a.sx/a.dpx+a.Cx;
250% Y=Yd/a.dpy+a.Cy;
251%
252% %transform coeff for differentiels
253% a.C11=R(1)*R(8)-R(2)*R(7);
254% a.C12=R(2)*R(7)-R(1)*R(8);
255% a.C21=R(4)*R(8)-R(5)*R(7);
256% a.C22=R(5)*R(7)-R(4)*R(8);
257% a.C1x=R(3)*R(7)-R(9)*R(1);
258% a.C1y=R(3)*R(8)-R(9)*R(2);
259% a.C2x=R(6)*R(7)-R(9)*R(4);
260% a.C2y=R(6)*R(8)-R(9)*R(5);
261%
262%
263% %dependence in x,y
264% denom=(R(7)*x+R(8)*y+a.Tz).*(R(7)*x+R(8)*y+a.Tz);
265% A11=(a.f*a.sx*(a.C11*y+R(1)*a.Tz-R(7)*a.Tx)./denom)/a.dpx;
266% A12=(a.f*a.sx*(a.C12*x+R(2)*a.Tz-R(8)*a.Tx)./denom)/a.dpx;
267% A21=(a.f*a.sx*(a.C21*y+R(4)*a.Tz-R(7)*a.Ty)./denom)/a.dpy;
268% A22=(a.f*(a.C22*x+R(5)*a.Tz-R(8)*a.Ty)./denom)/a.dpy;
269% A13=(a.f*(a.C1x*x+a.C1y*y+R(3)*a.Tz-R(9)*a.Tx)./denom)/a.dpx;
270% A23=(a.f*(a.C2x*x+a.C2y*y+R(6)*a.Tz-R(9)*a.Ty)./denom)/a.dpy;
271%
272
273
274
275
276
277
278
279
280
281
Note: See TracBrowser for help on using the repository browser.