Ignore:
Timestamp:
Apr 2, 2013, 9:13:42 AM (8 years ago)
Author:
sommeria
Message:

various bugs repaired . civ_series further developed

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/src/series/aver_stat.m

    r596 r598  
    166166end
    167167
    168 %% MAIN LOOP ON SLICES
    169168% for i_slice=1:NbSlice
    170    % index_slice=i_slice:NbSlice:nbfield;% select file indices of the slice
    171     nbfiles=0;
    172     nbmissing=0;
    173 
    174     %%%%%%%%%%%%%%%% loop on field indices %%%%%%%%%%%%%%%%
    175     for index=1:nbfield
    176     %for index=index_slice
    177           if checkrun
    178                 stopstate=get(Param.RUNHandle,'BusyAction');
    179                 update_waitbar(Param.WaitbarHandle,index/nbfield)
    180           else
    181                 stopstate='queue';
    182           end
    183         if isequal(stopstate,'queue')% enable STOP command
    184            
    185         %%%%%%%%%%%%%%%% loop on views (input lines) %%%%%%%%%%%%%%%%
    186         for iview=1:nbview
    187             % reading input file(s)
    188             [Data{iview},tild,errormsg] = read_field(filecell{iview,index},FileType{iview},InputFields{iview},frame_index{iview}(index));
     169% index_slice=i_slice:NbSlice:nbfield;% select file indices of the slice
     170nbfiles=0;
     171nbmissing=0;
     172
     173%%%%%%%%%%%%%%%% loop on field indices %%%%%%%%%%%%%%%%
     174for index=1:nbfield
     175%for index=index_slice
     176      if checkrun
     177            stopstate=get(Param.RUNHandle,'BusyAction');
     178            update_waitbar(Param.WaitbarHandle,index/nbfield)
     179      else
     180            stopstate='queue';
     181      end
     182    if isequal(stopstate,'queue')% enable STOP command
     183
     184    %%%%%%%%%%%%%%%% loop on views (input lines) %%%%%%%%%%%%%%%%
     185    for iview=1:nbview
     186        % reading input file(s)
     187        [Data{iview},tild,errormsg] = read_field(filecell{iview,index},FileType{iview},InputFields{iview},frame_index{iview}(index));
     188        if ~isempty(errormsg)
     189            errormsg=['error of input reading: ' errormsg];
     190            break
     191        end
     192        if ~isempty(NbSlice_calib)
     193            Data{iview}.ZIndex=mod(i1_series{iview}(index)-1,NbSlice_calib{iview})+1;%Zindex for phys transform
     194        end
     195    end
     196    else
     197        errormsg='stop';
     198    end
     199    %%%%%%%%%%%%%%%% end loop on views (input lines) %%%%%%%%%%%%%%%%
     200    %%%%%%%%%%%% END STANDARD PART  %%%%%%%%%%%%
     201    % EDIT FROM HERE
     202
     203    if isempty(errormsg)
     204        Field=Data{1}; % default input field structure
     205        %% coordinate transform (or other user defined transform)
     206        if ~isempty(transform_fct)
     207            switch nargin(transform_fct)
     208                case 4
     209                    if length(Data)==2
     210                        Field=transform_fct(Data{1},XmlData{1},Data{2},XmlData{2});
     211                    else
     212                        Field=transform_fct(Data{1},XmlData{1});
     213                    end
     214                case 3
     215                    if length(Data)==2
     216                        Field=transform_fct(Data{1},XmlData{1},Data{2});
     217                    else
     218                        Field=transform_fct(Data{1},XmlData{1});
     219                    end
     220                case 2
     221                    Field=transform_fct(Data{1},XmlData{1});
     222                case 1
     223                    Field=transform_fct(Data{1});
     224            end
     225        end
     226
     227        %% calculate tps coefficients if needed
     228        if isfield(Param,'ProjObject')&&isfield(Param.ProjObject,'ProjMode')&& strcmp(Param.ProjObject.ProjMode,'interp_tps')
     229            Field=tps_coeff_field(Field,check_proj_tps);
     230        end
     231
     232        %field projection on an object
     233        if Param.CheckObject
     234            [Field,errormsg]=proj_field(Field,Param.ProjObject);
    189235            if ~isempty(errormsg)
    190                 errormsg=['error of input reading: ' errormsg];
    191                 break
    192             end
    193             if ~isempty(NbSlice_calib)
    194                 Data{iview}.ZIndex=mod(i1_series{iview}(index)-1,NbSlice_calib{iview})+1;%Zindex for phys transform
    195             end
    196         end
    197         else
    198             errormsg='stop';
    199         end
    200         %%%%%%%%%%%%%%%% end loop on views (input lines) %%%%%%%%%%%%%%%%
    201         %%%%%%%%%%%% END STANDARD PART  %%%%%%%%%%%%
    202         % EDIT FROM HERE
    203    
    204         if isempty(errormsg)
    205             Field=Data{1}; % default input field structure
    206             %% coordinate transform (or other user defined transform)
    207             if ~isempty(transform_fct)
    208                 switch nargin(transform_fct)
    209                     case 4
    210                         if length(Data)==2
    211                             Field=transform_fct(Data{1},XmlData{1},Data{2},XmlData{2});
    212                         else
    213                             Field=transform_fct(Data{1},XmlData{1});
    214                         end
    215                     case 3
    216                         if length(Data)==2
    217                             Field=transform_fct(Data{1},XmlData{1},Data{2});
    218                         else
    219                             Field=transform_fct(Data{1},XmlData{1});
    220                         end
    221                     case 2
    222                         Field=transform_fct(Data{1},XmlData{1});
    223                     case 1
    224                         Field=transform_fct(Data{1});
     236                msgbox_uvmat('ERROR',['error in aver_stat/proj_field:' errormsg])
     237                return
     238            end
     239        end
     240        nbfiles=nbfiles+1;
     241
     242        %%%%%%%%%%%% MAIN RUNNING OPERATIONS  %%%%%%%%%%%%
     243        %update sum
     244        if nbfiles==1 %first field
     245            time_1=[];
     246            if isfield(Field,'Time')
     247                time_1=Field.Time(1);
     248            end
     249            DataOut=Field;%default
     250            for ivar=1:length(Field.ListVarName)
     251                VarName=Field.ListVarName{ivar};
     252                DataOut.(VarName)=double(DataOut.(VarName));
     253            end
     254        else   %current field
     255            for ivar=1:length(Field.ListVarName)
     256                VarName=Field.ListVarName{ivar};
     257                sizmean=size(DataOut.(VarName));
     258                siz=size(Field.(VarName));
     259                if ~isequal(DataOut.(VarName),0)&& ~isequal(siz,sizmean)
     260                    msgbox_uvmat('ERROR',['unequal size of input field ' VarName ', need to project  on a grid'])
     261                    return
     262                else
     263                    DataOut.(VarName)=DataOut.(VarName)+ double(Field.(VarName)); % update the sum
    225264                end
    226265            end
    227            
    228             %% calculate tps coefficients if needed
    229             if isfield(Param,'ProjObject')&&isfield(Param.ProjObject,'ProjMode')&& strcmp(Param.ProjObject.ProjMode,'interp_tps')
    230                 Field=tps_coeff_field(Field,check_proj_tps);
    231             end
    232 
    233             %field projection on an object
    234             if Param.CheckObject
    235                 [Field,errormsg]=proj_field(Field,Param.ProjObject);
    236                 if ~isempty(errormsg)
    237                     msgbox_uvmat('ERROR',['error in aver_stat/proj_field:' errormsg])
    238                     return
    239                 end
    240             end
    241             nbfiles=nbfiles+1;
    242            
    243             %%%%%%%%%%%% MAIN RUNNING OPERATIONS  %%%%%%%%%%%%
    244             %update sum
    245             if nbfiles==1 %first field
    246                 time_1=[];
    247                 if isfield(Field,'Time')
    248                     time_1=Field.Time(1);
    249                 end
    250                 DataOut=Field;%default
    251                 for ivar=1:length(Field.ListVarName)
    252                     VarName=Field.ListVarName{ivar};
    253                     DataOut.(VarName)=double(DataOut.(VarName));
    254                 end
    255             else   %current field
    256                 for ivar=1:length(Field.ListVarName)
    257                     VarName=Field.ListVarName{ivar};
    258                     sizmean=size(DataOut.(VarName));
    259                     siz=size(Field.(VarName));
    260                     if ~isequal(DataOut.(VarName),0)&& ~isequal(siz,sizmean)
    261                         msgbox_uvmat('ERROR',['unequal size of input field ' VarName ', need to project  on a grid'])
    262                         return
    263                     else
    264                         DataOut.(VarName)=DataOut.(VarName)+ double(Field.(VarName)); % update the sum
    265                     end
    266                 end
    267             end
    268             %%%%%%%%%%%%   END MAIN RUNNING OPERATIONS  %%%%%%%%%%%%
    269         else
    270             display(errormsg)
    271         end
    272     end
    273     %%%%%%%%%%%%%%%% end loop on field indices %%%%%%%%%%%%%%%%
    274    
    275     for ivar=1:length(Field.ListVarName)
    276         VarName=Field.ListVarName{ivar};
    277         DataOut.(VarName)=DataOut.(VarName)/nbfiles; % normalize the mean
    278     end
    279     if nbmissing~=0
    280         msgbox_uvmat('WARNING',[num2str(nbmissing) ' input files are missing or skipted'])
    281     end
    282     if isempty(time) % time is read from files
    283         if isfield(Field,'Time')
    284             time_end=Field.Time(1);%last time read
    285             if ~isempty(time_1)
    286                 DataOut.Time=time_1;
    287                 DataOut.Time_end=time_end;
    288             end
    289         end
    290     else  % time from ImaDoc prevails if it exists
     266        end
     267        %%%%%%%%%%%%   END MAIN RUNNING OPERATIONS  %%%%%%%%%%%%
     268    else
     269        display(errormsg)
     270    end
     271end
     272%%%%%%%%%%%%%%%% end loop on field indices %%%%%%%%%%%%%%%%
     273
     274for ivar=1:length(Field.ListVarName)
     275    VarName=Field.ListVarName{ivar};
     276    DataOut.(VarName)=DataOut.(VarName)/nbfiles; % normalize the mean
     277end
     278if nbmissing~=0
     279    msgbox_uvmat('WARNING',[num2str(nbmissing) ' input files are missing or skipted'])
     280end
     281if isempty(time) % time is read from files
     282    if isfield(Field,'Time')
     283        time_end=Field.Time(1);%last time read
     284        if ~isempty(time_1)
     285            DataOut.Time=time_1;
     286            DataOut.Time_end=time_end;
     287        end
     288    end
     289else  % time from ImaDoc prevails if it exists
    291290%         j1=1;%default
    292291%         if ~isempty(j1_series{1})
    293292%             j1=j1_series{1};
    294293%         end
    295         %DataOut.Time=time(1,i1_series{1}(1),j1);
    296         %DataOut.Time_end=time(end,i1_series{end}(end),j1_series{end}(end));
    297         DataOut.Time=time(1);
    298         DataOut.Time_end=time(end);
    299     end
    300    
    301     %writting the result file
    302     OutputFile=fullfile_uvmat(RootPath{1},OutputDir,RootFile{1},FileExtOut,NomTypeOut,i1_series{1}(1),i1_series{1}(end),j1_series{1}(1),j1_series{1}(end));
    303     if CheckImage{1} %case of images
    304         if isequal(FileInfo{1}.BitDepth,16)||(numel(FileInfo)==2 &&isequal(FileInfo{2}.BitDepth,16))
    305             DataOut.A=uint16(DataOut.A);
    306             imwrite(DataOut.A,OutputFile,'BitDepth',16); % case of 16 bit images
    307         else
    308             DataOut.A=uint8(DataOut.A);
    309             imwrite(DataOut.A,OutputFile,'BitDepth',8); % case of 16 bit images
    310         end
     294    %DataOut.Time=time(1,i1_series{1}(1),j1);
     295    %DataOut.Time_end=time(end,i1_series{end}(end),j1_series{end}(end));
     296    DataOut.Time=time(1);
     297    DataOut.Time_end=time(end);
     298end
     299
     300%writting the result file
     301OutputFile=fullfile_uvmat(RootPath{1},OutputDir,RootFile{1},FileExtOut,NomTypeOut,i1_series{1}(1),i1_series{1}(end),j1_series{1}(1),j1_series{1}(end));
     302if CheckImage{1} %case of images
     303    if isequal(FileInfo{1}.BitDepth,16)||(numel(FileInfo)==2 &&isequal(FileInfo{2}.BitDepth,16))
     304        DataOut.A=uint16(DataOut.A);
     305        imwrite(DataOut.A,OutputFile,'BitDepth',16); % case of 16 bit images
     306    else
     307        DataOut.A=uint8(DataOut.A);
     308        imwrite(DataOut.A,OutputFile,'BitDepth',8); % case of 16 bit images
     309    end
     310    display([OutputFile ' written']);
     311else %case of netcdf input file , determine global attributes
     312    errormsg=struct2nc(OutputFile,DataOut); %save result file
     313    if isempty(errormsg)
    311314        display([OutputFile ' written']);
    312     else %case of netcdf input file , determine global attributes
    313         errormsg=struct2nc(OutputFile,DataOut); %save result file
    314         if isempty(errormsg)
    315             display([OutputFile ' written']);
    316         else
    317             msgbox_uvmat('ERROR',['error in writting result file: ' errormsg])
    318             display(errormsg)
    319         end
    320     end  % end averaging  loop
     315    else
     316        msgbox_uvmat('ERROR',['error in writting result file: ' errormsg])
     317        display(errormsg)
     318    end
     319end  % end averaging  loop
    321320% end
    322321%%%%%%%%%%%%%%%% end loop on slices %%%%%%%%%%%%%%%%
Note: See TracChangeset for help on using the changeset viewer.