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

Last change on this file since 839 was 839, checked in by sommeria, 9 years ago

extract_rdvision corrected, civ2vel_3C introduced (combined vel fields to get 3D components)

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