source: trunk/src/series/civ2vel_3C.m @ 840

Last change on this file since 840 was 840, checked in by sommeria, 10 years ago
File size: 17.5 KB
Line 
1%'civ2vel_3C': combine velocity fields from two camerasto get three velocity components
2%------------------------------------------------------------------------
3% function ParamOut=civ2vel_3C(Param)
4%
5%OUTPUT
6% ParamOut: sets options in the GUI series.fig needed for the function
7%
8%INPUT:
9% In run mode, the input parameters are given as a Matlab structure Param copied from the GUI series.
10% In batch mode, Param is the name of the corresponding xml file containing the same information
11% when Param.Action.RUN=0 (as activated when the current Action is selected
12% in series), the function ouput paramOut set the activation of the needed GUI elements
13%
14% Param contains the elements:(use the menu bar command 'export/GUI config' in series to
15% see the current structure Param)
16%    .InputTable: cell of input file names, (several lines for multiple input)
17%                      each line decomposed as {RootPath,SubDir,Rootfile,NomType,Extension}
18%    .OutputSubDir: name of the subdirectory for data outputs
19%    .OutputDirExt: directory extension for data outputs
20%    .Action: .ActionName: name of the current activated function
21%             .ActionPath:   path of the current activated function
22%             .ActionExt: fct extension ('.m', Matlab fct, '.sh', compiled   Matlab fct
23%             .RUN =0 for GUI input, =1 for function activation
24%             .RunMode='local','background', 'cluster': type of function  use
25%             
26%    .IndexRange: set the file or frame indices on which the action must be performed
27%    .InputFields: sub structure describing the input fields withfields
28%              .FieldName: name(s) of the field
29%              .VelType: velocity type
30%              .FieldName_1: name of the second field in case of two input series
31%              .VelType_1: velocity type of the second field in case of two input series
32%              .Coord_y: name of y coordinate variable
33%              .Coord_x: name of x coordinate variable'
34
35%=======================================================================
36% Copyright 2008-2014, LEGI UMR 5519 / CNRS UJF G-INP, Grenoble, France
37%   http://www.legi.grenoble-inp.fr
38%   Joel.Sommeria - Joel.Sommeria (A) legi.cnrs.fr
39%
40%     This file is part of the toolbox UVMAT.
41%
42%     UVMAT is free software; you can redistribute it and/or modify
43%     it under the terms of the GNU General Public License as published
44%     by the Free Software Foundation; either version 2 of the license,
45%     or (at your option) any later version.
46%
47%     UVMAT is distributed in the hope that it will be useful,
48%     but WITHOUT ANY WARRANTY; without even the implied warranty of
49%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
50%     GNU General Public License (see LICENSE.txt) for more details.
51%=======================================================================
52
53function ParamOut=civ2vel_3C(Param)
54
55%% set the input elements needed on the GUI series when the function is selected in the menu ActionName or InputTable refreshed
56if isstruct(Param) && isequal(Param.Action.RUN,0)
57    ParamOut.AllowInputSort='on';% allow alphabetic sorting of the list of input file SubDir (options 'off'/'on', 'off' by default)
58    ParamOut.WholeIndexRange='off';% prescribes the file index ranges from min to max (options 'off'/'on', 'off' by default)
59    ParamOut.NbSlice='off'; %nbre of slices ('off' by default)
60    ParamOut.VelType='one';% menu for selecting the velocity type (options 'off'/'one'/'two',  'off' by default)
61    ParamOut.FieldName='off';% menu for selecting the field (s) in the input file(options 'off'/'one'/'two', 'off' by default)
62    ParamOut.FieldTransform = 'off';%use the phys  transform function without choice
63    %ParamOut.TransformPath=fullfile(fileparts(which('uvmat')),'transform_field');% path to transform functions (needed for compilation only)
64    ParamOut.ProjObject='off';%can use projection object(option 'off'/'on',
65    ParamOut.Mask='off';%can use mask option   (option 'off'/'on', 'off' by default)
66    ParamOut.OutputDirExt='.vel3C';%set the output dir extension
67    ParamOut.OutputFileMode='NbInput';% '=NbInput': 1 output file per input file index, '=NbInput_i': 1 file per input file index i, '=NbSlice': 1 file per slice
68      %check the input files
69    first_j=[];
70    if isfield(Param.IndexRange,'first_j'); first_j=Param.IndexRange.first_j; end
71    PairString='';
72    if isfield(Param.IndexRange,'PairString'); PairString=Param.IndexRange.PairString; end
73    [i1,i2,j1,j2] = get_file_index(Param.IndexRange.first_i,first_j,PairString);
74    FirstFileName=fullfile_uvmat(Param.InputTable{1,1},Param.InputTable{1,2},Param.InputTable{1,3},...
75        Param.InputTable{1,5},Param.InputTable{1,4},i1,i2,j1,j2);
76    if ~exist(FirstFileName,'file')
77        msgbox_uvmat('WARNING',['the first input file ' FirstFileName ' does not exist'])
78    elseif isequal(size(Param.InputTable,1),1) && ~isfield(Param,'ProjObject')
79        msgbox_uvmat('WARNING','You may need a projection object of type plane for merge_proj')
80    end
81    return
82end
83
84%%%%%%%%%%%% STANDARD PART (DO NOT EDIT) %%%%%%%%%%%%
85ParamOut=[]; %default output
86%% read input parameters from an xml file if input is a file name (batch mode)
87checkrun=1;
88if ischar(Param)
89    Param=xml2struct(Param);% read Param as input file (batch case)
90    checkrun=0;
91end
92hseries=findobj(allchild(0),'Tag','series');
93RUNHandle=findobj(hseries,'Tag','RUN');%handle of RUN button in GUI series
94WaitbarHandle=findobj(hseries,'Tag','Waitbar');%handle of waitbar in GUI series
95
96%% define the directory for result file (with path=RootPath{1})
97OutputDir=[Param.OutputSubDir Param.OutputDirExt];% subdirectory for output files
98%
99% if ~isfield(Param,'InputFields')
100%     Param.InputFields.FieldName='';
101% end
102
103%% calibration data and timing: read the ImaDoc files
104[XmlData,NbSlice_calib,time,errormsg]=read_multimadoc(RootPath,SubDir,RootFile,FileExt,i1_series,i2_series,j1_series,j2_series);
105if size(time,1)>1
106    diff_time=max(max(diff(time)));
107    if diff_time>0
108        disp_uvmat('WARNING',['times of series differ by (max) ' num2str(diff_time) ': the mean time is chosen in result'],checkrun)
109    end   
110end
111if ~isempty(errormsg)
112    disp_uvmat('WARNING',errormsg,checkrun)
113end
114time=mean(time,1); %averaged time taken for the merged field         
115if isfield(XmlData{1},'GeometryCalib')
116     tsaiA=XmlData{1}.GeometryCalib;
117 else
118     msgbox_uvmat('ERROR','no geometric calibration available for image A')
119     return
120 end
121 if isfield(XmlData{2},'GeometryCalib')
122     tsaiB=XmlData{2}.GeometryCalib;
123 else
124     msgbox_uvmat('ERROR','no geometric calibration available for image B')
125     return
126 end
127[filecell,i1_series,i2_series,j1_series,j2_series]=get_file_series(Param);
128
129%% grid of physical positions
130[x,y]=meshgrid(ObjectData.RangeX(1):ObjectData.DX:ObjectData.RangeX(2),ObjectData.RangeY(1):ObjectData.DY:ObjectData.RangeY(2));
131%camera coordinates:initialisation 2 cameras
132NbCamera=2;
133X=zeros(NbCamera,size(x,1),size(y,1));
134Y=zeros(NbCamera,size(x,1),size(y,1));
135for icamera=1:NbCamera
136%camera coordinates
137Calib=XmlData{1}.GeometryCalib;
138    xc=R(1)*Xphys+R(2)*Yphys+R(3)*Zphys+Calib.Tx_Ty_Tz(1);
139    yc=R(4)*Xphys+R(5)*Yphys+R(6)*Zphys+Calib.Tx_Ty_Tz(2);
140    zc=R(7)*Xphys+R(8)*Yphys+R(9)*Zphys+Calib.Tx_Ty_Tz(3);
141   
142    %undistorted image coordinates
143    X(icamera,:,:)=xc./zc;
144    Y(icamera,:,:)=yc./zc;
145end
146
147%% case of fixed planes (for moving interface, coeff need to be calculated for each field
148if isequal(size(InputTable),2)
149[S,D]=get_coeff(XmlData,X,Y);
150end
151
152 %% MAIN LOOP ON FIELDS
153for index=1:NbField
154        update_waitbar(WaitbarHandle,index/NbField)
155    if ~isempty(RUNHandle) && ~strcmp(get(RUNHandle,'BusyAction'),'queue')
156        disp('program stopped by user')
157        return
158    end
159   
160    %%%%%%%%%%%%%%%% loop on views (input lines) %%%%%%%%%%%%%%%%
161    Data=cell(1,NbView);%initiate the set Data
162    timeread=zeros(1,NbView);
163    for iview=1:NbView
164        %% reading input file(s)
165        [Data{iview},tild,errormsg] = read_field(filecell{iview,index},FileType{iview},Param.InputFields,frame_index{iview}(index));
166        if ~isempty(errormsg)
167            disp_uvmat('ERROR',['ERROR in civ2vel_3C/read_field/' errormsg],checkrun)
168            return
169        end
170        % get the time defined in the current file if not already defined from the xml file
171        if ~isempty(time) && isfield(Data{iview},'Time')
172            timeread(iview)=Data{iview}.Time;
173        end
174        if ~isempty(NbSlice_calib)
175            Data{iview}.ZIndex=mod(i1_series{iview}(index)-1,NbSlice_calib{iview})+1;%Zindex for phys transform
176        end
177       
178        %% transform the input field (e.g; phys) if requested (no transform involving two input fields)
179            %camera coordinates
180        xc=R(1)*Xphys+R(2)*Yphys+R(3)*Zphys+Calib.Tx_Ty_Tz(1);
181        yc=R(4)*Xphys+R(5)*Yphys+R(6)*Zphys+Calib.Tx_Ty_Tz(2);
182        zc=R(7)*Xphys+R(8)*Yphys+R(9)*Zphys+Calib.Tx_Ty_Tz(3);
183
184        %undistorted image coordinates
185        Xu=xc./zc;
186        Yu=yc./zc;
187        Data{iview}=phys(Data{iview},XmlData{iview});
188       
189        %% projection on object (gridded plane)
190%         if Param.CheckObject
191            [Data{iview},errormsg]=proj_field(Data{iview},Param.ProjObject);
192            if ~isempty(errormsg)
193                disp_uvmat('ERROR',['ERROR in merge_proge/proj_field: ' errormsg],checkrun)
194                return
195            end
196
197    end
198    %%%%%%%%%%%%%%%% END LOOP ON VIEWS %%%%%%%%%%%%%%%%
199
200    %% merge the NbView fields
201    [MergeData,errormsg]=merge_field(Data);
202    if ~isempty(errormsg)
203        disp_uvmat('ERROR',errormsg,checkrun);
204        return
205    end
206
207    %% time of the merged field: take the average of the different views
208    if ~isempty(time)
209        timeread=time(index);   
210    elseif ~isempty(find(timeread))% time defined from ImaDoc
211        timeread=mean(timeread(timeread~=0));% take average over times form the files (when defined)
212    else
213        timeread=index;% take time=file index
214    end
215
216    %% generating the name of the merged field
217    i1=i1_series{1}(index);
218    if ~isempty(i2_series{end})
219        i2=i2_series{end}(index);
220    else
221        i2=i1;
222    end
223    j1=1;
224    j2=1;
225    if ~isempty(j1_series{1})
226        j1=j1_series{1}(index);
227        if ~isempty(j2_series{end})
228            j2=j2_series{end}(index);
229        else
230            j2=j1;
231        end
232    end
233    OutputFile=fullfile_uvmat(RootPath{1},OutputDir,RootFileOut,FileExtOut,NomType{1},i1,i2,j1,j2);
234
235    %% recording the merged field
236   
237        MergeData.ListGlobalAttribute={'Conventions','Project','InputFile_1','InputFile_end','nb_coord','nb_dim'};
238        MergeData.Conventions='uvmat';
239        MergeData.nb_coord=2;
240        MergeData.nb_dim=2;
241        dt=[];
242        if isfield(Data{1},'dt')&& isnumeric(Data{1}.dt)
243            dt=Data{1}.dt;
244        end
245        for iview =2:numel(Data)
246            if ~(isfield(Data{iview},'dt')&& isequal(Data{iview}.dt,dt))
247                dt=[];%dt not the same for all fields
248            end
249        end
250        if ~isempty(timeread)
251            MergeData.ListGlobalAttribute=[MergeData.ListGlobalAttribute {'Time'}];
252            MergeData.Time=timeread;
253        end
254        if ~isempty(dt)
255            MergeData.ListGlobalAttribute=[MergeData.ListGlobalAttribute {'dt'}];
256            MergeData.dt=dt;
257        end
258        error=struct2nc(OutputFile,MergeData);%save result file
259        if isempty(error)
260            disp(['output file ' OutputFile ' written'])
261        else
262            disp(error)
263        end
264end
265
266% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
267% %read the velocity fields
268% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
269%read field A
270[Field,VelTypeOut]=read_civxdata(file_A,[],vel_type);
271%removes false vectors
272if isfield(Field,'FF')
273    Field.X=Field.X(find(Field.FF==0));
274    Field.Y=Field.Y(find(Field.FF==0));
275    Field.U=Field.U(find(Field.FF==0));
276    Field.V=Field.V(find(Field.FF==0));
277end
278%interpolate on the grid common to both images in phys coordinates
279dXa= griddata_uvmat(Field.X,Field.Y,Field.U,XimaA,YimaA);
280dYa= griddata_uvmat(Field.X,Field.Y,Field.V,XimaA,YimaA);
281dt=Field.dt;
282time=Field.Time;
283
284%read field B
285[Field,VelTypeOut]=read_civxdata(file_B,[],vel_type);
286if ~isequal(Field.dt,dt)
287    msgbox_uvmat('ERROR','different time intervals for the two velocity fields ')
288     return
289end
290if ~isequal(Field.Time,time)
291    msgbox_uvmat('ERROR','different times for the two velocity fields ')
292     return
293end
294%removes false vectors
295if isfield(Field,'FF')
296Field.X=Field.X(find(Field.FF==0));
297Field.Y=Field.Y(find(Field.FF==0));
298Field.U=Field.U(find(Field.FF==0));
299Field.V=Field.V(find(Field.FF==0));
300end
301%interpolate on XimaB
302dXb=griddata_uvmat(Field.X,Field.Y,Field.U,XimaB,YimaB);
303dYb=griddata_uvmat(Field.X,Field.Y,Field.V,XimaB,YimaB);
304%eliminate Not-a-Number
305ind_Nan=find(and(~isnan(dXa),~isnan(dXb)));
306dXa=dXa(ind_Nan);
307dYa=dYa(ind_Nan);
308dXb=dXb(ind_Nan);
309dYb=dYb(ind_Nan);
310grid_phys1(:,1)=grid_real_x(ind_Nan);
311grid_phys1(:,2)=grid_real_y(ind_Nan);
312 
313%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
314%compute the differential coefficients of the geometric calibration
315%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
316[A11,A12,A13,A21,A22,A23]=pxcm_tsai(tsaiA,grid_phys1);
317[B11,B12,B13,B21,B22,B23]=pxcm_tsai(tsaiB,grid_phys1);
318
319C1=A11.*A22-A12.*A21;
320C2=A13.*A22-A12.*A23;
321C3=A13.*A21-A11.*A23;
322D1=B11.*B22-B12.*B21;
323D2=B13.*B22-B12.*B23;
324D3=B13.*B21-B11.*B23;
325A1=(A22.*D1.*(C1.*D3-C3.*D1)+A21.*D1.*(C2.*D1-C1.*D2));
326A2=(A12.*D1.*(C3.*D1-C1.*D3)+A11.*D1.*(C1.*D2-C2.*D1));
327B1=(B22.*C1.*(C3.*D1-C1.*D3)+B21.*C1.*(C1.*D2-C2.*D1));
328B2=(B12.*C1.*(C1.*D3-C3.*D1)+B11.*C1.*(C2.*D1-C1.*D2));
329Lambda=(A1.*dXa+A2.*dYa+B1.*dXb+B2.*dYb)./(A1.*A1+A2.*A2+B1.*B1+B2.*B2);
330
331%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
332%Projection for compatible displacements
333%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
334Ua=dXa-Lambda.*A1;
335Va=dYa-Lambda.*A2;
336Ub=dXb-Lambda.*B1;
337Vb=dYb-Lambda.*B2;
338
339%%%%%%%%%%%%%%%%%%%%%%%%%%%%
340%Calculations of displacements and error
341%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
342U=(A22.*D2.*Ua-A12.*D2.*Va-B22.*C2.*Ub+B12.*C2.*Vb)./(C1.*D2-C2.*D1);
343V=(A21.*D3.*Ua-A11.*D3.*Va-B21.*C3.*Ub+B11.*C3.*Vb)./(C3.*D1-C1.*D3);
344W=(A22.*D1.*Ua-A12.*D1.*Va-B22.*C1.*Ub+B12.*C1.*Vb)./(C2.*D1-C1.*D2);
345W1=(-A21.*D1.*Ua+A11.*D1.*Va+B21.*C1.*Ub-B11.*C1.*Vb)./(C1.*D3-C3.*D1);
346
347error=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));
348
349ind_error=(find(error<thresh_patch));
350U=U(ind_error);
351V=V(ind_error);
352W=W(ind_error);%correction for water interface
353error=error(ind_error);
354
355%create nc grid file
356Result.ListGlobalAttribute={'nb_coord','nb_dim','constant_pixcm','absolut_time_T0','hart','dt','civ'};
357Result.nb_coord=3;%grid file, no velocity
358Result.nb_dim=2;
359Result.constant_pixcm=0;%no linear correspondance with images
360Result.absolut_time_T0=time;%absolute time of the field
361Result.hart=0;
362Result.dt=dt;%time interval for image correlation (put  by default)
363% cte.title='grid';
364Result.civ=0;%not a civ file (no direct correspondance with an image)
365% Result.ListDimName={'nb_vectors'}
366% Result.DimValue=length(U);
367Result.ListVarName={'vec_X','vec_Y','vec_U','vec_V','vec_W','vec_E'};
368Result.VarDimName={'nb_vectors','nb_vectors','nb_vectors','nb_vectors','nb_vectors','nb_vectors'}
369Result.vec_X= grid_phys1(ind_error,1);
370Result.vec_Y= grid_phys1(ind_error,2);
371Result.vec_U=U/dt;
372Result.vec_V=V/dt;
373Result.vec_W=W/dt;
374Result.vec_E=error;
375% error=write_netcdf(file_st,cte,fieldlabels,grid_phys);
376error=struct2nc(file_st,Result);
377display([file_st ' written'])
378
379function [S,D]=get_coeff(XmlData,X,Y);
380Calib_a=XmlData{1}.GeometryCalib;
381R=(Calib_a.R)';%rotation matrix
382A(1,1)=R(1)-R(7)*X;
383A(1,2)=R(2)-R(8)*X;
384A(1,3)=R(3)-R(9)*X;
385A(2,1)=R(4)-R(7)*Y;
386A(2,2)=R(5)-R(8)*Y;
387A(2,3)=R(6)-R(9)*Y;
388
389%'pxcm_tsai': find differentials of the Tsai calibration
390function [A11,A12,A13,A21,A22,A23]=pxcm_tsai(a,var_phys)
391R=(a.R)';
392
393x=var_phys(:,1);
394y=var_phys(:,2);
395
396if isfield(a,'PlanePos')
397    prompt={'Plane 1 Index','Plane 2 Index'};
398    Rep=inputdlg(prompt,'Target displacement test');
399    Z1=str2double(Rep(1));
400    Z2=str2double(Rep(2));
401    z=(a.PlanePos(Z2,3)+a.PlanePos(Z1,3))/2
402else
403    z=0;
404end
405
406%transform coeff for differentiels
407a.C11=R(1)*R(8)-R(2)*R(7);
408a.C12=R(2)*R(7)-R(1)*R(8);
409a.C21=R(4)*R(8)-R(5)*R(7);
410a.C22=R(5)*R(7)-R(4)*R(8);
411a.C1x=R(3)*R(7)-R(9)*R(1);
412a.C1y=R(3)*R(8)-R(9)*R(2);
413a.C2x=R(6)*R(7)-R(9)*R(4);
414a.C2y=R(6)*R(8)-R(9)*R(5);
415
416% %dependence in x,y
417% denom=(R(7)*x+R(8)*y+R(9)*z+a.Tz).*(R(7)*x+R(8)*y+R(9)*z+a.Tz);
418% A11=(a.f*a.sx*(a.C11*y-a.C1x*z+R(1)*a.Tz-R(7)*a.Tx)./denom)/a.dpx;
419% A12=(a.f*a.sx*(a.C12*x-a.C1y*z+R(2)*a.Tz-R(8)*a.Tx)./denom)/a.dpx;
420% A21=(a.f*a.sx*(a.C21*y-a.C2x*z+R(4)*a.Tz-R(7)*a.Ty)./denom)/a.dpy;
421% A22=(a.f*(a.C22*x-a.C2y*z+R(5)*a.Tz-R(8)*a.Ty)./denom)/a.dpy;
422% A13=(a.f*(a.C1x*x+a.C1y*y+R(3)*a.Tz-R(9)*a.Tx)./denom)/a.dpx;
423% A23=(a.f*(a.C2x*x+a.C2y*y+R(6)*a.Tz-R(9)*a.Ty)./denom)/a.dpy;
424
425%dependence in x,y
426denom=(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));
427A11=(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);
428A12=(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);
429A21=(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);
430A22=(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);
431A13=(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);
432A23=(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);
433
434
435
436
437
438
439
440
Note: See TracBrowser for help on using the repository browser.