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

Last change on this file since 504 was 494, checked in by sommeria, 12 years ago

various bugs corrected after testing in Windows OS. Introduction
of filter tps

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