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

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

fgunctions under series adapted to the renovated GUI series

File size: 20.2 KB
Line 
1%'aver_stat': calculate field average, used with series.fig
2%------------------------------------------------------------------------
3% function GUI_input=aver_stat(Param)
4%
5%OUTPUT
6% GUI_input=list of options in the GUI series.fig needed for the function
7%
8%INPUT:
9% Param: structure containing all the parameters read on the GUI series
10%  or name of the xml file containing these parameters (BATCH case)
11%
12function GUI_input=aver_stat(Param)
13%----------------------------------------------------------------------
14% --- make average on a series of files
15%----------------------------------------------------------------------
16%INPUT:
17%i1_series: series of first indices i (given from the series interface as first_i:incr_i:last_i, mode and list_pair_civ)
18%num_i2: series of second indices i (given from the series interface as first_i:incr_i:last_i, mode and list_pair_civ)
19%j1_series: series of first indices j (given from the series interface as first_j:incr_j:last_j, mode and list_pair_civ )
20%num_j2: series of second indices j (given from the series interface as first_j:incr_j:last_j, mode and list_pair_civ)
21%OTHER INPUTS given by the structure Series
22%  Series.Time:
23%  Series.GeometryCalib:%requests for the visibility of input windows in the GUI series  (activated directly by the selection in the menu ACTION)
24if ~exist('Param','var')
25    GUI_input={'RootPath';'two';...%nbre of possible input series (options 'on'/'two'/'many', default:'one')
26        'SubDir';'on';... % subdirectory of derived files (PIV fields), ('on' by default)
27        'RootFile';'on';... %root input file name ('on' by default)
28        'FileExt';'on';... %input file extension ('on' by default)
29        'NomType';'on';...%type of file indexing ('on' by default)
30        'NbSlice';'on'; ...%nbre of slices ('off' by default)
31        'VelTypeMenu';'two';...% menu for selecting the velocity type (options 'off'/'one'/'two',  'off' by default)
32        'FieldMenu';'two';...% menu for selecting the field (s) in the input file(options 'off'/'one'/'two', 'off' by default)
33        'CoordType'; 'on';...%can use a transform function
34        'GetObject';'on';...%can use projection object(option 'off'/'one'/'two',
35        %'GetMask';'on'...%can use mask option   
36        %'PARAMETER'; %options: name of the user defined parameter',repeat a line for each parameter
37               ''};
38        return
39end
40
41%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
42
43%% input parameters
44% read the xml file for batch case
45if ischar(Param) && ~isempty(find(regexp('Param','.xml$')))
46    Param=xml2struct(Param);
47else %  RUN case: parameters introduced as the input structure Param
48    hseries=guidata(Param.hseries);%handles of the GUI series
49    WaitbarPos=get(hseries.waitbar_frame,'Position');
50end
51[filecell,i1_series,i2_series,j1_series,j2_series]=get_file_series(Param);
52
53%% coordinate transform or other user defined transform (TODO: case BATCH ?)
54transform_fct='';%default
55if isfield(Param,'transform_fct') % transform function handle
56    transform_fct=Param.transform_fct;
57end
58
59%% projection object
60test_object=get(hseries.GetObject,'Value');
61if test_object%isfield(Series,'sethandles')
62    hset_object=findobj(allchild(0),'tag','set_object');
63    ProjObject=read_set_object(guidata(hset_object));
64    %answeryes=questdlg({['field series projected on ' Series.ProjObject.Style]});
65    answeryes=msgbox_uvmat('INPUT_Y-N',['field series projected on ' ProjObject.Style ' before averaging']);
66    if ~isequal(answeryes,'Yes')
67        return
68    end
69end
70
71%% root input file and type
72    RootPath=Param.InputTable(:,1);
73    RootFile=Param.InputTable(:,3);
74    SubDir=Param.InputTable(:,2);
75    NomType=Param.InputTable(:,4);
76    FileExt=Param.InputTable(:,5);
77ext=FileExt{1};
78form=imformats(ext([2:end]));%test valid Matlab image formats
79testima=0;
80if ~isempty(form)||isequal(lower(ext),'.avi')||isequal(lower(ext),'.vol')
81    testima(1)=1;
82end
83if length(FileExt)>=2
84    ext_1=FileExt{2};
85    form=imformats(ext_1([2:end]));%test valid Matlab image formats
86    if ~isempty(form)||isequal(lower(ext_1),'.avi')||isequal(lower(ext_1),'.vol')
87        testima(2)=1;
88    end
89    if testima(2)~=testima(1)
90        msgbox_uvmat('ERROR','images and netcdf files cannot be compared')
91        return
92    end
93end
94
95
96%% Number of input series: this function  accepts two input file series at most (then it operates on the difference of fields)
97nbview=length(RootPath);
98if nbview>2 
99    RootPath=RootPath(1:2);
100    set(hseries.RootPath,'String',RootPath)
101    SubDir=SubDir(1:2);
102    set(hseries.SubDir,'String',SubDir)
103    RootFile=RootFile(1:2);
104    set(hseries.RootFile,'String',RootFile)
105    NomType=NomType(1:2);
106    FileExt=FileExt(1:2);
107    set(hseries.FileExt,'String',FileExt)
108    nbview=2;
109end
110
111%% determine image type
112hhh=which('mmreader');
113for iview=1:nbview
114    if isequal(FileExt{iview},'.nc')||isequal(FileExt{iview},'.cdf')
115        FileType{iview}='netcdf';
116    elseif isequal(lower(FileExt{iview}),'.avi')
117        if ~isequal(hhh,'')&& mmreader.isPlatformSupported()
118            MovieObject{iview}=mmreader(fullfile(RootPath{iview},[RootFile{iview} FileExt{iview}]));
119            FileType{iview}='movie';
120        else
121            FileType{iview}='avi';
122        end
123    elseif isequal(lower(FileExt{iview}),'.vol')
124        FileType{iview}='vol';
125    else
126       form=imformats(FileExt{iview}(2:end));
127       if ~isempty(form)% if the extension corresponds to an image format recognized by Matlab
128           if isequal(NomType{iview},'*');
129               FileType{iview}='multimage';
130           else
131               FileType{iview}='image';
132           end
133       end
134    end
135end
136
137%% number of slices
138NbSlice=Param.IndexRange.NbSlice;
139NbSlice_name=num2str(NbSlice);
140
141%% Field and velocity type (the same for the two views)
142Field_str=get(hseries.FieldMenu,'String');
143FieldName=[]; %default
144testfield=get(hseries.FieldMenu,'Visible');
145if isequal(testfield,'on')
146    val=get(hseries.FieldMenu,'Value');
147    FieldName=Field_str(val);%the same set of fields for all views
148    if isequal(FieldName,{'get_field...'})
149        hget_field=findobj(allchild(0),'name','get_field');%find the get_field... GUI
150        if length(hget_field)>1
151            delete(hget_field(2:end))
152        elseif isempty(hget_field)
153           filename=...
154                 name_generator(fullfile(RootPath{1},RootFile{1}),i1_series{1}(1),j1_series{1}(1),FileExt{1},NomType{1},1,i2_series{1}(1),num_j2{1}(1),SubDir{1});
155           get_field(filename);
156           return
157        end
158        %hhget_field=guidata(hget_field);%handles of GUI elements in get_field
159        SubField=read_get_field(hget_field); %read the names of the variables to plot in the get_field GUI
160    end
161end
162
163%% get the velocity type
164testcivx=0;
165% FileExt=get(hseries.FileExt,'String');
166if ~isequal(FieldName,{'get_field...'})
167    testcivx=isequal(FileType{1},'netcdf');
168end
169if testcivx
170    VelType_str=get(hseries.VelTypeMenu,'String');
171    VelType_val=get(hseries.VelTypeMenu,'Value');
172    VelType{1}=VelType_str{VelType_val};
173    if nbview==2
174        VelType_str=get(hseries.VelTypeMenu_1,'String');
175        VelType_val=get(hseries.VelTypeMenu_1,'Value');
176        VelType{2}=VelType_str{VelType_val};
177    end
178end
179
180%% Calibration data and timing: read the ImaDoc files
181mode=''; %default
182timecell={};
183itime=0;
184NbSlice_calib={};
185for iview=1:nbview%Loop on views
186    XmlData{iview}=[];%default
187    filebase{iview}=fullfile(RootPath{iview},RootFile{iview});
188    if exist([filebase{iview} '.xml'],'file')
189        [XmlData{iview},error]=imadoc2struct([filebase{iview} '.xml']);
190        if isfield(XmlData{iview},'Time')
191            itime=itime+1;
192            timecell{itime}=XmlData{iview}.Time;
193        end
194        if isfield(XmlData{iview},'GeometryCalib') && isfield(XmlData{iview}.GeometryCalib,'SliceCoord')
195            NbSlice_calib{iview}=size(XmlData{iview}.GeometryCalib.SliceCoord,1);%nbre of slices for Zindex in phys transform
196            if ~isequal(NbSlice_calib{iview},NbSlice_calib{1})
197                msgbox_uvmat('WARNING','inconsistent number of Z indices for the two field series');
198            end
199        end
200    elseif exist([filebase{iview} '.civ'],'file')
201        [error,time,TimeUnit,mode,npx,npy,pxcmx,pxcmy]=read_imatext([filebase{iview} '.civ']);
202        itime=itime+1;
203        timecell{itime}=time;
204        XmlData{iview}.Time=time;
205        GeometryCalib.R=[pxcmx 0 0; 0 pxcmy 0;0 0 0];
206        GeometryCalib.Tx=0;
207        GeometryCalib.Ty=0;
208        GeometryCalib.Tz=1;
209        GeometryCalib.dpx=1;
210        GeometryCalib.dpy=1;
211        GeometryCalib.sx=1;
212        GeometryCalib.Cx=0;
213        GeometryCalib.Cy=0;
214        GeometryCalib.f=1;
215        GeometryCalib.kappa1=0;
216        GeometryCalib.CoordUnit='cm';
217        XmlData{iview}.GeometryCalib=GeometryCalib;
218        if error==1
219            msgbox_uvmat('WARNING','inconsistent number of fields in the .civ file');
220        end
221    end
222end
223
224%% check coincidence in time
225multitime=0;
226if isempty(timecell)
227    time=[];
228elseif length(timecell)==1
229    time=timecell{1};
230elseif length(timecell)>1
231    multitime=1;
232    for icell=1:length(timecell)
233        if ~isequal(size(timecell{icell}),size(timecell{1}))
234            msgbox_uvmat('WARNING','inconsistent time array dimensions in ImaDoc fields, the time for the first series is used')
235            time=timecell{1};
236            multitime=0;
237            break
238        end
239    end
240end
241if multitime
242    for icell=1:length(timecell)
243        time(icell,:,:)=timecell{icell};
244    end
245    diff_time=max(max(diff(time)));
246    if diff_time>0
247        msgbox_uvmat('WARNING',['times of series differ by more than ' num2str(diff_time)])
248    end   
249end
250if size(time,2) < i2_series{1}(end) || size(time,3) < num_j2{1}(end)% ime array absent or too short in ImaDoc xml file'
251    time=[];
252end
253
254%% Name(s) of output file(s)
255filebase_out=filebase{1};% the result file has the same root name as the input file series (and the first one is chosen in case of two input series)
256%file extension of the result 
257if testima %case of images
258    ext_out='.png';
259else
260    ext_out='.nc';
261end
262subdir_result='aver_stat';%subdirectory for the results
263pathdir=fullfile(RootPath{1},subdir_result);% full subdirectory name, including path
264testexist=1;
265while testexist
266    pathdir=fullfile(RootPath{1},subdir_result);% full subdirectory name, including path
267    if NbSlice==1% keep track of the first and lsat indices of the input files
268        %NomTypeOut=nomtype2pair(Param.InputTable{1,4},i2_series{end}(end)-i1_series{1}(1),j2_series{end}(end)-j1_series{1}(1));
269        NomTypeOut='_1-2';
270%         RootPath,SubDir,RootFile,FileExt,NomType,i1,i2,j1,j2)
271        fileresult{1}=fullfile_uvmat(RootPath{1},subdir_result,RootFile{1},ext_out,NomTypeOut,i1_series{1}(1),[],[],[]);
272%         fileresult{1}=name_generator(filebase_out,i1_series{1}(1),j1_series{1}(1),ext_out,NomTypeOut,1,i2_series{end}(end),j2_series{end}(end),subdir_result);
273        testexist=exist(fileresult{1},'file');
274    else % simplified indexing with i_slice for multiple slices
275        testexist=0;
276        for i_slice=1:NbSlice
277            fileresult{1}=fullfile_uvmat(RootPath{1},subdir_result,RootFile{1},ext_out,NomTypeOut,i_slice,[],[],[]);
278            %fileresult{i_slice}=name_generator(filebase_out,i_slice,[],ext_out,'_1',1,i_slice,[],subdir_result);
279            if exist(fileresult{i_slice},'file')
280                testexist=1;
281                break
282            end
283        end
284    end
285    if testexist
286        subdir_result=[subdir_result '.0'];
287    end
288end
289% create result directory if needed
290if ~exist(pathdir,'dir')
291    [m1,m2,m3]=mkdir(pathdir);
292    if ~isequal(m2,'')
293        msgbox_uvmat('CONFIRMATION',m2);%error message for directory creation
294    end
295end
296[xx,msg2] = fileattrib(pathdir,'+w','g'); %yield writing access (+w) to user group (g)
297if ~strcmp(msg2,'')
298    msgbox_uvmat('ERROR',['pb of permission for ' pathdir ': ' msg2])%error message for writting access
299    return
300end
301
302%% coordinate transform or other user defined transform
303transform_fct='';%default
304if isfield(Param,'FieldTransform')&&isfield(Param.FieldTransform,'transform_fct')
305    transform_fct=Param.FieldTransform.transform_fct;
306end
307
308%% main loop
309siz=size(i1_series{1});
310nbfield2=siz(1); %nb of consecutive fields at each level(burst
311nbfield=siz(1)*siz(2);
312nbfield=floor(nbfield/(nbfield2*NbSlice));%total number of i indexes (adjusted to an integer number of slices)
313
314% loop on slices
315for i_slice=1:NbSlice
316    for ifield=1:nbfield
317         indselect(:,ifield)=((ifield-1)*NbSlice+(i_slice-1))*nbfield2+[1:nbfield2]';%selected indices on the list of files of a slice
318    end
319    S=0; %initiate the image sum S
320    nbfiles=0;
321    nbmissing=0;
322    % averaging loop
323    for index=1:nbfield*nbfield2
324 %       stopstate=get(hseries.RUN,'BusyAction');
325 %       if isequal(stopstate,'queue') % enable STOP command
326         %   update_waitbar(hseries.waitbar,WaitbarPos,index/(nbfield*nbfield2))
327            ifile=indselect(index);
328            % reading input file(s)
329            for iview=1:nbview
330%                 [filename]=...
331%                     name_generator(filebase{iview},num_i1{iview}(ifile),num_j1{iview}(ifile),FileExt{iview},NomType{iview},1,num_i2{iview}(ifile),j2_series{iview}(ifile),SubDir{iview});
332                    filename=filecell{iview,index};
333                    if ~isequal(FileType{iview},'netcdf')
334                    Data{iview}.ListVarName={'A'};
335                    Data{iview}.AName='image';
336                    switch FileType{iview}
337                        case 'movie'
338                            A=read(MovieObject{iview},i1_series{iview}(ifile));
339                        case 'avi'
340                            mov=aviread(filename,i1_series{iview}(ifile));
341                            A=frame2im(mov(1));
342                        case 'vol'
343                            A=imread(filename);
344                        case 'multimage'
345                            A=imread(filename,i1_series{iview}(ifile));
346                        case 'image'
347                            A=imread(filename);
348                    end
349                    Data{iview}.ListVarName={'AY','AX','A'}; %
350                    Atype{iview}=class(A);
351                    npy=size(A,1);
352                    npx=size(A,2);
353                    nbcolor=size(A,3);
354                    if nbcolor==3
355                        Data{iview}.VarDimName={'AY','AX',{'AY','AX','rgb'}};
356                    else
357                        Data{iview}.VarDimName={'AY','AX',{'AY','AX'}};
358                    end
359                    Data{iview}.AY=[npy-0.5 0.5];
360                    Data{iview}.AX=[0.5 npx-0.5];
361                    Data{iview}.A=double(A);
362                    Data{iview}.CoordUnit='pixel';
363                elseif testcivx
364                    [Data{iview},VelTypeOut]=read_civxdata(filename,FieldName,VelType);
365                else
366                    [Data{iview},var_detect]=nc2struct(filename,SubField.ListVarName); %read the corresponding input data
367                    Data{iview}.VarAttribute=SubField.VarAttribute;
368                end
369                if isfield(Data{iview},'Txt')
370                    msgbox_uvmat('ERROR',['error of input reading: ' Data{iview}.Txt])
371                    return
372                end
373            end
374           
375            % coordinate transform (or other user defined transform)
376            if ~isempty(transform_fct)
377                if ~isempty(NbSlice_calib)
378                    Data{iview}.ZIndex=mod(i1_series{iview}(ifile)-1,NbSlice_calib{1})+1;%Zindex for phys transform
379                end
380                if nbview==2
381                    [Data{1},Data{2}]=transform_fct(Data{1},XmlData{1},Data{2},XmlData{2});
382                    if isempty(Data{2})
383                        Data(2)=[];
384                    end
385                else
386                    Data{1}=transform_fct(Data{1},XmlData{1});
387                end
388            end
389           
390            % field calculation (vort, div...)
391            if testcivx
392                Data{iview}=calc_field(FieldName,Data{iview});%calculate field (vort..)
393            end
394           
395            % field substration (for two input file series)
396            if length(Data)==2
397                [Field,errormsg]=sub_field(Data{1},Data{2}); %substract the two fields
398                if ~isempty(errormsg)
399                    msgbox_uvmat('ERROR',['error in aver_stat/sub_field:' errormsg])
400                    return
401                end
402            else
403                Field=Data{1};
404            end
405            if test_object
406                [Field,errormsg]=proj_field(Field,ProjObject);
407                if ~isempty(errormsg)
408                    msgbox_uvmat('ERROR',['error in aver_stat/proj_field:' errormsg])
409                    return
410                end
411            end
412            nbfiles=nbfiles+1;
413            if nbfiles==1 %first field
414                time_1=[];
415                if isfield(Field,'Time')
416                    time_1=Field.Time(1);
417                end
418                DataMean=Field;%default
419            else
420                for ivar=1:length(Field.ListVarName)
421                    VarName=Field.ListVarName{ivar};
422                    eval(['sizmean=size(DataMean.' VarName ');']);
423                    eval(['siz=size(Field.' VarName ');']);
424                    if ~isequal(siz,sizmean)
425                        msgbox_uvmat('ERROR',['unequal size of input field ' VarName ', need to project  on a grid'])
426                        return
427                    else
428                        eval(['DataMean.' VarName '=DataMean.' VarName '+ Field.' VarName ';']); % update the sum
429                    end
430                end
431            end
432%         end
433    end
434    %end averaging loop
435    for ivar=1:length(Field.ListVarName)
436        VarName=Field.ListVarName{ivar};
437        eval(['DataMean.' VarName '=DataMean.' VarName '/nbfiles;']); % normalize the mean
438    end
439    if nbmissing~=0
440        msgbox_uvmat('WARNING',[num2str(nbmissing) ' input files are missing or skipted'])
441    end
442    if isempty(time) % time read from files  prevails
443        time_end=[];
444        if isfield(Field,'Time')
445            time_end=Field.Time(1);%last time read
446            if ~isempty(time_1)
447                DataMean.Time=time_1;
448                DataMean.Time_end=time_end;
449            end
450        end
451    else  % time from ImaDoc prevails
452        DataMean.Time=time(1,i1_series{1}(1),j1_series{1}(1));
453        DataMean.Time_end=time(end,i1_series{end}(end),j1_series{end}(end));
454    end
455   
456    %writing the result file
457    if testima %case of images
458        if isequal(Atype{1},'uint16')
459            imwrite(uint16(DataMean.A),fileresult{i_slice},'BitDepth',16); % case of 16 bit images
460        else
461            imwrite(uint8(DataMean.A),fileresult{i_slice},'BitDepth',8); % case of 8 bit images
462        end
463        display([fileresult{i_slice} ' written']);
464    else %case of netcdf input file , determine global attributes
465        if isempty(strcmp('Conventions',DataMean.ListGlobalAttribute))
466            DataMean.ListGlobalAttribute=['Conventions' DataMean.ListGlobalAttribute];
467        end
468        DataMean.Conventions='uvmat';
469        DataMean.ListGlobalAttribute=[DataMean.ListGlobalAttribute {Param.Action}];
470        ActionKey='Action';
471        while isfield(DataMean,ActionKey)
472            ActionKey=[ActionKey '_1'];
473        end
474        DataMean.(ActionKey)=Param.Action;
475        DataMean.ListGlobalAttribute=[DataMean.ListGlobalAttribute {ActionKey}];
476        if isfield(DataMean,'Time')
477            DataMean.ListGlobalAttribute=[DataMean.ListGlobalAttribute {'Time','Time_end'}];
478        end
479        errormsg=struct2nc(fileresult{i_slice},DataMean); %save result file
480        if isempty(errormsg)
481            display([fileresult{i_slice} ' written']);
482        else
483            msgbox_uvmat('ERROR',['error in writting result file: ' errormsg])
484            display(errormsg)
485        end
486    end  % end averaging  loop
487end % end loop on slices
488
489%% reproduce ImaDoc/GeometryCalib for image series
490if isfield(XmlData{1},'GeometryCalib') && ~isempty(XmlData{1}.GeometryCalib)
491    [pp,RootFile]=fileparts(filebase_out);
492    outputxml=fullfile(pathdir,[RootFile '.xml'])
493    errormsg=update_imadoc(XmlData{1}.GeometryCalib,outputxml);% introduce the calibration data in the xml file
494    if strcmp(errormsg,'')
495        display(['GeometryCalib transferred to ' outputxml])
496    else
497        msgbox_uvmat('ERROR',errormsg);
498    end
499end
500
501%% open the result file with uvmat
502hget_field=findobj(allchild(0),'name','get_field');%find the get_field... GUI
503delete(hget_field)
504uvmat(fileresult{end})
Note: See TracBrowser for help on using the repository browser.