Ignore:
Timestamp:
Jun 24, 2016, 8:49:08 PM (5 years ago)
Author:
sommeria
Message:

various

File:
1 edited

Legend:

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

    r924 r959  
    111111        end
    112112    end
     113    % determine volume scan mode
     114        prompt = {'volume scan mode (Yes/No)'};
     115    dlg_title = 'determine volume scan';
     116    num_lines= 1;
     117    def     = { 'No'};
     118     answer=msgbox_uvmat('INPUT_Y-N','volume scan mode (OK/No)?');
     119%     answer = inputdlg(prompt,dlg_title,num_lines,def);
     120    if isempty(answer)
     121        return
     122    end
     123    %check input consistency
     124    if strcmp(answer,'Yes')
     125            ParamOut.NbSlice=1;% set NbSlice to 1 ( for i index)
     126            ParamOut.ActionInput.CheckVolume=1;
     127    end
    113128    return
    114129end
     
    248263    end
    249264end
    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
     265
     266%% set volume scan
     267NbSlice_j=1;
     268index_series=1:size(filecell,2);
     269first_j_out=first_j;
     270last_j_out=last_j;
     271if isfield(Param,'ActionInput') && isfield(Param.ActionInput,'CheckVolume') ...
     272        && Param.ActionInput.CheckVolume
     273   index_j=Param.IndexRange.first_j:Param.IndexRange.incr_j:Param.IndexRange.last_j;
     274   NbSlice_j=numel(index_j);
     275   index_i=index_series(1:NbSlice_j:end);
     276end
     277
     278%%% loop on slices (volume scan)
     279for islice=index_j
     280    if NbSlice_j>1
     281    first_j_out=islice;
     282    last_j_out=islice;
    257283   
    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];
     284    end
     285    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     286    %%%%%%%%%%%%%%%% loop on field indices %%%%%%%%%%%%%%%%
     287    for index=index_i+index_j(islice)
     288        update_waitbar(WaitbarHandle,index/NbField)
     289        if ~isempty(RUNHandle)&& ~strcmp(get(RUNHandle,'BusyAction'),'queue')
     290            disp('program stopped by user')
    264291            break
    265292        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});
     293       
     294        %%%%%%%%%%%%%%%% loop on views (input lines) %%%%%%%%%%%%%%%%
     295        for iview=1:NbView
     296            % reading input file(s)
     297            %InputFile=fullfile_uvmat(RootPath{iview},SubDir{iview},RootFile{iview},FileExt{iview},NomType{iview},first_i,last_i,first_j_out,last_j_out);
     298            [Data{iview},tild,errormsg] = read_field(filecell{iview,index},FileType{iview},InputFields{iview},frame_index{iview}(index));
     299            if ~isempty(errormsg)
     300                errormsg=['error of input reading: ' errormsg];
     301                break
     302            end
     303            if ~isempty(NbSlice_calib)
     304                Data{iview}.ZIndex=mod(i1_series{iview}(index)-1,NbSlice_calib{iview})+1;%Zindex for phys transform
     305            end
     306        end
     307        %%%%%%%%%%%%%%%% end loop on views (input lines) %%%%%%%%%%%%%%%%
     308       
     309        if isempty(errormsg)
     310            Field=Data{1}; % default input field structure
     311            nbfiles=nbfiles+1; %increment the file counter
     312            %% coordinate transform (or other user defined transform)
     313            if ~isempty(transform_fct)
     314                switch nargin(transform_fct)
     315                    case 4
     316                        if length(Data)==2
     317                            Field=transform_fct(Data{1},XmlData{1},Data{2},XmlData{2});
     318                        else
     319                            Field=transform_fct(Data{1},XmlData{1});
     320                        end
     321                    case 3
     322                        if length(Data)==2
     323                            Field=transform_fct(Data{1},XmlData{1},Data{2});
     324                        else
     325                            Field=transform_fct(Data{1},XmlData{1});
     326                        end
     327                    case 2
     328                        Field=transform_fct(Data{1},XmlData{1});
     329                    case 1
     330                        Field=transform_fct(Data{1});
     331                end
     332            end
     333           
     334            %% field projection on an object
     335            if Param.CheckObject
     336                if strcmp(Param.ProjObject.ProjMode,'interp_tps')
     337                    Field=tps_coeff_field(Field,check_proj_tps);% calculate tps coefficients if needed
     338                end
     339                [Field,errormsg]=proj_field(Field,Param.ProjObject,VarMesh);
     340                if ~isempty(errormsg)
     341                    disp_uvmat('ERROR',['error in aver_stat/proj_field:' errormsg],checkrun)
     342                    return
     343                end
     344            end
     345           
     346            %%%%%%%%%%%% MAIN RUNNING OPERATIONS  %%%%%%%%%%%%
     347            if nbfiles==1 %first field
     348                time_1=[];
     349                if isfield(Field,'Time')
     350                    time_1=Field.Time(1);
     351                end
     352                DataOut=Field;%outcome reproduces the first (projected) field by default
     353                DataOut.Conventions='uvmat'; %suppress Conventions='uvmat/civdata' for civ input files
     354                if isfield(Param,'ProjObject')&& ismember(Param.ProjObject.ProjMode,{'inside','outside'})%case of histograms
     355                    for ivar=1:numel(Field.ListVarName)% list of variable names before projection (histogram)
     356                        VarName=Field.ListVarName{ivar};
     357                        if isfield(Data{1},VarName)
     358                            DataOut.(VarName)=Field.(VarName);
     359                            DataOut.([VarName 'Histo'])=zeros(size(DataOut.(VarName)));
     360                            VarMesh=DataOut.(VarName)(2)-DataOut.(VarName)(1);
     361                        end
     362                    end
     363                    disp(['mesh for histogram = ' num2str(VarMesh)])
     364                else
     365                    errorvar=zeros(numel(Field.ListVarName));%index of errorflag associated to each variable
     366                    if isfield(Field,'VarAttribute')
     367                        for ivar=1:numel(Field.ListVarName)
     368                            VarName=Field.ListVarName{ivar};
     369                            DataOut.(VarName)=zeros(size(DataOut.(VarName)));% initiate each field to zero
     370                            NbData.(VarName)=zeros(size(DataOut.(VarName)));% initiate the nbre of good data to zero
     371                           
     372                            for iivar=1:length(Field.VarAttribute)
     373                                if isequal(Field.VarDimName{iivar},Field.VarDimName{ivar})&& isfield(Field.VarAttribute{iivar},'Role')...
     374                                        && strcmp(Field.VarAttribute{iivar}.Role,'errorflag')
     375                                    errorvar(ivar)=iivar; % index of the errorflag variable corresponding to ivar
     376                                end
     377                            end
     378                        end
     379                        DataOut.ListVarName(errorvar(errorvar~=0))=[]; %remove errorflag from result
     380                        DataOut.VarDimName(errorvar(errorvar~=0))=[]; %remove errorflag from result
     381                        DataOut.VarAttribute(errorvar(errorvar~=0))=[]; %remove errorflag from result
    283382                    else
    284                         Field=transform_fct(Data{1},XmlData{1});
     383                        for ivar=1:numel(Field.ListVarName)
     384                            VarName=Field.ListVarName{ivar};
     385                            DataOut.(VarName)=zeros(size(DataOut.(VarName)));% initiate each field to zero
     386                            NbData.(VarName)=zeros(size(DataOut.(VarName)));% initiate the nbre of good data to zero
     387                        end
    285388                    end
    286                 case 3
    287                     if length(Data)==2
    288                         Field=transform_fct(Data{1},XmlData{1},Data{2});
     389                   
     390                end
     391            end   %current field
     392            for ivar=1:length(DataOut.ListVarName)
     393                VarName=DataOut.ListVarName{ivar};
     394                sizmean=size(DataOut.(VarName));
     395                siz=size(Field.(VarName));
     396                if isfield(Param,'ProjObject') && ismember(Param.ProjObject.ProjMode,{'inside','outside'})
     397                    if isfield(Data{1},VarName)
     398                        MaxValue=max(DataOut.(VarName));% current max of histogram absissa
     399                        MinValue=min(DataOut.(VarName));% current min of histogram absissa
     400                        %                     VarMesh=Field.VarAttribute{ivar}.Mesh;
     401                        MaxIndex=round(MaxValue/VarMesh);
     402                        MinIndex=round(MinValue/VarMesh);
     403                        MaxIndex_new=round(max(Field.(VarName)/VarMesh));% max of the current field
     404                        MinIndex_new=round(min(Field.(VarName)/VarMesh));
     405                        if MaxIndex_new>MaxIndex% the variable max for the current field exceeds the previous one
     406                            DataOut.(VarName)=[DataOut.(VarName) VarMesh*(MaxIndex+1:MaxIndex_new)];% append the new variable values
     407                            DataOut.([VarName 'Histo'])=[DataOut.([VarName 'Histo']) zeros(1,MaxIndex_new-MaxIndex)]; % append the new histo values
     408                        end
     409                        if MinIndex_new <= MinIndex-1
     410                            DataOut.(VarName)=[VarMesh*(MinIndex_new:MinIndex-1) DataOut.(VarName)];% insert the new variable values
     411                            DataOut.([VarName 'Histo'])=[zeros(1,MinIndex-MinIndex_new) DataOut.([VarName 'Histo'])];% insert the new histo values
     412                            ind_start=1;
     413                        else
     414                            ind_start=MinIndex_new-MinIndex+1;
     415                        end
     416                        DataOut.([VarName 'Histo'])(ind_start:ind_start+MaxIndex_new-MinIndex_new)=...
     417                            DataOut.([VarName 'Histo'])(ind_start:ind_start+MaxIndex_new-MinIndex_new)+Field.([VarName 'Histo']);
     418                    end
     419                elseif ~isequal(DataOut.(VarName),0)&& ~isequal(siz,sizmean)
     420                    disp_uvmat('ERROR',['unequal size of input field ' VarName ', need to project  on a grid'],checkrun)
     421                    return
     422                else
     423                    if errorvar(ivar)==0
     424                        check_bad=isnan(Field.(VarName));%=0 for NaN data values, 1 else
    289425                    else
    290                         Field=transform_fct(Data{1},XmlData{1});
     426                        check_bad=isnan(Field.(VarName)) | Field.(Field.ListVarName{errorvar(ivar)})~=0;%=0 for NaN or error flagged data values, 1 else
    291427                    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_new-MaxIndex)]; % append the new histo values                   
    373                     end
    374                     if MinIndex_new <= MinIndex-1
    375                         DataOut.(VarName)=[VarMesh*(MinIndex_new:MinIndex-1) DataOut.(VarName)];% insert the new variable values
    376                         DataOut.([VarName 'Histo'])=[zeros(1,MinIndex-MinIndex_new) DataOut.([VarName 'Histo'])];% insert the new histo values
    377                         ind_start=1;
    378                     else
    379                         ind_start=MinIndex_new-MinIndex+1;
    380                     end
    381                     DataOut.([VarName 'Histo'])(ind_start:ind_start+MaxIndex_new-MinIndex_new)=...
    382                         DataOut.([VarName 'Histo'])(ind_start:ind_start+MaxIndex_new-MinIndex_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 %%%%%%%%%%%%%%%%
     428                    Field.(VarName)(check_bad)=0; %set to zero NaN or data marked by error flag
     429                    DataOut.(VarName)=DataOut.(VarName)+ double(Field.(VarName)); % update the sum
     430                    NbData.(VarName)=NbData.(VarName)+ ~check_bad;% records the number of data for each point
     431                end
     432            end
     433            %%%%%%%%%%%%   END MAIN RUNNING OPERATIONS  %%%%%%%%%%%%
     434        else
     435            disp(errormsg)
     436        end
     437    end
     438    %%%%%%%%%%%%%%%% end loop on field indices %%%%%%%%%%%%%%%%
     439
    404440if ~(isfield(Param,'ProjObject') && ismember(Param.ProjObject.ProjMode,{'inside','outside'}))
    405441    for ivar=1:length(Field.ListVarName)
     
    426462
    427463%% writing the result file
    428 OutputFile=fullfile_uvmat(RootPath{1},OutputDir,RootFile{1},FileExtOut,NomTypeOut,first_i,last_i,first_j,last_j);
     464OutputFile=fullfile_uvmat(RootPath{1},OutputDir,RootFile{1},FileExtOut,NomTypeOut,first_i,last_i,first_j_out,last_j_out);
    429465if strcmp(FileExtOut,'.png') %case of images
    430466    if isequal(FileInfo{1}.BitDepth,16)||(numel(FileInfo)==2 &&isequal(FileInfo{2}.BitDepth,16))
     
    444480    end
    445481end  % end averaging  loop
     482end
    446483
    447484%% open the result file with uvmat (in RUN mode)
Note: See TracChangeset for help on using the changeset viewer.