source: trunk/src/series/aver_stat.m @ 449

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

bug corrected for BATCH mode

File size: 17.0 KB
Line 
1%'aver_stat': calculate field average, used with series.fig
2%------------------------------------------------------------------------
3% function GUI_input=aver_stat(Param)
4%
5%%%%%%%%%%% GENERAL TO ALL SERIES ACTION FCTS %%%%%%%%%%%%%%%%%%%%%%%%%%%
6%OUTPUT
7% GUI_input=list of options in the GUI series.fig needed for the function
8%
9%INPUT:
10% In run mode, the input parameters are given as a Matlab structure Param copied from the GUI series.
11% In batch mode, Param is the name of the corresponding xml file containing the same information
12% In the absence of input (as activated when the current Action is selected
13% in series), the function ouput GUI_input set the activation of the needed GUI elements
14%
15% Param contains the elements:(use the menu bar command 'export/GUI config' in series to 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%    .OutputDir: directory for data outputs, including path
20%    .Action: .ActionName: name of the current activated function
21%             .ActionPath:   path of the current activated function
22%    .IndexRange: set the file or frame indices on which the action must be performed
23%    .FieldTransform: .TransformName: name of the selected transform function
24%                     .TransformPath:   path  of the selected transform function
25%                     .TransformHandle: corresponding function handle
26%    .InputFields: sub structure describing the input fields withfields
27%              .FieldName: name of the field
28%              .VelType: velocity type
29%              .FieldName_1: name of the second field in case of two input series
30%              .VelType_1: velocity type of the second field in case of two input series
31%    .ProjObject: %sub structure describing a projection object (read from ancillary GUI set_object)
32%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
33
34function GUI_input=aver_stat(Param)
35
36%% set the input elements needed on the GUI series when the action is selected in the menu ActionName
37if ~exist('Param','var') % case with no input parameter
38    GUI_input={'NbViewMax';2;...% max nbre of input file series (default='' , no limitation)
39        'AllowInputSort';'off';...% allow alphabetic sorting of the list of input files (options 'off'/'on', 'off' by default)
40        'NbSlice';'on'; ...%nbre of slices ('off' by default)
41        'VelType';'two';...% menu for selecting the velocity type (options 'off'/'one'/'two',  'off' by default)
42        'FieldName';'two';...% menu for selecting the field (s) in the input file(options 'off'/'one'/'two', 'off' by default)
43        'FieldTransform'; 'on';...%can use a transform function
44        'ProjObject';'on';...%can use projection object(option 'off'/'on',
45        'Mask';'off';...%can use mask option   (option 'off'/'on', 'off' by default)
46        'OutputDirExt';'.stat';...%set the output dir extension
47               ''};
48        return
49end
50
51%%%%%%%%%%%% STANDARD PART (DO NOT EDIT) %%%%%%%%%%%%
52%% get input parameters, file names and indices
53% BATCH  case: read the xml file for batch case
54if ischar(Param) && ~isempty(find(regexp(Param,'.xml$')))
55    Param=xml2struct(Param);
56    checkrun=0;
57% RUN case: parameters introduced as the input structure Param 
58else
59    hseries=guidata(Param.hseries);%handles of the GUI series
60    WaitbarPos=get(hseries.waitbar_frame,'Position');%position of the waitbar on the GUI series
61    checkrun=1; % indicate the RUN option is used
62end
63% get the set of input file names (cell array filecell), and the lists of
64% input file or frame indices i1_series,i2_series,j1_series,j2_series
65[filecell,i1_series,i2_series,j1_series,j2_series]=get_file_series(Param);
66% filecell{iview,fileindex}: cell array representing the list of file names
67%        iview: line in the table corresponding to a given file series
68%        fileindex: file index within  the file series,
69% 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
70% i1_series(iview,fileindex) expresses the same indices as a 1D array in file indices
71% set of frame indices used for movie or multimage input
72if ~isempty(j1_series)
73    frame_index=j1_series;
74else
75    frame_index=i1_series;
76end
77
78%% root input file(s) and type
79RootPath=Param.InputTable(:,1);
80RootFile=Param.InputTable(:,3);
81SubDir=Param.InputTable(:,2);
82NomType=Param.InputTable(:,4);
83FileExt=Param.InputTable(:,5);
84
85% numbers of slices and file indices
86NbSlice=1;%default
87if isfield(Param.IndexRange,'NbSlice')
88    NbSlice=Param.IndexRange.NbSlice;
89end
90nbview=size(i1_series,1);%number of input file series (lines in InputTable)
91nbfield_j=size(i1_series,2); %nb of consecutive fields at each level(burst
92nbfield=nbfield_j*size(i1_series,3); %total number of files or frames
93nbfield_i=floor(nbfield/NbSlice);%total number of i indexes (adjusted to an integer number of slices)
94nbfield=nbfield_i*nbfield_j; %total number of fields after adjustement
95
96%determine the file type on each line from the first input file
97ImageTypeOptions={'image','multimage','mmreader','video'};
98NcTypeOptions={'netcdf','civx','civdata'};
99for iview=1:nbview
100    [FileType{iview},FileInfo{iview},Object{iview}]=get_file_type(filecell{iview,1});
101    CheckImage{iview}=~isempty(find(strcmp(FileType{iview},ImageTypeOptions)));% =1 for images
102    CheckNc{iview}=~isempty(find(strcmp(FileType{iview},NcTypeOptions)));% =1 for netcdf files
103end
104
105%% Calibration data and timing: read the ImaDoc files
106mode=''; %default
107timecell={};
108itime=0;
109NbSlice_calib={};
110for iview=1:nbview%Loop on views
111    XmlData{iview}=[];%default
112    filebase{iview}=fullfile(RootPath{iview},RootFile{iview});
113    if exist([filebase{iview} '.xml'],'file')
114        [XmlData{iview},error]=imadoc2struct([filebase{iview} '.xml']);
115        if isfield(XmlData{iview},'Time')
116            itime=itime+1;
117            timecell{itime}=XmlData{iview}.Time;
118        end
119        if isfield(XmlData{iview},'GeometryCalib') && isfield(XmlData{iview}.GeometryCalib,'SliceCoord')
120            NbSlice_calib{iview}=size(XmlData{iview}.GeometryCalib.SliceCoord,1);%nbre of slices for Zindex in phys transform
121            if ~isequal(NbSlice_calib{iview},NbSlice_calib{1})
122                msgbox_uvmat('WARNING','inconsistent number of Z indices for the two field series');
123            end
124        end
125    elseif exist([filebase{iview} '.civ'],'file')
126        [error,time,TimeUnit,mode,npx,npy,pxcmx,pxcmy]=read_imatext([filebase{iview} '.civ']);
127        itime=itime+1;
128        timecell{itime}=time;
129        XmlData{iview}.Time=time;
130        GeometryCalib.R=[pxcmx 0 0; 0 pxcmy 0;0 0 0];
131        GeometryCalib.Tx=0;
132        GeometryCalib.Ty=0;
133        GeometryCalib.Tz=1;
134        GeometryCalib.dpx=1;
135        GeometryCalib.dpy=1;
136        GeometryCalib.sx=1;
137        GeometryCalib.Cx=0;
138        GeometryCalib.Cy=0;
139        GeometryCalib.f=1;
140        GeometryCalib.kappa1=0;
141        GeometryCalib.CoordUnit='cm';
142        XmlData{iview}.GeometryCalib=GeometryCalib;
143        if error==1
144            msgbox_uvmat('WARNING','inconsistent number of fields in the .civ file');
145        end
146    end
147end
148
149%% check coincidence in time for several input file series
150multitime=0;
151if isempty(timecell)
152    time=[];
153elseif length(timecell)==1
154    time=timecell{1};
155elseif length(timecell)>1
156    multitime=1;
157    for icell=1:length(timecell)
158        if ~isequal(size(timecell{icell}),size(timecell{1}))
159            msgbox_uvmat('WARNING','inconsistent time array dimensions in ImaDoc fields, the time for the first series is used')
160            time=timecell{1};
161            multitime=0;
162            break
163        end
164    end
165end
166if multitime
167    for icell=1:length(timecell)
168        time(icell,:,:)=timecell{icell};
169    end
170    diff_time=max(max(diff(time)));
171    if diff_time>0
172        msgbox_uvmat('WARNING',['times of series differ by (max) ' num2str(diff_time)])
173    end   
174end
175if size(time,2) < i2_series(1,end) || size(time,3) < j2_series(1,end)% time array absent or too short in ImaDoc xml file'
176    time=[];
177end
178
179%% coordinate transform or other user defined transform
180transform_fct='';%default
181if isfield(Param,'FieldTransform')&&isfield(Param.FieldTransform,'TransformHandle')
182    transform_fct=Param.FieldTransform.TransformHandle;
183end
184%%%%%%%%%%%% END STANDARD PART  %%%%%%%%%%%%
185 % EDIT FROM HERE
186
187%% check the validity of  input file types
188if CheckImage{1}
189    FileExtOut='.png'; % write result as .png images for image inputs
190elseif CheckNc{1}
191    FileExtOut='.nc';% write result as .nc files for netcdf inputs
192else
193    msgbox_uvmat('ERROR',['invalid file type input ' FileType{1}])
194    return
195end
196if nbview==2 && ~isequal(CheckImage{1},CheckImage{2})
197        msgbox_uvmat('ERROR','input must be two image series or two netcdf file series')
198    return
199end
200NomTypeOut='_1-2_1';% output file index will indicate the first and last ref index in the series
201if NbSlice~=nbfield_j
202    answer=msgbox_uvmat('INPUT_Y-N',['will not average slice by slice: for so cancel and set NbSlice= ' num2str(nbfield_j)]);
203    if ~strcmp(answer,'Yes')
204        return
205    end
206end
207
208%% Set field names and velocity types
209InputFields{1}=[];%default (case of images)
210if isfield(Param,'InputFields')
211    InputFields{1}=Param.InputFields;
212end
213if nbview==2
214    InputFields{2}=[];%default (case of images)
215    if isfield(Param,'InputFields')
216        InputFields{2}=Param.InputFields{1};%default
217        if isfield(Param.InputFields,'FieldName_1')
218            InputFields{2}.FieldName=Param.InputFields.FieldName_1;
219            if isfield(Param.InputFields,'VelType_1')
220                InputFields{2}.VelType=Param.InputFields.VelType_1;
221            end
222        end
223    end
224end
225
226%% Initiate output fields
227%initiate the output structure as a copy of the first input one (reproduce fields)
228[DataOut,ParamOut,errormsg] = read_field(filecell{1,1},FileType{1},InputFields{1},1);
229if ~isempty(errormsg)
230    msgbox_uvmat('ERROR',['error reading ' filecell{1,1} ': ' errormsg])
231    return
232end
233time_1=[];
234if isfield(DataOut,'Time')
235    time_1=DataOut.Time(1);
236end
237if CheckNc{iview}
238    if isempty(strcmp('Conventions',DataOut.ListGlobalAttribute))
239        DataOut.ListGlobalAttribute=['Conventions' DataOut.ListGlobalAttribute];
240    end
241    DataOut.Conventions='uvmat';
242    DataOut.ListGlobalAttribute=[DataOut.ListGlobalAttribute {Param.Action}];
243    ActionKey='Action';
244    while isfield(DataOut,ActionKey)
245        ActionKey=[ActionKey '_1'];
246    end
247    DataOut.(ActionKey)=Param.Action;
248    DataOut.ListGlobalAttribute=[DataOut.ListGlobalAttribute {ActionKey}];
249    if isfield(DataOut,'Time')
250        DataOut.ListGlobalAttribute=[DataOut.ListGlobalAttribute {'Time','Time_end'}];
251    end
252end
253
254%% MAIN LOOP ON SLICES
255%%%%%%%%%%%%% STANDARD PART (DO NOT EDIT) %%%%%%%%%%%%
256for i_slice=1:NbSlice
257    index_slice=i_slice:NbSlice:nbfield;% select file indices of the slice
258    nbfiles=0;
259    nbmissing=0;
260   
261   %initiate result fields
262   for ivar=1:length(DataOut.ListVarName)
263       DataOut.(DataOut.ListVarName{ivar})=0; % initialise all fields to zero
264   end
265
266    %%%%%%%%%%%%%%%% loop on field indices %%%%%%%%%%%%%%%%
267    for index=index_slice
268        if checkrun
269            update_waitbar(hseries.waitbar_frame,WaitbarPos,index/(nbfield))
270            stopstate=get(hseries.RUN,'BusyAction');
271        else
272            stopstate='queue';
273        end
274       
275        %%%%%%%%%%%%%%%% loop on views (input lines) %%%%%%%%%%%%%%%%
276        for iview=1:nbview
277            % reading input file(s)
278            [Data{iview},ParamOut,errormsg] = read_field(filecell{iview,index},FileType{iview},InputFields{iview},frame_index(iview,index));
279            if ~isempty(errormsg)
280                errormsg=['error of input reading: ' errormsg];
281                break
282            end
283            if ~isempty(NbSlice_calib)
284                Data{iview}.ZIndex=mod(i1_series(iview,index)-1,NbSlice_calib{1})+1;%Zindex for phys transform
285            end
286        end
287        Field=[]; % initiate the current input field structure
288        %%%%%%%%%%%%%%%% end loop on views (input lines) %%%%%%%%%%%%%%%%
289        %%%%%%%%%%%% END STANDARD PART  %%%%%%%%%%%%
290        % EDIT FROM HERE
291
292        if isempty(errormsg)     
293            % coordinate transform (or other user defined transform)
294            if ~isempty(transform_fct)
295                if nbview==2
296                    [Data{1},Data{2}]=transform_fct(Data{1},XmlData{1},Data{2},XmlData{2});
297                    if isempty(Data{2})
298                        Data(2)=[];
299                    end
300                else
301                    Data{1}=transform_fct(Data{1},XmlData{1});
302                end
303            end
304           
305            % field calculation (vort, div...)
306            if strcmp(FileType{1},'civx')||strcmp(FileType{1},'civ')
307                Data{1}=calc_field(InputFields{1}.FieldName,Data{1});%calculate field (vort..)
308            end
309           
310            % field substration (for two input file series)
311            if length(Data)==2
312                if strcmp(FileType{2},'civx')||strcmp(FileType{2},'civ')
313                    Data{2}=calc_field(InputFields{2}.FieldName,Data{2});%calculate field (vort..)
314                end
315                [Field,errormsg]=sub_field(Data{1},Data{2}); %substract the two fields
316                if ~isempty(errormsg)
317                    msgbox_uvmat('ERROR',['error in aver_stat/sub_field:' errormsg])
318                    return
319                end
320            else
321                Field=Data{1};
322            end
323            if Param.CheckObject
324                [Field,errormsg]=proj_field(Field,ProjObject);
325                if ~isempty(errormsg)
326                    msgbox_uvmat('ERROR',['error in aver_stat/proj_field:' errormsg])
327                    return
328                end
329            end
330            nbfiles=nbfiles+1;
331           
332            %%%%%%%%%%%% MAIN RUNNING OPERATIONS  %%%%%%%%%%%%
333            %update sum
334            for ivar=1:length(Field.ListVarName)
335                VarName=Field.ListVarName{ivar};
336                sizmean=size(DataOut.(VarName));
337                siz=size(Field.(VarName));
338                if ~isequal(DataOut.(VarName),0)&& ~isequal(siz,sizmean)
339                    msgbox_uvmat('ERROR',['unequal size of input field ' VarName ', need to project  on a grid'])
340                    return
341                else
342                    DataOut.(VarName)=DataOut.(VarName)+ double(Field.(VarName)); % update the sum
343                end
344            end
345            %%%%%%%%%%%%   END MAIN RUNNING OPERATIONS  %%%%%%%%%%%%
346        else
347            display(errormsg) 
348        end
349    end
350    %%%%%%%%%%%%%%%% end loop on field indices %%%%%%%%%%%%%%%%
351   
352    for ivar=1:length(Field.ListVarName)
353        VarName=Field.ListVarName{ivar};
354        DataOut.(VarName)=DataOut.(VarName)/nbfiles; % normalize the mean
355    end
356    if nbmissing~=0
357        msgbox_uvmat('WARNING',[num2str(nbmissing) ' input files are missing or skipted'])
358    end
359    if isempty(time) % time read from files  prevails
360        if isfield(Field,'Time')
361            time_end=Field.Time(1);%last time read
362            if ~isempty(time_1)
363                DataOut.Time=time_1;
364                DataOut.Time_end=time_end;
365            end
366        end
367    else  % time from ImaDoc prevails
368        DataOut.Time=time(1,i1_series(1,1),j1_series(1,1));
369        DataOut.Time_end=time(end,i1_series(end,end),j1_series(end,end));
370    end
371   
372    %writing the result file
373    OutputFile=fullfile_uvmat(RootPath{1},Param.OutputSubDir,RootFile{1},FileExtOut,NomTypeOut,i1_series(1,1),i1_series(1,end),i_slice,[]);
374    if CheckImage{1} %case of images
375        if isequal(FileInfo{1}.BitDepth,16)||(numel(FileInfo)==2 &&isequal(FileInfo{2}.BitDepth,16))
376            DataOut.A=uint16(DataOut.A);
377            imwrite(DataOut.A,OutputFile,'BitDepth',16); % case of 16 bit images
378        else
379            DataOut.A=uint8(DataOut.A);
380            imwrite(DataOut.A,OutputFile,'BitDepth',8); % case of 16 bit images
381        end
382        display([OutputFile ' written']);
383    else %case of netcdf input file , determine global attributes
384        errormsg=struct2nc(OutputFile,DataOut); %save result file
385        if isempty(errormsg)
386            display([OutputFile ' written']);
387        else
388            msgbox_uvmat('ERROR',['error in writting result file: ' errormsg])
389            display(errormsg)
390        end
391    end  % end averaging  loop
392end
393%%%%%%%%%%%%%%%% end loop on slices %%%%%%%%%%%%%%%%
394
395%% reproduce ImaDoc/GeometryCalib for image series
396% if isfield(XmlData{1},'GeometryCalib') && ~isempty(XmlData{1}.GeometryCalib)
397%     [tild,RootFile]=fileparts(filebase_out);
398%     outputxml=fullfile(pathdir,[RootFile '.xml']);
399%     errormsg=update_imadoc(XmlData{1}.GeometryCalib,outputxml);% introduce the calibration data in the xml file
400%     if strcmp(errormsg,'')
401%         display(['GeometryCalib transferred to ' outputxml])
402%     else
403%         msgbox_uvmat('ERROR',errormsg);
404%     end
405% end
406
407%% open the result file with uvmat
408hget_field=findobj(allchild(0),'name','get_field');%find the get_field... GUI
409delete(hget_field)
410uvmat(OutputFile)% open the last result file with uvmat
Note: See TracBrowser for help on using the repository browser.