source: trunk/src/series/merge_proj.m @ 605

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

bugs corrected to get an advancement bar with status

File size: 17.1 KB
Line 
1%'merge_proj': concatene several fields from series, can project them on a regular grid in phys coordinates
2%------------------------------------------------------------------------
3% function ParamOut=merge_proj(Param)
4%------------------------------------------------------------------------
5%%%%%%%%%%% GENERAL TO ALL SERIES ACTION FCTS %%%%%%%%%%%%%%%%%%%%%%%%%%%
6%
7%OUTPUT
8% ParamOut: sets options in the GUI series.fig needed for the function
9%
10%INPUT:
11% In run mode, the input parameters are given as a Matlab structure Param copied from the GUI series.
12% In batch mode, Param is the name of the corresponding xml file containing the same information
13% when Param.Action.RUN=0 (as activated when the current Action is selected
14% in series), the function ouput paramOut set the activation of the needed GUI elements
15%
16% Param contains the elements:(use the menu bar command 'export/GUI config' in series to
17% see the current structure Param)
18%    .InputTable: cell of input file names, (several lines for multiple input)
19%                      each line decomposed as {RootPath,SubDir,Rootfile,NomType,Extension}
20%    .OutputSubDir: name of the subdirectory for data outputs
21%    .OutputDirExt: directory extension for data outputs
22%    .Action: .ActionName: name of the current activated function
23%             .ActionPath:   path of the current activated function
24%             .ActionExt: fct extension ('.m', Matlab fct, '.sh', compiled   Matlab fct
25%             .RUN =0 for GUI input, =1 for function activation
26%             .RunMode='local','background', 'cluster': type of function  use
27%             
28%    .IndexRange: set the file or frame indices on which the action must be performed
29%    .FieldTransform: .TransformName: name of the selected transform function
30%                     .TransformPath:   path  of the selected transform function
31%    .InputFields: sub structure describing the input fields withfields
32%              .FieldName: name(s) of the field
33%              .VelType: velocity type
34%              .FieldName_1: name of the second field in case of two input series
35%              .VelType_1: velocity type of the second field in case of two input series
36%              .Coord_y: name of y coordinate variable
37%              .Coord_x: name of x coordinate variable
38%    .ProjObject: %sub structure describing a projection object (read from ancillary GUI set_object)
39%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
40
41function ParamOut=merge_proj(Param)
42
43%% set the input elements needed on the GUI series when the function is selected in the menu ActionName
44if isstruct(Param) && isequal(Param.Action.RUN,0)
45    ParamOut.AllowInputSort='off';...% allow alphabetic sorting of the list of input file SubDir (options 'off'/'on', 'off' by default)
46    ParamOut.WholeIndexRange='on';...% prescribes the file index ranges from min to max (options 'off'/'on', 'off' by default)
47    ParamOut.NbSlice='off'; ...%nbre of slices ('off' by default)
48    ParamOut.VelType='one';...% menu for selecting the velocity type (options 'off'/'one'/'two',  'off' by default)
49    ParamOut.FieldName='one';...% menu for selecting the field (s) in the input file(options 'off'/'one'/'two', 'off' by default)
50    ParamOut.FieldTransform = 'on';...%can use a transform function
51    ParamOut.ProjObject='on';...%can use projection object(option 'off'/'on',
52    ParamOut.Mask='off';...%can use mask option   (option 'off'/'on', 'off' by default)
53    ParamOut.OutputDirExt='.mproj';%set the output dir extension
54    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
55    filecell=get_file_series(Param);%check existence of the first input file
56    if ~exist(filecell{1,1},'file')
57        msgbox_uvmat('WARNING','the first input file does not exist')
58    elseif isequal(size(Param.InputTable,1),1) && ~isfield(Param,'ProjObject')
59         msgbox_uvmat('WARNING','a projection object of type plane needs to be introduced for merge_proj')
60    end
61return
62end
63
64%%%%%%%%%%%% STANDARD PART (DO NOT EDIT) %%%%%%%%%%%%
65
66
67%% read input parameters from an xml file if input is a file name (batch mode)
68checkrun=1;
69if ischar(Param)
70    Param=xml2struct(Param);% read Param as input file (batch case)
71    checkrun=0;
72end
73
74ParamOut=Param; %default output
75if ~isfield(Param,'InputFields')
76    Param.InputFields.FieldName='';
77end
78OutputSubDir=[Param.OutputSubDir Param.OutputDirExt];% subdirectory for output files
79
80%% root input file type
81RootPath=Param.InputTable(:,1);
82RootFile=Param.InputTable(:,3);
83SubDir=Param.InputTable(:,2);
84NomType=Param.InputTable(:,4);
85FileExt=Param.InputTable(:,5);
86[filecell,i1_series,i2_series,j1_series,j2_series]=get_file_series(Param);
87%%%%%%%%%%%%
88% The cell array filecell is the list of input file names, while
89% filecell{iview,fileindex}:
90%        iview: line in the table corresponding to a given file series
91%        fileindex: file index within  the file series,
92% i1_series(iview,ref_j,ref_i)... are the corresponding arrays of indices i1,i2,j1,j2, depending on the input line iview and the two reference indices ref_i,ref_j
93% i1_series(iview,fileindex) expresses the same indices as a 1D array in file indices
94%%%%%%%%%%%%
95NbSlice=1;%default
96if isfield(Param.IndexRange,'NbSlice')&&~isempty(Param.IndexRange.NbSlice)
97    NbSlice=Param.IndexRange.NbSlice;
98end
99nbview=numel(i1_series);%number of input file series (lines in InputTable)
100nbfield_j=size(i1_series{1},1); %nb of fields for the j index (bursts or volume slices)
101nbfield_i=size(i1_series{1},2); %nb of fields for the i index
102nbfield=nbfield_j*nbfield_i; %total number of fields
103nbfield_i=floor(nbfield/NbSlice);%total number of  indexes in a slice (adjusted to an integer number of slices)
104nbfield=nbfield_i*NbSlice; %total number of fields after adjustement
105
106%determine the file type on each line from the first input file
107ImageTypeOptions={'image','multimage','mmreader','video'};
108NcTypeOptions={'netcdf','civx','civdata'};
109for iview=1:nbview
110    if ~exist(filecell{iview,1}','file')
111        displ_uvmat('ERROR',['the first input file ' filecell{iview,1} ' does not exist'],checkrun)
112        return
113    end
114    [FileType{iview},FileInfo{iview},MovieObject{iview}]=get_file_type(filecell{iview,1});
115    CheckImage{iview}=~isempty(find(strcmp(FileType{iview},ImageTypeOptions)));% =1 for images
116    CheckNc{iview}=~isempty(find(strcmp(FileType{iview},NcTypeOptions)));% =1 for netcdf files
117    if ~isempty(j1_series{iview})
118        frame_index{iview}=j1_series{iview};
119    else
120        frame_index{iview}=i1_series{iview};
121    end
122end
123
124
125%% calibration data and timing: read the ImaDoc files
126[XmlData,NbSlice_calib,time,errormsg]=read_multimadoc(RootPath,SubDir,RootFile,FileExt,i1_series,i2_series,j1_series,j2_series);
127if size(time,1)>1
128    diff_time=max(max(diff(time)));
129    if diff_time>0
130        displ_uvmat('WARNING',['times of series differ by (max) ' num2str(diff_time)],checkrun)
131    end   
132end
133
134%% coordinate transform or other user defined transform
135transform_fct='';%default
136if isfield(Param,'FieldTransform')&&~isempty(Param.FieldTransform.TransformName)
137    addpath(Param.FieldTransform.TransformPath)
138    transform_fct=str2func(Param.FieldTransform.TransformName);
139    rmpath(Param.FieldTransform.TransformPath)
140end
141
142%%%%%%%%%%%% END STANDARD PART  %%%%%%%%%%%%
143 % EDIT FROM HERE
144
145%% check the validity of  input file types
146if CheckImage{1}
147    FileExtOut='.png'; % write result as .png images for image inputs
148elseif CheckNc{1}
149    FileExtOut='.nc';% write result as .nc files for netcdf inputs
150else
151    displ_uvmat('ERROR',['invalid file type input ' FileType{1}],checkrun)
152    return
153end
154for iview=1:nbview
155        if ~isequal(CheckImage{iview},CheckImage{1})||~isequal(CheckNc{iview},CheckNc{1})
156        displ_uvmat('ERROR','input set of input series: need  either netcdf either image series',checkrun)
157    return
158    end
159end
160NomTypeOut=NomType;% output file index will indicate the first and last ref index in the series
161% if checkrun==1
162%     ParamOut.Specific=[];%no specific parameter
163%     return %stop here for interactive input (option Param.Specific='?')
164% end
165
166%% Set field names and velocity types
167%use Param.InputFields for all views
168
169%% MAIN LOOP ON SLICES
170%%%%%%%%%%%%% STANDARD PART (DO NOT EDIT) %%%%%%%%%%%%
171for i_slice=1:NbSlice
172    index_slice=i_slice:NbSlice:nbfield;% select file indices of the slice
173    nbfiles=0;
174    nbmissing=0;
175
176    %%%%%%%%%%%%%%%% loop on field indices %%%%%%%%%%%%%%%%
177    for index=index_slice
178        if checkrun
179            stopstate=get(Param.RUNHandle,'BusyAction');
180            update_waitbar(Param.WaitbarHandle,index/nbfield)
181        else
182            stopstate='queue';
183        end
184        if ~isequal(stopstate,'queue')% enable STOP command
185            return
186        end
187        %%%%%%%%%%%%%%%% loop on views (input lines) %%%%%%%%%%%%%%%%
188        Data=cell(1,nbview);%initiate the set Data
189        nbtime=0;
190        for iview=1:nbview
191            %% reading input file(s)
192            [Data{iview},tild,errormsg] = read_field(filecell{iview,index},FileType{iview},Param.InputFields,frame_index{iview}(index));
193            if ~isempty(errormsg)
194                errormsg=['merge_proj/read_field/' errormsg];
195                display(errormsg)
196                break
197            end
198            timeread(iview)=0;
199            if isfield(Data{iview},'Time')
200                    timeread(iview)=Data{iview}.Time;
201                    nbtime=nbtime+1;
202                end
203            if ~isempty(NbSlice_calib)
204                Data{iview}.ZIndex=mod(i1_series{iview}(index)-1,NbSlice_calib{iview})+1;%Zindex for phys transform
205            end
206           
207            %% transform the input field (e.g; phys) if requested
208            if ~isempty(transform_fct)
209                if nargin(transform_fct)>=2
210                    Data{iview}=transform_fct(Data{iview},XmlData{iview});
211                else
212                    Data{iview}=transform_fct(Data{iview});
213                end
214            end
215           
216            %% check whether tps is needed, then calculate tps coefficients if needed
217            check_tps=0;
218            if isfield(Param.InputFields,'FieldName')
219                if ischar(Param.InputFields.FieldName)
220                    Param.InputFields.FieldName={Param.InputFields.FieldName};
221                end
222            else
223                Param.InputFields.FieldName={};
224            end
225            for ilist=1:numel(Param.InputFields.FieldName)
226                switch Param.InputFields.FieldName{ilist}
227                    case {'vort','div','strain'}
228                        check_tps=1;
229                end
230            end
231               
232            %% calculate tps coeff if needed
233             check_proj_tps= ~isempty(Param.ProjObject)&& strcmp(Param.ProjObject.ProjMode,'filter')&&~isfield(Data{iview},'Coord_tps');
234            Data{iview}=tps_coeff_field(Data{iview},check_proj_tps);
235
236            %% projection on object (gridded plane)
237            if Param.CheckObject
238                [Data{iview},errormsg]=proj_field(Data{iview},Param.ProjObject);
239                if ~isempty(errormsg)
240                    displ_uvmat('ERROR',['error in merge_proge/proj_field: ' errormsg],checkrun)
241                    return
242                end
243            end
244        end
245        %----------END LOOP ON VIEWS----------------------
246       
247        %% merge the nbview fields
248        MergeData=merge_field(Data);
249        if isfield(MergeData,'Txt')
250            displ_uvmat('ERROR',MergeData.Txt,checkrun)
251            return
252        end
253       
254        % time of the merged field:
255        if ~isempty(time)% time defined from ImaDoc
256            timeread=time(:,index);
257        end
258        timeread=mean(timeread);
259       
260        % generating the name of the merged field
261        i1=i1_series{iview}(index);
262        if ~isempty(i2_series{iview})
263            i2=i2_series{iview}(index);
264        else
265            i2=i1;
266        end
267        j1=1;
268        j2=1;
269        if ~isempty(j1_series{iview})
270            j1=j1_series{iview}(index);
271            if ~isempty(j2_series{iview})
272                j2=j2_series{iview}(index);
273            else
274                j2=j1;
275            end
276        end
277        OutputFile=fullfile_uvmat(RootPath{1},OutputSubDir,RootFile{1},FileExtOut,NomType{1},i1,i2,j1,j2);
278       
279        % recording the merged field
280        if CheckImage{1}    %in case of input images an image is produced
281            if isa(MergeData.A,'uint8')
282                bitdepth=8;
283            elseif isa(MergeData.A,'uint16')
284                bitdepth=16;
285            end
286            imwrite(MergeData.A,OutputFile,'BitDepth',bitdepth);
287            %write xml calibration file
288            siz=size(MergeData.A);
289            npy=siz(1);
290            npx=siz(2);
291            if isfield(MergeData,'VarAttribute')&&isfield(MergeData.VarAttribute{1},'Coord_2')&&isfield(MergeData.VarAttribute{1},'Coord_1')
292                Rangx=MergeData.VarAttribute{1}.Coord_2;
293                Rangy=MergeData.VarAttribute{1}.Coord_1;
294            elseif isfield(MergeData,'AX')&& isfield(MergeData,'AY')
295                Rangx=[MergeData.AX(1) MergeData.AX(end)];
296                Rangy=[MergeData.AY(1) MergeData.AY(end)];
297            else
298                Rangx=[0.5 npx-0.5];
299                Rangy=[npy-0.5 0.5];%default
300            end
301            pxcmx=(npx-1)/(Rangx(2)-Rangx(1));
302            pxcmy=(npy-1)/(Rangy(1)-Rangy(2));
303            T_x=-pxcmx*Rangx(1)+0.5;
304            T_y=-pxcmy*Rangy(2)+0.5;
305            GeometryCal.focal=1;
306            GeometryCal.R=[pxcmx,0,0;0,pxcmy,0;0,0,1];
307            GeometryCal.Tx_Ty_Tz=[T_x T_y 1];
308            ImaDoc.GeometryCalib=GeometryCal;
309%             t=struct2xml(ImaDoc);
310%             t=set(t,1,'name','ImaDoc');
311%             save(t,[filebase_merge '.xml'])
312%             display([filebase_merge '.xml saved'])
313        else
314            MergeData.ListGlobalAttribute={'Conventions','Project','InputFile_1','InputFile_end','nb_coord','nb_dim','dt','Time','civ'};
315            MergeData.Conventions='uvmat';
316            MergeData.nb_coord=2;
317            MergeData.nb_dim=2;
318            dt=[];
319            if isfield(Data{1},'dt')&& isnumeric(Data{1}.dt)
320                dt=Data{1}.dt;
321            end
322            for iview =2:numel(Data)
323                if ~(isfield(Data{iview},'dt')&& isequal(Data{iview}.dt,dt))
324                    dt=[];%dt not the same for all fields
325                end
326            end
327            if isempty(dt)
328                MergeData.ListGlobalAttribute(6)=[];
329            else
330                MergeData.dt=dt;
331            end
332            MergeData.Time=timeread;
333            error=struct2nc(OutputFile,MergeData);%save result file
334            if isempty(error)
335                display(['output file ' OutputFile ' written'])
336            else
337                display(error)
338            end
339        end
340    end
341end
342
343%'merge_field': concatene fields
344%------------------------------------------------------------------------
345function MergeData=merge_field(Data)
346%% default output
347if isempty(Data)||~iscell(Data)
348    MergeData=[];
349    return
350end
351MergeData=Data{1};%default
352error=0;
353nbview=length(Data);
354if nbview==1
355    return
356end
357
358%% group the variables (fields of 'Data') in cells of variables with the same dimensions
359[CellVarIndex,NbDim,VarTypeCell]=find_field_cells(Data{1});
360%LOOP ON GROUPS OF VARIABLES SHARING THE SAME DIMENSIONS
361% CellVarIndex=cells of variable index arrays
362for icell=1:length(CellVarIndex)
363    if NbDim(icell)==1
364        continue
365    end
366    VarIndex=CellVarIndex{icell};%  indices of the selected variables in the list FieldData.ListVarName
367    VarType=VarTypeCell{icell};
368    ivar_X=VarType.coord_x;
369    ivar_Y=VarType.coord_y;
370    ivar_FF=VarType.errorflag;
371    if isempty(ivar_X)
372        test_grid=1;%test for input data on regular grid (e.g. image)coordinates
373    else
374        if length(ivar_Y)~=1
375                displ_uvmat('ERROR','y coordinate missing in proj_field.m',checkrun)
376                return
377        end
378        test_grid=0;
379    end
380    %case of input fields with unstructured coordinates
381    if ~test_grid
382        for ivar=VarIndex
383            VarName=MergeData.ListVarName{ivar};
384            for iview=1:nbview
385                MergeData.(VarName)=[MergeData.(VarName); Data{iview}.(VarName)];
386            end
387        end
388    %case of fields defined on a structured  grid
389    else 
390        testFF=0;
391        for iview=2:nbview
392            for ivar=VarIndex
393                VarName=MergeData.ListVarName{ivar};
394                if isfield(MergeData,'VarAttribute')
395                    if length(MergeData.VarAttribute)>=ivar && isfield(MergeData.VarAttribute{ivar},'Role') && isequal(MergeData.VarAttribute{ivar}.Role,'errorflag')
396                        testFF=1;
397                    end
398                end
399                MergeData.(VarName)=MergeData.(VarName) + Data{iview}.(VarName);
400            end
401        end
402        if testFF
403            nbaver=nbview-MergeData.FF;
404            indgood=find(nbaver>0);
405            for ivar=VarIndex
406                VarName=MergeData.ListVarName{ivar};
407                MergeData.(VarName)(indgood)=double(MergeData.(VarName)(indgood))./nbaver(indgood);
408            end
409        else
410            for ivar=VarIndex
411                VarName=MergeData.ListVarName{ivar};
412                MergeData.(VarName)=double(MergeData.(VarName))./nbview;
413            end   
414        end
415    end
416end
417
418   
Note: See TracBrowser for help on using the repository browser.