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

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

functions updated in series for the new file configuration

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