Ignore:
Timestamp:
Feb 16, 2015, 12:15:23 AM (6 years ago)
Author:
sommeria
Message:

histo improved

File:
1 edited

Legend:

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

    r869 r871  
    7171    ParamOut.OutputDirExt='.stat';%set the output dir extension
    7272    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 the existence of the first file in the series
     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_Y-N','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 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=str2double(answer);
     91        end
     92    end
     93    % check the existence of the first and last file in the series
    7494    first_j=[];
    7595    if isfield(Param.IndexRange,'first_j'); first_j=Param.IndexRange.first_j; end
     
    87107        LastFileName=fullfile_uvmat(Param.InputTable{1,1},Param.InputTable{1,2},Param.InputTable{1,3},...
    88108        Param.InputTable{1,5},Param.InputTable{1,4},i1,i2,j1,j2);
    89         if ~exist(FirstFileName,'file')
     109        if ~exist(LastFileName,'file')
    90110             msgbox_uvmat('WARNING',['the last input file ' LastFileName ' does not exist'])
    91111        end
     
    170190 % EDIT FROM HERE
    171191
    172 %% check the validity of  input file types
     192%% check the validity of  input file types and set the output file type
    173193if CheckImage{1}
    174194    FileExtOut='.png'; % write result as .png images for image inputs
     
    182202    disp_uvmat('ERROR','input must be two image series or two netcdf file series',checkrun)
    183203    return
     204end
     205if isfield(Param,'ProjObject') && ~strcmp(Param.ProjObject.Type,'plane')
     206      FileExtOut='.nc';% write result as .nc files (even for image input)
    184207end
    185208
     
    212235    end
    213236end
    214 % for i_slice=1:NbSlice
    215 % index_slice=i_slice:NbSlice:nbfield;% select file indices of the slice
    216 nbfiles=0;
    217 % nbmissing=0;
    218 
     237nbfiles=0;%counter of the successfully read files (bad files are skipped)
     238VarMesh=NaN;
     239if isfield(Param,'ProjObject') && ismember(Param.ProjObject.ProjMode,{'inside','outside'})&& isfield(Param.ActionInput,'VarMesh')%case of histograms
     240     VarMesh=Param.ActionInput.VarMesh;
     241end
    219242%%%%%%%%%%%%%%%% loop on field indices %%%%%%%%%%%%%%%%
    220243for index=1:NbField
     
    243266    if isempty(errormsg)
    244267        Field=Data{1}; % default input field structure
     268        nbfiles=nbfiles+1; %increment the file counter
    245269        %% coordinate transform (or other user defined transform)
    246270        if ~isempty(transform_fct)
     
    265289        end
    266290       
    267         %% calculate tps coefficients if needed
    268         if isfield(Param,'ProjObject')&&isfield(Param.ProjObject,'ProjMode')&& strcmp(Param.ProjObject.ProjMode,'interp_tps')
    269             Field=tps_coeff_field(Field,check_proj_tps);
    270         end
    271        
    272         %field projection on an object
     291        %% field projection on an object
    273292        if Param.CheckObject
    274             [Field,errormsg]=proj_field(Field,Param.ProjObject);
     293            if strcmp(Param.ProjObject.ProjMode,'interp_tps')
     294                Field=tps_coeff_field(Field,check_proj_tps);% calculate tps coefficients if needed
     295            end
     296            [Field,errormsg]=proj_field(Field,Param.ProjObject,VarMesh);
    275297            if ~isempty(errormsg)
    276298                disp_uvmat('ERROR',['error in aver_stat/proj_field:' errormsg],checkrun)
     
    278300            end
    279301        end
    280         nbfiles=nbfiles+1;
    281        
     302             
    282303        %%%%%%%%%%%% MAIN RUNNING OPERATIONS  %%%%%%%%%%%%
    283304        if nbfiles==1 %first field
     
    286307                time_1=Field.Time(1);
    287308            end
    288             DataOut=Field;%default
    289             DataOut.Conventions='uvmat'; %suppress Conventions='uvmat/civdata' for civ input files
    290             errorvar=zeros(numel(Field.ListVarName));%index of errorflag associated to each variable
    291             for ivar=1:numel(Field.ListVarName)
    292                 VarName=Field.ListVarName{ivar};
    293                 DataOut.(VarName)=zeros(size(DataOut.(VarName)));% initiate each field to zero
    294                 NbData.(VarName)=zeros(size(DataOut.(VarName)));% initiate the nbre of good data to zero
    295                 for iivar=1:length(Field.VarAttribute)
    296                     if isequal(Field.VarDimName{iivar},Field.VarDimName{ivar})&& isfield(Field.VarAttribute{iivar},'Role')...
    297                             && strcmp(Field.VarAttribute{iivar}.Role,'errorflag')
    298                         errorvar(ivar)=iivar; % index of the errorflag variable corresponding to ivar
     309            DataOut=Field;%outcome reproduces the first (projected) field by default
     310            DataOut.Conventions='uvmat'; %suppress Conventions='uvmat/civdata' for civ input files         
     311            if ismember(Param.ProjObject.ProjMode,{'inside','outside'})%case of histograms
     312                for ivar=1:numel(Field.ListVarName)% list of variable names before projection (histogram)
     313                    VarName=Field.ListVarName{ivar};
     314                    if isfield(Data{1},VarName)
     315                        DataOut.(VarName)=Field.(VarName);
     316                        DataOut.([VarName 'Histo'])=zeros(size(DataOut.(VarName)));
     317                        VarMesh=DataOut.(VarName)(2)-DataOut.(VarName)(1);
    299318                    end
    300319                end
    301             end
    302             DataOut.ListVarName(errorvar(errorvar~=0))=[]; %remove errorflag from result
    303             DataOut.VarDimName(errorvar(errorvar~=0))=[]; %remove errorflag from result
    304             DataOut.VarAttribute(errorvar(errorvar~=0))=[]; %remove errorflag from result
     320                disp(['mesh for histogram = ' num2str(VarMesh)])
     321            else
     322                errorvar=zeros(numel(Field.ListVarName));%index of errorflag associated to each variable
     323                for ivar=1:numel(Field.ListVarName)
     324                    VarName=Field.ListVarName{ivar};
     325                    DataOut.(VarName)=zeros(size(DataOut.(VarName)));% initiate each field to zero
     326                    NbData.(VarName)=zeros(size(DataOut.(VarName)));% initiate the nbre of good data to zero
     327                    for iivar=1:length(Field.VarAttribute)
     328                        if isequal(Field.VarDimName{iivar},Field.VarDimName{ivar})&& isfield(Field.VarAttribute{iivar},'Role')...
     329                                && strcmp(Field.VarAttribute{iivar}.Role,'errorflag')
     330                            errorvar(ivar)=iivar; % index of the errorflag variable corresponding to ivar
     331                        end
     332                    end
     333                end
     334                DataOut.ListVarName(errorvar(errorvar~=0))=[]; %remove errorflag from result
     335                DataOut.VarDimName(errorvar(errorvar~=0))=[]; %remove errorflag from result
     336                DataOut.VarAttribute(errorvar(errorvar~=0))=[]; %remove errorflag from result
     337            end
    305338        end   %current field
    306339        for ivar=1:length(DataOut.ListVarName)
    307             VarName=Field.ListVarName{ivar};
     340            VarName=DataOut.ListVarName{ivar};
    308341            sizmean=size(DataOut.(VarName));
    309342            siz=size(Field.(VarName));
    310             if ~isequal(DataOut.(VarName),0)&& ~isequal(siz,sizmean)
     343            if ismember(Param.ProjObject.ProjMode,{'inside','outside'})
     344                if isfield(Data{1},VarName)
     345                    MaxValue=max(DataOut.(VarName));% current max of histogram absissa
     346                    MinValue=min(DataOut.(VarName));% current min of histogram absissa
     347%                     VarMesh=Field.VarAttribute{ivar}.Mesh;
     348                    MaxIndex=round(MaxValue/VarMesh);
     349                    MinIndex=round(MinValue/VarMesh);
     350                    MaxIndex_new=round(max(Field.(VarName)/VarMesh));% max of the current field
     351                    MinIndex_new=round(min(Field.(VarName)/VarMesh));
     352                    if MaxIndex_new>MaxIndex% the variable max for the current field exceeds the previous one
     353                        DataOut.(VarName)=[DataOut.(VarName) VarMesh*(MaxIndex+1:MaxIndex_new)];% append the new variable values
     354                        DataOut.([VarName 'Histo'])=[DataOut.([VarName 'Histo']) zeros(1,MaxIndex_new-MaxIndex)]; % append the new histo values                   
     355                    end
     356                    if MinIndex_new <= MinIndex-1
     357                        DataOut.(VarName)=[VarMesh*(MinIndex_new:MinIndex-1) DataOut.(VarName)];% insert the new variable values
     358                        DataOut.([VarName 'Histo'])=[zeros(1,MinIndex-MinIndex_new) DataOut.([VarName 'Histo'])];% insert the new histo values
     359                        ind_start=1;
     360                    else
     361                        ind_start=MinIndex_new-MinIndex+1;
     362                    end
     363                    DataOut.([VarName 'Histo'])(ind_start:ind_start+MaxIndex_new-MinIndex_new)=...
     364                        DataOut.([VarName 'Histo'])(ind_start:ind_start+MaxIndex_new-MinIndex_new)+Field.([VarName 'Histo']);   
     365                end
     366            elseif ~isequal(DataOut.(VarName),0)&& ~isequal(siz,sizmean)
    311367                disp_uvmat('ERROR',['unequal size of input field ' VarName ', need to project  on a grid'],checkrun)
    312368                return
     
    328384end
    329385%%%%%%%%%%%%%%%% end loop on field indices %%%%%%%%%%%%%%%%
    330 
    331 for ivar=1:length(Field.ListVarName)
    332     VarName=Field.ListVarName{ivar};
    333     DataOut.(VarName)=DataOut.(VarName)./NbData.(VarName); % normalize the mean
     386if ~ismember(Param.ProjObject.ProjMode,{'inside','outside'})
     387    for ivar=1:length(Field.ListVarName)
     388        VarName=Field.ListVarName{ivar};
     389        DataOut.(VarName)=DataOut.(VarName)./NbData.(VarName); % normalize the mean
     390    end
    334391end
    335392nbmissing=NbField-nbfiles;
    336393if nbmissing~=0
    337     disp_uvmat('WARNING',[num2str(nbmissing) ' input files are missing or skipted'],checkrun)
     394    disp_uvmat('WARNING',[num2str(nbmissing) ' input files are missing or skipped'],checkrun)
    338395end
    339396if isempty(time) % time is read from files
     
    352409%% writing the result file
    353410OutputFile=fullfile_uvmat(RootPath{1},OutputDir,RootFile{1},FileExtOut,NomTypeOut,first_i,last_i,first_j,last_j);
    354 if CheckImage{1} %case of images
     411if strcmp(FileExtOut,'.png') %case of images
    355412    if isequal(FileInfo{1}.BitDepth,16)||(numel(FileInfo)==2 &&isequal(FileInfo{2}.BitDepth,16))
    356413        DataOut.A=uint16(DataOut.A);
     
    361418    end
    362419    disp([OutputFile ' written']);
    363 else %case of netcdf input file , determine global attributes
     420else %case of netcdf file , determine global attributes
    364421    errormsg=struct2nc(OutputFile,DataOut); %save result file
    365422    if isempty(errormsg)
Note: See TracChangeset for help on using the changeset viewer.