source: trunk/src/RUN_STLIN.m @ 849

Last change on this file since 849 was 809, checked in by g7moreau, 7 years ago
  • Add license
File size: 10.2 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'
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
23function 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
111if 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));
116end
117%interpolate on the grid common to both images in phys coordinates
118dXa= griddata_uvmat(Field.X,Field.Y,Field.U,XimaA,YimaA);
119dYa= griddata_uvmat(Field.X,Field.Y,Field.V,XimaA,YimaA);
120dt=Field.dt;
121time=Field.Time;
122
123%read field B
124[Field,VelTypeOut]=read_civxdata(file_B,[],vel_type);
125if ~isequal(Field.dt,dt)
126    msgbox_uvmat('ERROR','different time intervals for the two velocity fields ')
127     return
128end
129if ~isequal(Field.Time,time)
130    msgbox_uvmat('ERROR','different times for the two velocity fields ')
131     return
132end
133%removes false vectors
134if isfield(Field,'FF')
135Field.X=Field.X(find(Field.FF==0));
136Field.Y=Field.Y(find(Field.FF==0));
137Field.U=Field.U(find(Field.FF==0));
138Field.V=Field.V(find(Field.FF==0));
139end
140%interpolate on XimaB
141dXb=griddata_uvmat(Field.X,Field.Y,Field.U,XimaB,YimaB);
142dYb=griddata_uvmat(Field.X,Field.Y,Field.V,XimaB,YimaB);
143%eliminate Not-a-Number
144ind_Nan=find(and(~isnan(dXa),~isnan(dXb)));
145dXa=dXa(ind_Nan);
146dYa=dYa(ind_Nan);
147dXb=dXb(ind_Nan);
148dYb=dYb(ind_Nan);
149grid_phys1(:,1)=grid_real_x(ind_Nan);
150grid_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
158C1=A11.*A22-A12.*A21;
159C2=A13.*A22-A12.*A23;
160C3=A13.*A21-A11.*A23;
161D1=B11.*B22-B12.*B21;
162D2=B13.*B22-B12.*B23;
163D3=B13.*B21-B11.*B23;
164A1=(A22.*D1.*(C1.*D3-C3.*D1)+A21.*D1.*(C2.*D1-C1.*D2));
165A2=(A12.*D1.*(C3.*D1-C1.*D3)+A11.*D1.*(C1.*D2-C2.*D1));
166B1=(B22.*C1.*(C3.*D1-C1.*D3)+B21.*C1.*(C1.*D2-C2.*D1));
167B2=(B12.*C1.*(C1.*D3-C3.*D1)+B11.*C1.*(C2.*D1-C1.*D2));
168Lambda=(A1.*dXa+A2.*dYa+B1.*dXb+B2.*dYb)./(A1.*A1+A2.*A2+B1.*B1+B2.*B2);
169
170%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
171%Projection for compatible displacements
172%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
173Ua=dXa-Lambda.*A1;
174Va=dYa-Lambda.*A2;
175Ub=dXb-Lambda.*B1;
176Vb=dYb-Lambda.*B2;
177
178%%%%%%%%%%%%%%%%%%%%%%%%%%%%
179%Calculations of displacements and error
180%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
181U=(A22.*D2.*Ua-A12.*D2.*Va-B22.*C2.*Ub+B12.*C2.*Vb)./(C1.*D2-C2.*D1);
182V=(A21.*D3.*Ua-A11.*D3.*Va-B21.*C3.*Ub+B11.*C3.*Vb)./(C3.*D1-C1.*D3);
183W=(A22.*D1.*Ua-A12.*D1.*Va-B22.*C1.*Ub+B12.*C1.*Vb)./(C2.*D1-C1.*D2);
184W1=(-A21.*D1.*Ua+A11.*D1.*Va+B21.*C1.*Ub-B11.*C1.*Vb)./(C1.*D3-C3.*D1);
185
186error=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
188ind_error=(find(error<thresh_patch));
189U=U(ind_error);
190V=V(ind_error);
191W=W(ind_error);%correction for water interface
192error=error(ind_error);
193
194%create nc grid file
195Result.ListGlobalAttribute={'nb_coord','nb_dim','constant_pixcm','absolut_time_T0','hart','dt','civ'};
196Result.nb_coord=3;%grid file, no velocity
197Result.nb_dim=2;
198Result.constant_pixcm=0;%no linear correspondance with images
199Result.absolut_time_T0=time;%absolute time of the field
200Result.hart=0;
201Result.dt=dt;%time interval for image correlation (put  by default)
202% cte.title='grid';
203Result.civ=0;%not a civ file (no direct correspondance with an image)
204% Result.ListDimName={'nb_vectors'}
205% Result.DimValue=length(U);
206Result.ListVarName={'vec_X','vec_Y','vec_U','vec_V','vec_W','vec_E'};
207Result.VarDimName={'nb_vectors','nb_vectors','nb_vectors','nb_vectors','nb_vectors','nb_vectors'}
208Result.vec_X= grid_phys1(ind_error,1);
209Result.vec_Y= grid_phys1(ind_error,2);
210Result.vec_U=U/dt;
211Result.vec_V=V/dt;
212Result.vec_W=W/dt;
213Result.vec_E=error;
214% error=write_netcdf(file_st,cte,fieldlabels,grid_phys);
215error=struct2nc(file_st,Result);
216display([file_st ' written'])
217
218
219
220%'pxcm_tsai': find differentials of the Tsai calibration
221function [A11,A12,A13,A21,A22,A23]=pxcm_tsai(a,var_phys)
222R=(a.R)';
223
224x=var_phys(:,1);
225y=var_phys(:,2);
226
227if 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
233else
234    z=0;
235end
236
237%transform coeff for differentiels
238a.C11=R(1)*R(8)-R(2)*R(7);
239a.C12=R(2)*R(7)-R(1)*R(8);
240a.C21=R(4)*R(8)-R(5)*R(7);
241a.C22=R(5)*R(7)-R(4)*R(8);
242a.C1x=R(3)*R(7)-R(9)*R(1);
243a.C1y=R(3)*R(8)-R(9)*R(2);
244a.C2x=R(6)*R(7)-R(9)*R(4);
245a.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
257denom=(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));
258A11=(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);
259A12=(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);
260A21=(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);
261A22=(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);
262A13=(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);
263A23=(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
Note: See TracBrowser for help on using the repository browser.