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

Last change on this file since 862 was 862, checked in by sommeria, 6 years ago

bugs repaired

File size: 17.5 KB
Line 
1%'civ2vel_3C': combine velocity fields from two cameras to 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     disp_uvmat('ERROR','no geometric calibration available for image A',checkrun)
119     return
120 end
121 if isfield(XmlData{2},'GeometryCalib')
122     tsaiB=XmlData{2}.GeometryCalib;
123 else
124     disp_uvmat('ERROR','no geometric calibration available for image B',checkrun)
125     return
126 end
127[filecell,i1_series,i2_series,j1_series,j2_series]=get_file_series(Param);
128
129%% grid of physical positions (given by projection plane)
130if ~Param.CheckObject
131         disp_uvmat('ERROR','a projection plane with interpolation is needed',checkrun)
132     return
133end
134ObjectData=Param.ProjObject;
135
136[x,y]=meshgrid(ObjectData.RangeX(1):ObjectData.DX:ObjectData.RangeX(2),ObjectData.RangeY(1):ObjectData.DY:ObjectData.RangeY(2));
137z=zeros(size(x));
138[Field,ParamOut,errormsg] = read_field(FileName,FileType,ParamIn,num)
139%camera coordinates:initialisation 2 cameras
140NbCamera=2;
141X=zeros(NbCamera,size(x,1),size(y,1));
142Y=zeros(NbCamera,size(x,1),size(y,1));
143for icamera=1:NbCamera
144%camera coordinates
145Calib=XmlData{1}.GeometryCalib;
146    xc=R(1)*Xphys+R(2)*Yphys+R(3)*Zphys+Calib.Tx_Ty_Tz(1);
147    yc=R(4)*Xphys+R(5)*Yphys+R(6)*Zphys+Calib.Tx_Ty_Tz(2);
148    zc=R(7)*Xphys+R(8)*Yphys+R(9)*Zphys+Calib.Tx_Ty_Tz(3);
149   
150    %undistorted image coordinates
151    X(icamera,:,:)=xc./zc;
152    Y(icamera,:,:)=yc./zc;
153end
154
155%% case of fixed planes (for moving interface, coeff need to be calculated for each field
156if isequal(size(InputTable),2)
157[S,D]=get_coeff(XmlData,X,Y);
158end
159
160 %% MAIN LOOP ON FIELDS
161for index=1:NbField
162        update_waitbar(WaitbarHandle,index/NbField)
163    if ~isempty(RUNHandle) && ~strcmp(get(RUNHandle,'BusyAction'),'queue')
164        disp('program stopped by user')
165        return
166    end
167   
168    %%%%%%%%%%%%%%%% loop on views (input lines) %%%%%%%%%%%%%%%%
169    Data=cell(1,NbView);%initiate the set Data
170    timeread=zeros(1,NbView);
171    for iview=1:NbView
172        %% reading input file(s)
173        [Data{iview},tild,errormsg] = read_field(filecell{iview,index},FileType{iview},Param.InputFields,frame_index{iview}(index));
174        if ~isempty(errormsg)
175            disp_uvmat('ERROR',['ERROR in civ2vel_3C/read_field/' errormsg],checkrun)
176            return
177        end
178        % get the time defined in the current file if not already defined from the xml file
179        if ~isempty(time) && isfield(Data{iview},'Time')
180            timeread(iview)=Data{iview}.Time;
181        end
182        if ~isempty(NbSlice_calib)
183            Data{iview}.ZIndex=mod(i1_series{iview}(index)-1,NbSlice_calib{iview})+1;%Zindex for phys transform
184        end
185       
186        %% transform the input field (e.g; phys) if requested (no transform involving two input fields)
187            %camera coordinates
188        Data{iview}=phys(Data{iview},XmlData{iview});
189       
190        %% projection on object (gridded plane)
191%         if Param.CheckObject
192            [Data{iview},errormsg]=proj_field(Data{iview},Param.ProjObject);
193            if ~isempty(errormsg)
194                disp_uvmat('ERROR',['ERROR in merge_proge/proj_field: ' errormsg],checkrun)
195                return
196            end
197
198    end
199    %%%%%%%%%%%%%%%% END LOOP ON VIEWS %%%%%%%%%%%%%%%%
200
201    %% merge the NbView fields
202    [MergeData,errormsg]=merge_field(Data);
203    if ~isempty(errormsg)
204        disp_uvmat('ERROR',errormsg,checkrun);
205        return
206    end
207
208    %% time of the merged field: take the average of the different views
209    if ~isempty(time)
210        timeread=time(index);   
211    elseif ~isempty(find(timeread))% time defined from ImaDoc
212        timeread=mean(timeread(timeread~=0));% take average over times form the files (when defined)
213    else
214        timeread=index;% take time=file index
215    end
216
217    %% generating the name of the merged field
218    i1=i1_series{1}(index);
219    if ~isempty(i2_series{end})
220        i2=i2_series{end}(index);
221    else
222        i2=i1;
223    end
224    j1=1;
225    j2=1;
226    if ~isempty(j1_series{1})
227        j1=j1_series{1}(index);
228        if ~isempty(j2_series{end})
229            j2=j2_series{end}(index);
230        else
231            j2=j1;
232        end
233    end
234    OutputFile=fullfile_uvmat(RootPath{1},OutputDir,RootFileOut,FileExtOut,NomType{1},i1,i2,j1,j2);
235
236    %% recording the merged field
237   
238        MergeData.ListGlobalAttribute={'Conventions','Project','InputFile_1','InputFile_end','nb_coord','nb_dim'};
239        MergeData.Conventions='uvmat';
240        MergeData.nb_coord=2;
241        MergeData.nb_dim=2;
242        dt=[];
243        if isfield(Data{1},'dt')&& isnumeric(Data{1}.dt)
244            dt=Data{1}.dt;
245        end
246        for iview =2:numel(Data)
247            if ~(isfield(Data{iview},'dt')&& isequal(Data{iview}.dt,dt))
248                dt=[];%dt not the same for all fields
249            end
250        end
251        if ~isempty(timeread)
252            MergeData.ListGlobalAttribute=[MergeData.ListGlobalAttribute {'Time'}];
253            MergeData.Time=timeread;
254        end
255        if ~isempty(dt)
256            MergeData.ListGlobalAttribute=[MergeData.ListGlobalAttribute {'dt'}];
257            MergeData.dt=dt;
258        end
259        error=struct2nc(OutputFile,MergeData);%save result file
260        if isempty(error)
261            disp(['output file ' OutputFile ' written'])
262        else
263            disp(error)
264        end
265end
266
267% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
268% %read the velocity fields
269% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
270%read field A
271[Field,VelTypeOut]=read_civxdata(file_A,[],vel_type);
272%removes false vectors
273if isfield(Field,'FF')
274    Field.X=Field.X(find(Field.FF==0));
275    Field.Y=Field.Y(find(Field.FF==0));
276    Field.U=Field.U(find(Field.FF==0));
277    Field.V=Field.V(find(Field.FF==0));
278end
279%interpolate on the grid common to both images in phys coordinates
280dXa= griddata_uvmat(Field.X,Field.Y,Field.U,XimaA,YimaA);
281dYa= griddata_uvmat(Field.X,Field.Y,Field.V,XimaA,YimaA);
282dt=Field.dt;
283time=Field.Time;
284
285%read field B
286[Field,VelTypeOut]=read_civxdata(file_B,[],vel_type);
287if ~isequal(Field.dt,dt)
288    msgbox_uvmat('ERROR','different time intervals for the two velocity fields ')
289     return
290end
291if ~isequal(Field.Time,time)
292    msgbox_uvmat('ERROR','different times for the two velocity fields ')
293     return
294end
295%removes false vectors
296if isfield(Field,'FF')
297Field.X=Field.X(find(Field.FF==0));
298Field.Y=Field.Y(find(Field.FF==0));
299Field.U=Field.U(find(Field.FF==0));
300Field.V=Field.V(find(Field.FF==0));
301end
302%interpolate on XimaB
303dXb=griddata_uvmat(Field.X,Field.Y,Field.U,XimaB,YimaB);
304dYb=griddata_uvmat(Field.X,Field.Y,Field.V,XimaB,YimaB);
305%eliminate Not-a-Number
306ind_Nan=find(and(~isnan(dXa),~isnan(dXb)));
307dXa=dXa(ind_Nan);
308dYa=dYa(ind_Nan);
309dXb=dXb(ind_Nan);
310dYb=dYb(ind_Nan);
311grid_phys1(:,1)=grid_real_x(ind_Nan);
312grid_phys1(:,2)=grid_real_y(ind_Nan);
313 
314%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
315%compute the differential coefficients of the geometric calibration
316%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
317[A11,A12,A13,A21,A22,A23]=pxcm_tsai(tsaiA,grid_phys1);
318[B11,B12,B13,B21,B22,B23]=pxcm_tsai(tsaiB,grid_phys1);
319
320C1=A11.*A22-A12.*A21;
321C2=A13.*A22-A12.*A23;
322C3=A13.*A21-A11.*A23;
323D1=B11.*B22-B12.*B21;
324D2=B13.*B22-B12.*B23;
325D3=B13.*B21-B11.*B23;
326A1=(A22.*D1.*(C1.*D3-C3.*D1)+A21.*D1.*(C2.*D1-C1.*D2));
327A2=(A12.*D1.*(C3.*D1-C1.*D3)+A11.*D1.*(C1.*D2-C2.*D1));
328B1=(B22.*C1.*(C3.*D1-C1.*D3)+B21.*C1.*(C1.*D2-C2.*D1));
329B2=(B12.*C1.*(C1.*D3-C3.*D1)+B11.*C1.*(C2.*D1-C1.*D2));
330Lambda=(A1.*dXa+A2.*dYa+B1.*dXb+B2.*dYb)./(A1.*A1+A2.*A2+B1.*B1+B2.*B2);
331
332%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
333%Projection for compatible displacements
334%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
335Ua=dXa-Lambda.*A1;
336Va=dYa-Lambda.*A2;
337Ub=dXb-Lambda.*B1;
338Vb=dYb-Lambda.*B2;
339
340%%%%%%%%%%%%%%%%%%%%%%%%%%%%
341%Calculations of displacements and error
342%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
343U=(A22.*D2.*Ua-A12.*D2.*Va-B22.*C2.*Ub+B12.*C2.*Vb)./(C1.*D2-C2.*D1);
344V=(A21.*D3.*Ua-A11.*D3.*Va-B21.*C3.*Ub+B11.*C3.*Vb)./(C3.*D1-C1.*D3);
345W=(A22.*D1.*Ua-A12.*D1.*Va-B22.*C1.*Ub+B12.*C1.*Vb)./(C2.*D1-C1.*D2);
346W1=(-A21.*D1.*Ua+A11.*D1.*Va+B21.*C1.*Ub-B11.*C1.*Vb)./(C1.*D3-C3.*D1);
347
348error=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));
349
350ind_error=(find(error<thresh_patch));
351U=U(ind_error);
352V=V(ind_error);
353W=W(ind_error);%correction for water interface
354error=error(ind_error);
355
356%create nc grid file
357Result.ListGlobalAttribute={'nb_coord','nb_dim','constant_pixcm','absolut_time_T0','hart','dt','civ'};
358Result.nb_coord=3;%grid file, no velocity
359Result.nb_dim=2;
360Result.constant_pixcm=0;%no linear correspondance with images
361Result.absolut_time_T0=time;%absolute time of the field
362Result.hart=0;
363Result.dt=dt;%time interval for image correlation (put  by default)
364% cte.title='grid';
365Result.civ=0;%not a civ file (no direct correspondance with an image)
366% Result.ListDimName={'nb_vectors'}
367% Result.DimValue=length(U);
368Result.ListVarName={'vec_X','vec_Y','vec_U','vec_V','vec_W','vec_E'};
369Result.VarDimName={'nb_vectors','nb_vectors','nb_vectors','nb_vectors','nb_vectors','nb_vectors'}
370Result.vec_X= grid_phys1(ind_error,1);
371Result.vec_Y= grid_phys1(ind_error,2);
372Result.vec_U=U/dt;
373Result.vec_V=V/dt;
374Result.vec_W=W/dt;
375Result.vec_E=error;
376% error=write_netcdf(file_st,cte,fieldlabels,grid_phys);
377error=struct2nc(file_st,Result);
378display([file_st ' written'])
379
380function [S,D]=get_coeff(XmlData,X,Y);
381Calib_a=XmlData{1}.GeometryCalib;
382R=(Calib_a.R)';%rotation matrix
383A(1,1)=R(1)-R(7)*X;
384A(1,2)=R(2)-R(8)*X;
385A(1,3)=R(3)-R(9)*X;
386A(2,1)=R(4)-R(7)*Y;
387A(2,2)=R(5)-R(8)*Y;
388A(2,3)=R(6)-R(9)*Y;
389
390%'pxcm_tsai': find differentials of the Tsai calibration
391function [A11,A12,A13,A21,A22,A23]=pxcm_tsai(a,var_phys)
392R=(a.R)';
393
394x=var_phys(:,1);
395y=var_phys(:,2);
396
397if isfield(a,'PlanePos')
398    prompt={'Plane 1 Index','Plane 2 Index'};
399    Rep=inputdlg(prompt,'Target displacement test');
400    Z1=str2double(Rep(1));
401    Z2=str2double(Rep(2));
402    z=(a.PlanePos(Z2,3)+a.PlanePos(Z1,3))/2
403else
404    z=0;
405end
406
407%transform coeff for differentiels
408a.C11=R(1)*R(8)-R(2)*R(7);
409a.C12=R(2)*R(7)-R(1)*R(8);
410a.C21=R(4)*R(8)-R(5)*R(7);
411a.C22=R(5)*R(7)-R(4)*R(8);
412a.C1x=R(3)*R(7)-R(9)*R(1);
413a.C1y=R(3)*R(8)-R(9)*R(2);
414a.C2x=R(6)*R(7)-R(9)*R(4);
415a.C2y=R(6)*R(8)-R(9)*R(5);
416
417% %dependence in x,y
418% denom=(R(7)*x+R(8)*y+R(9)*z+a.Tz).*(R(7)*x+R(8)*y+R(9)*z+a.Tz);
419% A11=(a.f*a.sx*(a.C11*y-a.C1x*z+R(1)*a.Tz-R(7)*a.Tx)./denom)/a.dpx;
420% A12=(a.f*a.sx*(a.C12*x-a.C1y*z+R(2)*a.Tz-R(8)*a.Tx)./denom)/a.dpx;
421% A21=(a.f*a.sx*(a.C21*y-a.C2x*z+R(4)*a.Tz-R(7)*a.Ty)./denom)/a.dpy;
422% A22=(a.f*(a.C22*x-a.C2y*z+R(5)*a.Tz-R(8)*a.Ty)./denom)/a.dpy;
423% A13=(a.f*(a.C1x*x+a.C1y*y+R(3)*a.Tz-R(9)*a.Tx)./denom)/a.dpx;
424% A23=(a.f*(a.C2x*x+a.C2y*y+R(6)*a.Tz-R(9)*a.Ty)./denom)/a.dpy;
425
426%dependence in x,y
427denom=(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));
428A11=(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);
429A12=(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);
430A21=(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);
431A22=(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);
432A13=(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);
433A23=(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);
434
435
436
437
438
439
440
441
Note: See TracBrowser for help on using the repository browser.