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

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

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

File size: 17.5 KB
Line 
1%'aver_stat': calculate field average, used with series.fig
2% this function can be used as a template for applying a global operation (here averaging) on a series of input fields
3%------------------------------------------------------------------------
4% function ParamOut=aver_stat(Param)
5%
6%%%%%%%%%%% GENERAL TO ALL SERIES ACTION FCTS %%%%%%%%%%%%%%%%%%%%%%%%%%%
7%
8% This function is used in four modes by the GUI series:
9%           1) config GUI: with no input argument, the function determine the suitable GUI configuration
10%           2) interactive input: the function is used to interactively introduce input parameters, and then stops
11%           3) RUN: the function itself runs, when an appropriate input  structure Param has been introduced.
12%           4) BATCH: the function itself proceeds in BATCH mode, using an xml file 'Param' as input.
13%
14% This function is used in four modes by the GUI series:
15%           1) config GUI: with no input argument, the function determine the suitable GUI configuration
16%           2) interactive input: the function is used to interactively introduce input parameters, and then stops
17%           3) RUN: the function itself runs, when an appropriate input  structure Param has been introduced.
18%           4) BATCH: the function itself proceeds in BATCH mode, using an xml file 'Param' as input.
19%
20%OUTPUT
21% GUI_input=list of options in the GUI series.fig needed for the function
22%
23%INPUT:
24% In run mode, the input parameters are given as a Matlab structure Param copied from the GUI series.
25% In batch mode, Param is the name of the corresponding xml file containing the same information
26% In the absence of input (as activated when the current Action is selected
27% in series), the function ouput GUI_input set the activation of the needed GUI elements
28%
29% Param contains the elements:(use the menu bar command 'export/GUI config' in series to see the current structure Param)
30%    .InputTable: cell of input file names, (several lines for multiple input)
31%                      each line decomposed as {RootPath,SubDir,Rootfile,NomType,Extension}
32%    .OutputSubDir: name of the subdirectory for data outputs
33%    .OutputDirExt: directory extension for data outputs
34%    .Action: .ActionName: name of the current activated function
35%             .ActionPath:   path of the current activated function
36%    .IndexRange: set the file or frame indices on which the action must be performed
37%    .FieldTransform: .TransformName: name of the selected transform function
38%                     .TransformPath:   path  of the selected transform function
39%                     .TransformHandle: corresponding function handle
40%    .InputFields: sub structure describing the input fields withfields
41%              .FieldName: name of the field
42%              .VelType: velocity type
43%              .FieldName_1: name of the second field in case of two input series
44%              .VelType_1: velocity type of the second field in case of two input series
45%    .ProjObject: %sub structure describing a projection object (read from ancillary GUI set_object)
46%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
47
48function ParamOut=aver_stat(Param)
49
50%% set the input elements needed on the GUI series when the action is selected in the menu ActionName
51if ~exist('Param','var') % case with no input parameter
52    ParamOut={'AllowInputSort';'off';...% allow alphabetic sorting of the list of input files (options 'off'/'on', 'off' by default)
53        'WholeIndexRange';'off';...% prescribes the file index ranges from min to max (options 'off'/'on', 'off' by default)
54        'NbSlice';'on'; ...%nbre of slices ('off' by default)
55        'VelType';'two';...% menu for selecting the velocity type (options 'off'/'one'/'two',  'off' by default)
56        'FieldName';'two';...% menu for selecting the field (s) in the input file(options 'off'/'one'/'two', 'off' by default)
57        'FieldTransform'; 'on';...%can use a transform function
58        'ProjObject';'on';...%can use projection object(option 'off'/'on',
59        'Mask';'off';...%can use mask option   (option 'off'/'on', 'off' by default)
60        'OutputDirExt';'.stat';...%set the output dir extension
61               ''};
62        return
63end
64
65%%%%%%%%%%%%  STANDARD PART  %%%%%%%%%%%%
66%% select different modes,  RUN, parameter input, BATCH
67% BATCH  case: read the xml file for batch case
68if ischar(Param)
69        Param=xml2struct(Param);
70        checkrun=0;
71% RUN case: parameters introduced as the input structure Param
72else
73    if isfield(Param,'Specific')&& strcmp(Param.Specific,'?')
74        checkrun=1;% will only search interactive input parameters (preparation of BATCH mode)
75    else
76        checkrun=2; % indicate the RUN option is used
77    end
78    hseries=guidata(Param.hseries);%handles of the GUI series
79end
80ParamOut=Param; %default output
81OutputDir=[Param.OutputSubDir Param.OutputDirExt];
82   
83%% root input file(s) and type
84RootPath=Param.InputTable(:,1);
85RootFile=Param.InputTable(:,3);
86SubDir=Param.InputTable(:,2);
87NomType=Param.InputTable(:,4);
88FileExt=Param.InputTable(:,5);
89[filecell,i1_series,i2_series,j1_series,j2_series]=get_file_series(Param);
90%%%%%%%%%%%%
91% The cell array filecell is the list of input file names, while
92% filecell{iview,fileindex}:
93%        iview: line in the table corresponding to a given file series
94%        fileindex: file index within  the file series,
95% 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
96% i1_series(iview,fileindex) expresses the same indices as a 1D array in file indices
97%%%%%%%%%%%%
98NbSlice=1;%default
99if isfield(Param.IndexRange,'NbSlice')&&~isempty(Param.IndexRange.NbSlice)
100    NbSlice=Param.IndexRange.NbSlice;
101end
102nbview=numel(i1_series);%number of input file series (lines in InputTable)
103nbfield_j=size(i1_series{1},1); %nb of fields for the j index (bursts or volume slices)
104nbfield_i=size(i1_series{1},2); %nb of fields for the i index
105nbfield=nbfield_j*nbfield_i; %total number of fields
106nbfield_i=floor(nbfield/NbSlice);%total number of  indexes in a slice (adjusted to an integer number of slices)
107nbfield=nbfield_i*NbSlice; %total number of fields after adjustement
108
109%determine the file type on each line from the first input file
110ImageTypeOptions={'image','multimage','mmreader','video'};
111NcTypeOptions={'netcdf','civx','civdata'};
112for iview=1:nbview
113    if ~exist(filecell{iview,1}','file')
114        msgbox_uvmat('ERROR',['the first input file ' filecell{iview,1} ' does not exist'])
115        return
116    end
117    [FileType{iview},FileInfo{iview},MovieObject{iview}]=get_file_type(filecell{iview,1});
118    CheckImage{iview}=~isempty(find(strcmp(FileType{iview},ImageTypeOptions)));% =1 for images
119    CheckNc{iview}=~isempty(find(strcmp(FileType{iview},NcTypeOptions)));% =1 for netcdf files
120    if ~isempty(j1_series{iview})
121        frame_index{iview}=j1_series{iview};
122    else
123        frame_index{iview}=i1_series{iview};
124    end
125end
126
127%% calibration data and timing: read the ImaDoc files
128[XmlData,NbSlice_calib,time,errormsg]=read_multimadoc(RootPath,SubDir,RootFile,FileExt,i1_series,i2_series,j1_series,j2_series);
129if size(time,1)>1
130    diff_time=max(max(diff(time)));
131    if diff_time>0
132        msgbox_uvmat('WARNING',['times of series differ by (max) ' num2str(diff_time)])
133    end   
134end
135
136%% coordinate transform or other user defined transform
137transform_fct='';%default
138if isfield(Param,'FieldTransform')&&~isempty(Param.FieldTransform.TransformName)
139    addpath(Param.FieldTransform.TransformPath)
140    transform_fct=str2func(Param.FieldTransform.TransformName);
141    rmpath(Param.FieldTransform.TransformPath)
142end
143
144%%%%%%%%%%%% END STANDARD PART  %%%%%%%%%%%%
145 % EDIT FROM HERE
146
147%% check the validity of  input file types
148if CheckImage{1}
149    FileExtOut='.png'; % write result as .png images for image inputs
150elseif CheckNc{1}
151    FileExtOut='.nc';% write result as .nc files for netcdf inputs
152else
153    msgbox_uvmat('ERROR',['invalid file type input ' FileType{1}])
154    return
155end
156if nbview==2 && ~isequal(CheckImage{1},CheckImage{2})
157        msgbox_uvmat('ERROR','input must be two image series or two netcdf file series')
158    return
159end
160NomTypeOut='_1-2_1';% output file index will indicate the first and last ref index in the series
161if checkrun==1
162    return % stop here for input checks
163end
164
165%% Set field names and velocity types
166InputFields{1}=[];%default (case of images)
167if isfield(Param,'InputFields')
168    InputFields{1}=Param.InputFields;
169end
170if nbview==2
171    InputFields{2}=[];%default (case of images)
172    if isfield(Param,'InputFields')
173        InputFields{2}=Param.InputFields{1};%default
174        if isfield(Param.InputFields,'FieldName_1')
175            InputFields{2}.FieldName=Param.InputFields.FieldName_1;
176            if isfield(Param.InputFields,'VelType_1')
177                InputFields{2}.VelType=Param.InputFields.VelType_1;
178            end
179        end
180    end
181end
182
183%% MAIN LOOP ON SLICES
184%%%%%%%%%%%%% STANDARD PART (DO NOT EDIT) %%%%%%%%%%%%
185for i_slice=1:NbSlice
186    index_slice=i_slice:NbSlice:nbfield;% select file indices of the slice
187    nbfiles=0;
188    nbmissing=0;
189
190    %%%%%%%%%%%%%%%% loop on field indices %%%%%%%%%%%%%%%%
191    for index=index_slice
192        if checkrun
193            update_waitbar(hseries.Waitbar,index/(nbfield))
194            stopstate=get(hseries.RUN,'BusyAction');
195        else
196            stopstate='queue';
197        end
198       
199        %%%%%%%%%%%%%%%% loop on views (input lines) %%%%%%%%%%%%%%%%
200        for iview=1:nbview
201            % reading input file(s)
202            [Data{iview},tild,errormsg] = read_field(filecell{iview,index},FileType{iview},InputFields{iview},frame_index{iview}(index));
203            if ~isempty(errormsg)
204                errormsg=['error of input reading: ' errormsg];
205                break
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        end
211        Field=[]; % initiate the current input field structure
212        %%%%%%%%%%%%%%%% end loop on views (input lines) %%%%%%%%%%%%%%%%
213        %%%%%%%%%%%% END STANDARD PART  %%%%%%%%%%%%
214        % EDIT FROM HERE
215       
216        if isempty(errormsg)
217            % coordinate transform (or other user defined transform)
218            if ~isempty(transform_fct)
219                if nbview==2
220                    [Data{1},Data{2}]=transform_fct(Data{1},XmlData{1},Data{2},XmlData{2});
221                    if isempty(Data{2})
222                        Data(2)=[];
223                    end
224                else
225                    Data{1}=transform_fct(Data{1},XmlData{1});
226                end
227            end
228           
229            %% check whether tps is needed, then calculate tps coefficients if needed
230            check_tps=0;
231            if ischar(Param.InputFields.FieldName)
232                Param.InputFields.FieldName={Param.InputFields.FieldName};
233            end
234            for ilist=1:numel(Param.InputFields.FieldName)
235                switch Param.InputFields.FieldName{ilist}
236                    case {'vort','div','strain'}
237                        check_tps=1;
238                end
239            end
240            if strcmp(Param.ProjObject.ProjMode,'filter')
241                check_tps=1;
242            end
243            if check_tps
244                SubDomain=1500; %default, estimated nbre of vectors in a subdomain used for tps
245                if isfield(Data{iview},'SubDomain')
246                    SubDomain=Data{iview}.SubDomain;%
247                end
248                [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] =...
249                    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);
250                nbvar=numel(Data{iview}.ListVarName);
251                Data{iview}.ListVarName=[Data{iview}.ListVarName {'SubRange','NbSites','Coord_tps','U_tps','V_tps'}];
252                Data{iview}.VarDimName=[Data{iview}.VarDimName {{'nb_coord','nb_bounds','nb_subdomain'},{'nb_subdomain'},...
253                    {'nb_tps','nb_coord','nb_subdomain'},{'nb_tps','nb_subdomain'},{'nb_tps','nb_subdomain'}}];
254                Data{iview}.VarAttribute{nbvar+3}.Role='coord_tps';
255                Data{iview}.VarAttribute{nbvar+4}.Role='vector_x';
256                Data{iview}.VarAttribute{nbvar+5}.Role='vector_y';
257                if isfield(Data{iview},'ListDimName')%cleaning
258                    Data{iview}=rmfield(Data{iview},'ListDimName');
259                end
260                if isfield(Data{iview},'DimValue')%cleaning
261                    Data{iview}=rmfield(Data{iview},'DimValue');
262                end
263            end
264                 
265            % field calculation (vort, div...)   
266            if strcmp(FileType{1},'civx')||strcmp(FileType{1},'civdata')
267                if isfield(Data{1},'Coord_tps')
268                    Data{1}.FieldList=Param.InputFields.FieldName;
269                else
270                    Data{1}=calc_field(Param.InputFields.FieldName,Data{1});%calculate field (vort..)
271                end
272            end
273         
274            % field substration (for two input file series)
275            if length(Data)==2
276                if strcmp(FileType{2},'civx')||strcmp(FileType{2},'civdata')
277                    if isfield(Data{2},'Coord_tps')
278                        Data{2}.FieldList=Param.InputFields.FieldName;
279                    else
280                        Data{2}=calc_field(Param.InputFields.FieldName,Data{2});%calculate field (vort..)
281                    end
282                end
283                [Field,errormsg]=sub_field(Data{1},Data{2}); %substract the two fields
284                if ~isempty(errormsg)
285                    msgbox_uvmat('ERROR',['error in aver_stat/sub_field:' errormsg])
286                    return
287                end
288            else
289                Field=Data{1};
290            end
291           
292            %field projection on an object
293            if Param.CheckObject
294                [Field,errormsg]=proj_field(Field,Param.ProjObject);
295                if ~isempty(errormsg)
296                    msgbox_uvmat('ERROR',['error in aver_stat/proj_field:' errormsg])
297                    return
298                end
299            end
300            nbfiles=nbfiles+1;
301           
302            %%%%%%%%%%%% MAIN RUNNING OPERATIONS  %%%%%%%%%%%%
303            %update sum
304            if nbfiles==1 %first field
305                time_1=[];
306                if isfield(Field,'Time')
307                    time_1=Field.Time(1);
308                end
309                DataOut=Field;%default
310                for ivar=1:length(Field.ListVarName)
311                    VarName=Field.ListVarName{ivar};
312                    DataOut.(VarName)=double(DataOut.(VarName));
313                end
314            else   %current field
315                for ivar=1:length(Field.ListVarName)
316                    VarName=Field.ListVarName{ivar};
317                    sizmean=size(DataOut.(VarName));
318                    siz=size(Field.(VarName));
319                    if ~isequal(DataOut.(VarName),0)&& ~isequal(siz,sizmean)
320                        msgbox_uvmat('ERROR',['unequal size of input field ' VarName ', need to project  on a grid'])
321                        return
322                    else
323                        DataOut.(VarName)=DataOut.(VarName)+ double(Field.(VarName)); % update the sum
324                    end
325                end
326            end
327            %%%%%%%%%%%%   END MAIN RUNNING OPERATIONS  %%%%%%%%%%%%
328        else
329            display(errormsg)
330        end
331    end
332    %%%%%%%%%%%%%%%% end loop on field indices %%%%%%%%%%%%%%%%
333   
334    for ivar=1:length(Field.ListVarName)
335        VarName=Field.ListVarName{ivar};
336        DataOut.(VarName)=DataOut.(VarName)/nbfiles; % normalize the mean
337    end
338    if nbmissing~=0
339        msgbox_uvmat('WARNING',[num2str(nbmissing) ' input files are missing or skipted'])
340    end
341    if isempty(time) % time is read from files
342        if isfield(Field,'Time')
343            time_end=Field.Time(1);%last time read
344            if ~isempty(time_1)
345                DataOut.Time=time_1;
346                DataOut.Time_end=time_end;
347            end
348        end
349    else  % time from ImaDoc prevails if it exists
350        j1=1;%default
351        if ~isempty(j1_series{1})
352            j1=j1_series{1};
353        end
354        DataOut.Time=time(1,i1_series{1}(1),j1filexml);
355        DataOut.Time_end=time(end,i1_series{end}(end),j1_series{end}(end));
356    end
357   
358    %writing the result file
359    OutputFile=fullfile_uvmat(RootPath{1},OutputDir,RootFile{1},FileExtOut,NomTypeOut,i1_series{1}(1),i1_series{1}(end),i_slice,[]);
360    if CheckImage{1} %case of images
361        if isequal(FileInfo{1}.BitDepth,16)||(numel(FileInfo)==2 &&isequal(FileInfo{2}.BitDepth,16))
362            DataOut.A=uint16(DataOut.A);
363            imwrite(DataOut.A,OutputFile,'BitDepth',16); % case of 16 bit images
364        else
365            DataOut.A=uint8(DataOut.A);
366            imwrite(DataOut.A,OutputFile,'BitDepth',8); % case of 16 bit images
367        end
368        display([OutputFile ' written']);
369    else %case of netcdf input file , determine global attributes
370        errormsg=struct2nc(OutputFile,DataOut); %save result file
371        if isempty(errormsg)
372            display([OutputFile ' written']);
373        else
374            msgbox_uvmat('ERROR',['error in writting result file: ' errormsg])
375            display(errormsg)
376        end
377    end  % end averaging  loop
378end
379%%%%%%%%%%%%%%%% end loop on slices %%%%%%%%%%%%%%%%
380
381%% open the result file with uvmat (in RUN mode)
382if checkrun
383    hget_field=findobj(allchild(0),'name','get_field');%find the get_field... GUI
384    delete(hget_field)
385    uvmat(OutputFile)% open the last result file with uvmat
386end
Note: See TracBrowser for help on using the repository browser.