[884]  1  %'aver_stat': calculate field average over a time series


 2  %


 3  % function ParamOut=aver_stat(Param)


 4  %


 5  %%%%%%%%%%% GENERAL TO ALL SERIES ACTION FCTS %%%%%%%%%%%%%%%%%%%%%%%%%%%


 6  %


 7  %OUTPUT


 8  % ParamOut: sets options in the GUI series.fig needed for the function


 9  %


 10  %INPUT:


 11  % In run mode, the input parameters are given as a Matlab structure Param copied from the GUI series.


 12  % In batch mode, Param is the name of the corresponding xml file containing the same information


 13  % when Param.Action.RUN=0 (as activated when the current Action is selected


 14  % in series), the function ouput paramOut set the activation of the needed GUI elements


 15  %


 16  % Param contains the elements:(use the menu bar command 'export/GUI config' in series to


 17  % see the current structure Param)


 18  % .InputTable: cell of input file names, (several lines for multiple input)


 19  % each line decomposed as {RootPath,SubDir,Rootfile,NomType,Extension}


 20  % .OutputSubDir: name of the subdirectory for data outputs


 21  % .OutputDirExt: directory extension for data outputs


 22  % .Action: .ActionName: name of the current activated function


 23  % .ActionPath: path of the current activated function


 24  % .ActionExt: fct extension ('.m', Matlab fct, '.sh', compiled Matlab fct


 25  % .RUN =0 for GUI input, =1 for function activation


 26  % .RunMode='local','background', 'cluster': type of function use


 27  %


 28  % .IndexRange: set the file or frame indices on which the action must be performed


 29  % .FieldTransform: .TransformName: name of the selected transform function


 30  % .TransformPath: path of the selected transform function


 31  % .InputFields: sub structure describing the input fields withfields


 32  % .FieldName: name(s) of the field


 33  % .VelType: velocity type


 34  % .FieldName_1: name of the second field in case of two input series


 35  % .VelType_1: velocity type of the second field in case of two input series


 36  % .Coord_y: name of y coordinate variable


 37  % .Coord_x: name of x coordinate variable


 38  % .ProjObject: %sub structure describing a projection object (read from ancillary GUI set_object)


 39  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


 40 


 41  %=======================================================================


 42  % Copyright 20082014, LEGI UMR 5519 / CNRS UJF GINP, Grenoble, France


 43  % http://www.legi.grenobleinp.fr


 44  % Joel.Sommeria  Joel.Sommeria (A) legi.cnrs.fr


 45  %


 46  % This file is part of the toolbox UVMAT.


 47  %


 48  % UVMAT is free software; you can redistribute it and/or modify


 49  % it under the terms of the GNU General Public License as published


 50  % by the Free Software Foundation; either version 2 of the license,


 51  % or (at your option) any later version.


 52  %


 53  % UVMAT is distributed in the hope that it will be useful,


 54  % but WITHOUT ANY WARRANTY; without even the implied warranty of


 55  % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the


 56  % GNU General Public License (see LICENSE.txt) for more details.


 57  %=======================================================================


 58 


 59  function ParamOut=aver_stat(Param)


 60 


 61  %% set the input elements needed on the GUI series when the action is selected in the menu ActionName


 62  if isstruct(Param) && isequal(Param.Action.RUN,0)% function activated from the GUI series but not RUN


 63  ParamOut.AllowInputSort='off';% allow alphabetic sorting of the list of input file SubDir (options 'off'/'on', 'off' by default)


 64  ParamOut.WholeIndexRange='off';% prescribes the file index ranges from min to max (options 'off'/'on', 'off' by default)


 65  ParamOut.NbSlice='on'; %nbre of slices ('off' by default)


 66  ParamOut.VelType='two';% menu for selecting the velocity type (options 'off'/'one'/'two', 'off' by default)


 67  ParamOut.FieldName='two';% menu for selecting the field (s) in the input file(options 'off'/'one'/'two', 'off' by default)


 68  ParamOut.FieldTransform = 'on';%can use a transform function


 69  ParamOut.ProjObject='on';%can use projection object(option 'off'/'on',


 70  ParamOut.Mask='off';%can use mask option (option 'off'/'on', 'off' by default)


 71  ParamOut.OutputDirExt='.stat';%set the output dir extension


 72  ParamOut.OutputFileMode='NbSlice';% '=NbInput': 1 output file per input file index, '=NbInput_i': 1 file per input file index i, '=NbSlice': 1 file per slice


 73  % check for selection of a projection object


 74  hseries=findobj(allchild(0),'Tag','series');% handles of the GUI series


 75  if ~isfield(Param,'ProjObject')


 76  answer=msgbox_uvmat('INPUT_YN','use a projection object?');


 77  if strcmp(answer,'Yes')


 78  hhseries=guidata(hseries);


 79  set(hhseries.CheckObject,'Visible','on')


 80  set(hhseries.CheckObject,'Value',1)


 81  Param.CheckObject=1;


 82  series('CheckObject_Callback',hseries,[],hhseries); %file input with xml reading in uvmat, show the image in phys coordinates


 83  end


 84  end


 85  % introduce bin size for histograms


 86  if isfield(Param,'CheckObject') && Param.CheckObject


 87  SeriesData=get(hseries,'UserData');


 88  if ismember(SeriesData.ProjObject.ProjMode,{'inside','outside'})


 89  answer=msgbox_uvmat('INPUT_TXT','set bin size for histograms (or keep ''auto'' by default)?','auto');


 90  ParamOut.ActionInput.VarMesh=str2num(answer);


 91  end


 92  end


 93  % check the existence of the first and last file in the series


 94  first_j=[];


 95  if isfield(Param.IndexRange,'first_j'); first_j=Param.IndexRange.first_j; end


 96  last_j=[];


 97  if isfield(Param.IndexRange,'last_j'); last_j=Param.IndexRange.last_j; end


 98  PairString='';


 99  if isfield(Param.IndexRange,'PairString'); PairString=Param.IndexRange.PairString; end


 100  [i1,i2,j1,j2] = get_file_index(Param.IndexRange.first_i,first_j,PairString);


 101  FirstFileName=fullfile_uvmat(Param.InputTable{1,1},Param.InputTable{1,2},Param.InputTable{1,3},...


 102  Param.InputTable{1,5},Param.InputTable{1,4},i1,i2,j1,j2);


 103  if ~exist(FirstFileName,'file')


 104  msgbox_uvmat('WARNING',['the first input file ' FirstFileName ' does not exist'])


 105  else


 106  [i1,i2,j1,j2] = get_file_index(Param.IndexRange.last_i,last_j,PairString);


 107  LastFileName=fullfile_uvmat(Param.InputTable{1,1},Param.InputTable{1,2},Param.InputTable{1,3},...


 108  Param.InputTable{1,5},Param.InputTable{1,4},i1,i2,j1,j2);


 109  if ~exist(LastFileName,'file')


 110  msgbox_uvmat('WARNING',['the last input file ' LastFileName ' does not exist'])


 111  end


 112  end


 113  return


 114  end


 115 


 116  %%%%%%%%%%%% STANDARD PART %%%%%%%%%%%%


 117  ParamOut=[];%default output


 118  %% read input parameters from an xml file if input is a file name (batch mode)


 119  checkrun=1;


 120  if ischar(Param)


 121  Param=xml2struct(Param);% read Param as input file (batch case)


 122  checkrun=0;


 123  end


 124  hseries=findobj(allchild(0),'Tag','series');


 125  RUNHandle=findobj(hseries,'Tag','RUN');%handle of RUN button in GUI series


 126  WaitbarHandle=findobj(hseries,'Tag','Waitbar');%handle of waitbar in GUI series


 127 


 128  %% define the directory for result file (with path=RootPath{1})


 129  OutputDir=[Param.OutputSubDir Param.OutputDirExt];


 130 


 131  %% root input file(s) name, type and index series


 132  RootPath=Param.InputTable(:,1);


 133  RootFile=Param.InputTable(:,3);


 134  SubDir=Param.InputTable(:,2);


 135  NomType=Param.InputTable(:,4);


 136  FileExt=Param.InputTable(:,5);


 137  hdisp=disp_uvmat('WAITING...','checking the file series',checkrun);


 138  [filecell,i1_series,i2_series,j1_series,j2_series]=get_file_series(Param);


 139  if ~isempty(hdisp),delete(hdisp),end;


 140  %%%%%%%%%%%%


 141  % The cell array filecell is the list of input file names, while


 142  % filecell{iview,fileindex}:


 143  % iview: line in the table corresponding to a given file series


 144  % fileindex: file index within the file series,


 145  % 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


 146  % i1_series(iview,fileindex) expresses the same indices as a 1D array in file indices


 147  %%%%%%%%%%%%


 148  NbView=numel(i1_series);%number of input file series (lines in InputTable)


 149  NbField_j=size(i1_series{1},1); %nb of fields for the j index (bursts or volume slices)


 150  NbField_i=size(i1_series{1},2); %nb of fields for the i index


 151  NbField=NbField_j*NbField_i; %total number of fields


 152 


 153  %% determine the file type on each line from the first input file


 154  ImageTypeOptions={'image','multimage','mmreader','video'};


 155  NcTypeOptions={'netcdf','civx','civdata'};


 156  for iview=1:NbView


 157  if ~exist(filecell{iview,1}','file')


 158  disp_uvmat('ERROR',['the first input file ' filecell{iview,1} ' does not exist'],checkrun)


 159  return


 160  end


 161  [FileInfo{iview},MovieObject{iview}]=get_file_info(filecell{iview,1});


 162  FileType{iview}=FileInfo{iview}.FileType;


 163  CheckImage{iview}=~isempty(find(strcmp(FileType{iview},ImageTypeOptions)));% =1 for images


 164  CheckNc{iview}=~isempty(find(strcmp(FileType{iview},NcTypeOptions)));% =1 for netcdf files


 165  if ~isempty(j1_series{iview})


 166  frame_index{iview}=j1_series{iview};


 167  else


 168  frame_index{iview}=i1_series{iview};


 169  end


 170  end


 171 


 172  %% calibration data and timing: read the ImaDoc files


 173  [XmlData,NbSlice_calib,time,errormsg]=read_multimadoc(RootPath,SubDir,RootFile,FileExt,i1_series,i2_series,j1_series,j2_series);


 174  if size(time,1)>1


 175  diff_time=max(max(diff(time)));


 176  if diff_time>0


 177  disp_uvmat('WARNING',['times of series differ by (max) ' num2str(diff_time)],checkrun)


 178  end


 179  end


 180 


 181  %% coordinate transform or other user defined transform


 182  transform_fct='';%default


 183  if isfield(Param,'FieldTransform')&&~isempty(Param.FieldTransform.TransformName)


 184  addpath(Param.FieldTransform.TransformPath)


 185  transform_fct=str2func(Param.FieldTransform.TransformName);


 186  rmpath(Param.FieldTransform.TransformPath)


 187  if isfield(Param,'TransformInput')


 188  XmlData{1}.TransformInput=Param.TransformInput;


 189  end


 190  end


 191 


 192  %%%%%%%%%%%% END STANDARD PART %%%%%%%%%%%%


 193  % EDIT FROM HERE


 194 


 195  %% check the validity of input file types and set the output file type


 196  if CheckImage{1}


 197  FileExtOut='.png'; % write result as .png images for image inputs


 198  elseif CheckNc{1}


 199  FileExtOut='.nc';% write result as .nc files for netcdf inputs


 200  else


 201  disp_uvmat('ERROR',['invalid file type input ' FileType{1}],checkrun)


 202  return


 203  end


 204  if NbView==2 && ~isequal(CheckImage{1},CheckImage{2})


 205  disp_uvmat('ERROR','input must be two image series or two netcdf file series',checkrun)


 206  return


 207  end


 208  if isfield(Param,'ProjObject') && ~strcmp(Param.ProjObject.Type,'plane')


 209  FileExtOut='.nc';% write result as .nc files (even for image input)


 210  end


 211 


 212  %% settings for the output file


 213  NomTypeOut=nomtype2pair(NomType{1});% determine the index nomenclature type for the output file


 214  first_i=i1_series{1}(1);


 215  last_i=i1_series{1}(end);


 216  if isempty(j1_series{1})% if there is no second index j


 217  first_j=1;last_j=1;


 218  else


 219  first_j=j1_series{1}(1);


 220  last_j=j1_series{1}(end);


 221  end


 222 


 223  %% Set field names and velocity types


 224  InputFields{1}=[];%default (case of images)


 225  if NbView==2


 226  InputFields{2}=[];%default (case of images)


 227  end


 228  if isfield(Param,'InputFields')


 229  InputFields{1}=Param.InputFields;


 230  if NbView==2


 231  InputFields{2}=Param.InputFields;%default


 232  if isfield(Param.InputFields,'FieldName_1')


 233  InputFields{2}.FieldName=Param.InputFields.FieldName_1;


 234  if isfield(Param.InputFields,'VelType_1')


 235  InputFields{2}.VelType=Param.InputFields.VelType_1;


 236  end


 237  end


 238  end


 239  end


 240  nbfiles=0;%counter of the successfully read files (bad files are skipped)


 241  VarMesh=[];


 242  if isfield(Param,'ProjObject') && ismember(Param.ProjObject.ProjMode,{'inside','outside'})


 243  if isfield(Param,'ActionInput') && isfield(Param.ActionInput,'VarMesh')%case of histograms


 244  VarMesh=Param.ActionInput.VarMesh;


 245  else


 246  VarMesh=[];


 247  disp_uvmat('WARNING','automatic bin size for histograms, select aver_stat again to set the value',checkrun)


 248  end


 249  end


 250  %%%%%%%%%%%%%%%% loop on field indices %%%%%%%%%%%%%%%%


 251  for index=1:NbField


 252  update_waitbar(WaitbarHandle,index/NbField)


 253  if ~isempty(RUNHandle)&& ~strcmp(get(RUNHandle,'BusyAction'),'queue')


 254  disp('program stopped by user')


 255  break


 256  end


 257 


 258  %%%%%%%%%%%%%%%% loop on views (input lines) %%%%%%%%%%%%%%%%


 259  for iview=1:NbView


 260  % reading input file(s)


 261  [Data{iview},tild,errormsg] = read_field(filecell{iview,index},FileType{iview},InputFields{iview},frame_index{iview}(index));


 262  if ~isempty(errormsg)


 263  errormsg=['error of input reading: ' errormsg];


 264  break


 265  end


 266  if ~isempty(NbSlice_calib)


 267  Data{iview}.ZIndex=mod(i1_series{iview}(index)1,NbSlice_calib{iview})+1;%Zindex for phys transform


 268  end


 269  end


 270  %%%%%%%%%%%%%%%% end loop on views (input lines) %%%%%%%%%%%%%%%%


 271  %%%%%%%%%%%% END STANDARD PART %%%%%%%%%%%%


 272  % EDIT FROM HERE


 273 


 274  if isempty(errormsg)


 275  Field=Data{1}; % default input field structure


 276  nbfiles=nbfiles+1; %increment the file counter


 277  %% coordinate transform (or other user defined transform)


 278  if ~isempty(transform_fct)


 279  switch nargin(transform_fct)


 280  case 4


 281  if length(Data)==2


 282  Field=transform_fct(Data{1},XmlData{1},Data{2},XmlData{2});


 283  else


 284  Field=transform_fct(Data{1},XmlData{1});


 285  end


 286  case 3


 287  if length(Data)==2


 288  Field=transform_fct(Data{1},XmlData{1},Data{2});


 289  else


 290  Field=transform_fct(Data{1},XmlData{1});


 291  end


 292  case 2


 293  Field=transform_fct(Data{1},XmlData{1});


 294  case 1


 295  Field=transform_fct(Data{1});


 296  end


 297  end


 298 


 299  %% field projection on an object


 300  if Param.CheckObject


 301  if strcmp(Param.ProjObject.ProjMode,'interp_tps')


 302  Field=tps_coeff_field(Field,check_proj_tps);% calculate tps coefficients if needed


 303  end


 304  [Field,errormsg]=proj_field(Field,Param.ProjObject,VarMesh);


 305  if ~isempty(errormsg)


 306  disp_uvmat('ERROR',['error in aver_stat/proj_field:' errormsg],checkrun)


 307  return


 308  end


 309  end


 310 


 311  %%%%%%%%%%%% MAIN RUNNING OPERATIONS %%%%%%%%%%%%


 312  if nbfiles==1 %first field


 313  time_1=[];


 314  if isfield(Field,'Time')


 315  time_1=Field.Time(1);


 316  end


 317  DataOut=Field;%outcome reproduces the first (projected) field by default


 318  DataOut.Conventions='uvmat'; %suppress Conventions='uvmat/civdata' for civ input files


 319  if isfield(Param,'ProjObject')&& ismember(Param.ProjObject.ProjMode,{'inside','outside'})%case of histograms


 320  for ivar=1:numel(Field.ListVarName)% list of variable names before projection (histogram)


 321  VarName=Field.ListVarName{ivar};


 322  if isfield(Data{1},VarName)


 323  DataOut.(VarName)=Field.(VarName);


 324  DataOut.([VarName 'Histo'])=zeros(size(DataOut.(VarName)));


 325  VarMesh=DataOut.(VarName)(2)DataOut.(VarName)(1);


 326  end


 327  end


 328  disp(['mesh for histogram = ' num2str(VarMesh)])


 329  else


 330  errorvar=zeros(numel(Field.ListVarName));%index of errorflag associated to each variable


 331  if isfield(Field,'VarAttribute')


 332  for ivar=1:numel(Field.ListVarName)


 333  VarName=Field.ListVarName{ivar};


 334  DataOut.(VarName)=zeros(size(DataOut.(VarName)));% initiate each field to zero


 335  NbData.(VarName)=zeros(size(DataOut.(VarName)));% initiate the nbre of good data to zero


 336 


 337  for iivar=1:length(Field.VarAttribute)


 338  if isequal(Field.VarDimName{iivar},Field.VarDimName{ivar})&& isfield(Field.VarAttribute{iivar},'Role')...


 339  && strcmp(Field.VarAttribute{iivar}.Role,'errorflag')


 340  errorvar(ivar)=iivar; % index of the errorflag variable corresponding to ivar


 341  end


 342  end


 343  end


 344  DataOut.ListVarName(errorvar(errorvar~=0))=[]; %remove errorflag from result


 345  DataOut.VarDimName(errorvar(errorvar~=0))=[]; %remove errorflag from result


 346  DataOut.VarAttribute(errorvar(errorvar~=0))=[]; %remove errorflag from result


 347  else


 348  for ivar=1:numel(Field.ListVarName)


 349  VarName=Field.ListVarName{ivar};


 350  DataOut.(VarName)=zeros(size(DataOut.(VarName)));% initiate each field to zero


 351  NbData.(VarName)=zeros(size(DataOut.(VarName)));% initiate the nbre of good data to zero


 352  end


 353  end


 354 


 355  end


 356  end %current field


 357  for ivar=1:length(DataOut.ListVarName)


 358  VarName=DataOut.ListVarName{ivar};


 359  sizmean=size(DataOut.(VarName));


 360  siz=size(Field.(VarName));


 361  if isfield(Param,'ProjObject') && ismember(Param.ProjObject.ProjMode,{'inside','outside'})


 362  if isfield(Data{1},VarName)


 363  MaxValue=max(DataOut.(VarName));% current max of histogram absissa


 364  MinValue=min(DataOut.(VarName));% current min of histogram absissa


 365  % VarMesh=Field.VarAttribute{ivar}.Mesh;


 366  MaxIndex=round(MaxValue/VarMesh);


 367  MinIndex=round(MinValue/VarMesh);


 368  MaxIndex_new=round(max(Field.(VarName)/VarMesh));% max of the current field


 369  MinIndex_new=round(min(Field.(VarName)/VarMesh));


 370  if MaxIndex_new>MaxIndex% the variable max for the current field exceeds the previous one


 371  DataOut.(VarName)=[DataOut.(VarName) VarMesh*(MaxIndex+1:MaxIndex_new)];% append the new variable values


 372  DataOut.([VarName 'Histo'])=[DataOut.([VarName 'Histo']) zeros(1,MaxIndex_newMaxIndex)]; % append the new histo values


 373  end


 374  if MinIndex_new <= MinIndex1


 375  DataOut.(VarName)=[VarMesh*(MinIndex_new:MinIndex1) DataOut.(VarName)];% insert the new variable values


 376  DataOut.([VarName 'Histo'])=[zeros(1,MinIndexMinIndex_new) DataOut.([VarName 'Histo'])];% insert the new histo values


 377  ind_start=1;


 378  else


 379  ind_start=MinIndex_newMinIndex+1;


 380  end


 381  DataOut.([VarName 'Histo'])(ind_start:ind_start+MaxIndex_newMinIndex_new)=...


 382  DataOut.([VarName 'Histo'])(ind_start:ind_start+MaxIndex_newMinIndex_new)+Field.([VarName 'Histo']);


 383  end


 384  elseif ~isequal(DataOut.(VarName),0)&& ~isequal(siz,sizmean)


 385  disp_uvmat('ERROR',['unequal size of input field ' VarName ', need to project on a grid'],checkrun)


 386  return


 387  else


 388  if errorvar(ivar)==0


 389  check_bad=isnan(Field.(VarName));%=0 for NaN data values, 1 else


 390  else


 391  check_bad=isnan(Field.(VarName))  Field.(Field.ListVarName{errorvar(ivar)})~=0;%=0 for NaN or error flagged data values, 1 else


 392  end


 393  Field.(VarName)(check_bad)=0; %set to zero NaN or data marked by error flag


 394  DataOut.(VarName)=DataOut.(VarName)+ double(Field.(VarName)); % update the sum


 395  NbData.(VarName)=NbData.(VarName)+ ~check_bad;% records the number of data for each point


 396  end


 397  end


 398  %%%%%%%%%%%% END MAIN RUNNING OPERATIONS %%%%%%%%%%%%


 399  else


 400  disp(errormsg)


 401  end


 402  end


 403  %%%%%%%%%%%%%%%% end loop on field indices %%%%%%%%%%%%%%%%


 404  if ~(isfield(Param,'ProjObject') && ismember(Param.ProjObject.ProjMode,{'inside','outside'}))


 405  for ivar=1:length(Field.ListVarName)


 406  VarName=Field.ListVarName{ivar};


 407  DataOut.(VarName)=DataOut.(VarName)./NbData.(VarName); % normalize the mean


 408  end


 409  end


 410  nbmissing=NbFieldnbfiles;


 411  if nbmissing~=0


 412  disp_uvmat('WARNING',[num2str(nbmissing) ' input files are missing or skipped'],checkrun)


 413  end


 414  if isempty(time) % time is read from files


 415  if isfield(Field,'Time')


 416  time_end=Field.Time(1);%last time read


 417  if ~isempty(time_1)


 418  DataOut.Time=time_1;


 419  DataOut.Time_end=time_end;


 420  end


 421  end


 422  else % time from ImaDoc prevails if it exists


 423  DataOut.Time=time(1);


 424  DataOut.Time_end=time(end);


 425  end


 426 


 427  %% writing the result file


 428  OutputFile=fullfile_uvmat(RootPath{1},OutputDir,RootFile{1},FileExtOut,NomTypeOut,first_i,last_i,first_j,last_j);


 429  if strcmp(FileExtOut,'.png') %case of images


 430  if isequal(FileInfo{1}.BitDepth,16)(numel(FileInfo)==2 &&isequal(FileInfo{2}.BitDepth,16))


 431  DataOut.A=uint16(DataOut.A);


 432  imwrite(DataOut.A,OutputFile,'BitDepth',16); % case of 16 bit images


 433  else


 434  DataOut.A=uint8(DataOut.A);


 435  imwrite(DataOut.A,OutputFile,'BitDepth',8); % case of 16 bit images


 436  end


 437  disp([OutputFile ' written']);


 438  else %case of netcdf file , determine global attributes


 439  errormsg=struct2nc(OutputFile,DataOut); %save result file


 440  if isempty(errormsg)


 441  disp([OutputFile ' written']);


 442  else


 443  disp(['error in writting result file: ' errormsg])


 444  end


 445  end % end averaging loop


 446 


 447  %% open the result file with uvmat (in RUN mode)


 448  if checkrun


 449  uvmat(OutputFile)% open the last result file with uvmat


 450  end

