Changeset 255 for trunk/src


Ignore:
Timestamp:
May 23, 2011, 9:27:19 AM (13 years ago)
Author:
sommeria
Message:

acver_stat corrected for multi slice case

File:
1 edited

Legend:

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

    r218 r255  
    4747%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    4848
    49 %projection object
     49%% projection object
    5050test_object=get(hseries.GetObject,'Value');
    5151if test_object%isfield(Series,'sethandles')
     
    5959end
    6060
    61 %root input file and type
     61%% root input file and type
    6262if ~iscell(Series.RootPath)% case of a single input field series
    6363    num_i1={num_i1};num_j1={num_j1};num_i2={num_i2};num_j2={num_j2};
     
    9292end
    9393
    94 %Number of input series: this function  accepts two input file series at most (then it operates on the difference of fields)
     94%% Number of input series: this function  accepts two input file series at most (then it operates on the difference of fields)
    9595nbview=length(RootPath);
    9696if nbview>2 
     
    107107end
    108108
    109 %determine image type
     109%% determine image type
    110110hhh=which('mmreader');
    111111for iview=1:nbview
     
    133133end
    134134
    135 % number of slices
     135%% number of slices
    136136NbSlice=str2num(get(hseries.NbSlice,'String'));
    137137if isempty(NbSlice)
     
    140140NbSlice_name=num2str(NbSlice);
    141141
    142 % Field and velocity type (the same for the two views)
     142%% Field and velocity type (the same for the two views)
    143143Field_str=get(hseries.FieldMenu,'String');
    144144FieldName=[]; %default
     
    161161    end
    162162end
    163 %detect whether the two files are 'images' or 'netcdf'
    164 % testima=0;
    165 % testvol=0;
     163
     164%% get the velocity type
    166165testcivx=0;
    167 % testnc=0;
    168166FileExt=get(hseries.FileExt,'String');
    169 % test_movie=0;
    170 % for iview=1:nbview
    171 %      ext=FileExt{iview};
    172 %      form=imformats(ext([2:end]));
    173 %      if isequal(lower(ext),'.vol')
    174 %          testvol=testvol+1;
    175 %      elseif ~isempty(form)||isequal(lower(ext),'.avi')% if the extension corresponds to an image format recognized by Matlab
    176 %          testima=testima+1;
    177 %      elseif isequal(ext,'.nc')
    178 %          testnc=testnc+1;
    179 %      end
    180 % end
    181 % if testvol
    182 %     msgbox_uvmat('ERROR','volume images not implemented yet')
    183 %     return
    184 % end
    185 % if testnc~=nbview && testima~=nbview && testvol~=nbview
    186 %     msgbox_uvmat('ERROR','compare two image series or two netcdf files with the same fields as input')
    187 %     return
    188 % end
    189167if ~isequal(FieldName,{'get_field...'})
    190168    testcivx=isequal(FileType{1},'netcdf');
    191169end
    192 % if ~isequal(FieldName,{'get_field...'})
    193 %     if isequal(FieldName,{''}) && ~testima
    194 %         msgbox_uvmat('ERROR','an input field needs to be selected')
    195 %         return
    196 %     end
    197 %     testcivx=testnc;
    198 % end
    199 
    200170if testcivx
    201171    VelType_str=get(hseries.VelTypeMenu,'String');
     
    209179end
    210180
    211 %Calibration data and timing: read the ImaDoc files
     181%% Calibration data and timing: read the ImaDoc files
    212182mode=''; %default
    213183timecell={};
     
    253223end
    254224
    255 %check coincidence in time
     225%% check coincidence in time
    256226multitime=0;
    257227if isempty(timecell)
     
    283253end
    284254
    285 %% Root name of output files (TO GENERALISE FOR TWO INPUT SERIES)
    286 subdir_result='aver_stat';
    287 pathdir=fullfile(RootPath{1},subdir_result);
    288 while exist(pathdir,'dir')
    289     subdir_result=[subdir_result '.0'];
    290     pathdir=fullfile(RootPath{1},subdir_result);
    291 end
    292 [m1,m2,m3]=mkdir(pathdir);
    293 if ~isequal(m2,'')
    294      msgbox_uvmat('CONFIRMATION',m2);%error message for directory creation
     255%% Name(s) of output file(s)
     256filebase_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)
     257%file extension of the result 
     258if testima %case of images
     259    ext_out='.png';
     260else
     261    ext_out='.nc';
     262end
     263subdir_result='aver_stat';%subdirectory for the results
     264pathdir=fullfile(RootPath{1},subdir_result);% full subdirectory name, including path
     265testexist=1;
     266while testexist
     267    pathdir=fullfile(RootPath{1},subdir_result);% full subdirectory name, including path
     268    if NbSlice==1% keep track of the first and lsat indices of the input files
     269        NomTypeOut=nomtype2pair(NomType{1},num_i2{end}(end)-num_i1{1}(1),num_j2{end}(end)-num_j1{1}(1));
     270        fileresult{1}=name_generator(filebase_out,num_i1{1}(1),num_j1{1}(1),ext_out,NomTypeOut,1,num_i2{end}(end),num_j2{end}(end),subdir_result);
     271        testexist=exist(fileresult{1},'file');
     272    else % simplified indexing with i_slice for multiple slices
     273        testexist=0;
     274        for i_slice=1:NbSlice
     275            fileresult{i_slice}=name_generator(filebase_out,i_slice,[],ext_out,'_1',1,i_slice,[],subdir_result);
     276            if exist(fileresult{i_slice},'file')
     277                testexist=1;
     278                break
     279            end
     280        end
     281    end
     282    if testexist
     283        subdir_result=[subdir_result '.0'];
     284    end
     285end
     286% create result directory if needed
     287if ~exist(pathdir,'dir')
     288    [m1,m2,m3]=mkdir(pathdir);
     289    if ~isequal(m2,'')
     290        msgbox_uvmat('CONFIRMATION',m2);%error message for directory creation
     291    end
    295292end
    296293[xx,msg2] = fileattrib(pathdir,'+w','g'); %yield writing access (+w) to user group (g)
    297294if ~strcmp(msg2,'')
    298     msgbox_uvmat('ERROR',['pb of permission for ' pathdir ': ' msg2])%error message for directory creation
     295    msgbox_uvmat('ERROR',['pb of permission for ' pathdir ': ' msg2])%error message for writting access
    299296    return
    300297end
    301 filebase_out=filebase{1};
    302 NomTypeOut=nomtype2pair(NomType{1},num_i2{end}(end)-num_i1{1}(1),num_j2{end}(end)-num_j1{1}(1));
    303 
    304 % coordinate transform or other user defined transform
     298
     299%% coordinate transform or other user defined transform
    305300transform_fct=[];%default
    306301if isfield(Series,'transform_fct')
     
    308303end
    309304
    310 %% slice loop
     305%% main loop
    311306siz=size(num_i1{1});
    312 lengthtot=siz(1)*siz(2);
    313 nbfield=floor(lengthtot/(siz(1)*NbSlice));%total number of i indexes (adjusted to an integer number of slices)
    314 nbfield_slice=nbfield*siz(1);% number of fields per slice
    315 
     307nbfield2=siz(1); %nb of consecutive fields at each level(burst
     308nbfield=siz(1)*siz(2);
     309nbfield=floor(nbfield/(nbfield2*NbSlice));%total number of i indexes (adjusted to an integer number of slices)
     310
     311% loop on slices
    316312for i_slice=1:NbSlice
    317    S=0; %initiate the image sum S
    318    nbfiles=0;
    319    nbmissing=0;
    320     %averaging loop
    321    for ifile=i_slice:NbSlice:lengthtot
     313    for ifield=1:nbfield
     314         indselect(:,ifield)=((ifield-1)*NbSlice+(i_slice-1))*nbfield2+[1:nbfield2]';%selected indices on the list of files of a slice
     315    end
     316    S=0; %initiate the image sum S
     317    nbfiles=0;
     318    nbmissing=0;
     319    % averaging loop
     320    for index=1:nbfield*nbfield2
    322321        stopstate=get(hseries.RUN,'BusyAction');
    323322        if isequal(stopstate,'queue') % enable STOP command
    324              update_waitbar(hseries.waitbar,WaitbarPos,ifile/lengthtot)
    325              for iview=1:nbview
     323            update_waitbar(hseries.waitbar,WaitbarPos,index/(nbfield*nbfield2))
     324            ifile=indselect(index);
     325            % reading input file(s)
     326            for iview=1:nbview
    326327                [filename]=...
    327                            name_generator(filebase{iview},num_i1{iview}(ifile),num_j1{iview}(ifile),FileExt{iview},NomType{iview},1,num_i2{iview}(ifile),num_j2{iview}(ifile),SubDir{iview});
    328                 if ~isequal(FileType{iview},'netcdf')               
     328                    name_generator(filebase{iview},num_i1{iview}(ifile),num_j1{iview}(ifile),FileExt{iview},NomType{iview},1,num_i2{iview}(ifile),num_j2{iview}(ifile),SubDir{iview});
     329                if ~isequal(FileType{iview},'netcdf')
    329330                    Data{iview}.ListVarName={'A'};
    330331                    Data{iview}.AName='image';
     
    341342                        case 'image'
    342343                            A=imread(filename);
    343                     end 
    344                     Data{iview}.ListVarName={'AY','AX','A'}; % 
     344                    end
     345                    Data{iview}.ListVarName={'AY','AX','A'}; %
    345346                    Atype{iview}=class(A);
    346347                    npy=size(A,1);
     
    348349                    nbcolor=size(A,3);
    349350                    if nbcolor==3
    350                          Data{iview}.VarDimName={'AY','AX',{'AY','AX','rgb'}};
     351                        Data{iview}.VarDimName={'AY','AX',{'AY','AX','rgb'}};
    351352                    else
    352                          Data{iview}.VarDimName={'AY','AX',{'AY','AX'}};
    353                     end 
     353                        Data{iview}.VarDimName={'AY','AX',{'AY','AX'}};
     354                    end
    354355                    Data{iview}.AY=[npy-0.5 0.5];
    355356                    Data{iview}.AX=[0.5 npx-0.5];
     
    359360                    [Data{iview},VelTypeOut]=read_civxdata(filename,FieldName,VelType);
    360361                else
    361                     [Data{iview},var_detect]=nc2struct(filename,SubField.ListVarName); %read the corresponding input data               
     362                    [Data{iview},var_detect]=nc2struct(filename,SubField.ListVarName); %read the corresponding input data
    362363                    Data{iview}.VarAttribute=SubField.VarAttribute;
    363                 end 
     364                end
    364365                if isfield(Data{iview},'Txt')
    365366                    msgbox_uvmat('ERROR',['error of input reading: ' Data{iview}.Txt])
    366367                    return
    367368                end
    368              end   
    369 
    370              % coordinate transform (or other user defined transform)
    371              if ~isempty(transform_fct)
    372                  % z index
     369            end
     370           
     371            % coordinate transform (or other user defined transform)
     372            if ~isempty(transform_fct)
    373373                if ~isempty(NbSlice_calib)
    374374                    Data{iview}.ZIndex=mod(num_i1{iview}(ifile)-1,NbSlice_calib{1})+1;%Zindex for phys transform
     
    382382                    Data{1}=transform_fct(Data{1},XmlData{1});
    383383                end
    384              end     
     384            end
     385           
     386            % field calculation (vort, div...)
    385387            if testcivx
    386                     Data{iview}=calc_field(FieldName,Data{iview});%calculate field (vort..)
    387             end
     388                Data{iview}=calc_field(FieldName,Data{iview});%calculate field (vort..)
     389            end
     390           
     391            % field substration (for two input file series)
    388392            if length(Data)==2
    389393                [Field,errormsg]=sub_field(Data{1},Data{2}); %substract the two fields
     
    397401            if test_object
    398402                [Field,errormsg]=proj_field(Field,ProjObject);
    399                  if ~isempty(errormsg)
     403                if ~isempty(errormsg)
    400404                    msgbox_uvmat('ERROR',['error in aver_stat/proj_field:' errormsg])
    401405                    return
    402406                end
    403              end                                                       
    404                 nbfiles=nbfiles+1;
    405                 if nbfiles==1 %first field
    406                     time_1=[];
    407                     if isfield(Field,'Time')
    408                         time_1=Field.Time(1);
     407            end
     408            nbfiles=nbfiles+1;
     409            if nbfiles==1 %first field
     410                time_1=[];
     411                if isfield(Field,'Time')
     412                    time_1=Field.Time(1);
     413                end
     414                DataMean=Field;%default
     415            else
     416                for ivar=1:length(Field.ListVarName)
     417                    VarName=Field.ListVarName{ivar};
     418                    eval(['sizmean=size(DataMean.' VarName ');']);
     419                    eval(['siz=size(Field.' VarName ');']);
     420                    if ~isequal(siz,sizmean)
     421                        msgbox_uvmat('ERROR',['unequal size of input field ' VarName ', need to project  on a grid'])
     422                        return
     423                    else
     424                        eval(['DataMean.' VarName '=DataMean.' VarName '+ Field.' VarName ';']); % update the sum
    409425                    end
    410                     DataMean=Field;%default
    411                 else
    412                     for ivar=1:length(Field.ListVarName)
    413                         VarName=Field.ListVarName{ivar};
    414                         eval(['sizmean=size(DataMean.' VarName ');']);
    415                         eval(['siz=size(Field.' VarName ');']);
    416                         if ~isequal(siz,sizmean)
    417                            msgbox_uvmat('ERROR',['unequal size of input field ' VarName ', need to interpolate on a grid'])
    418                            return
    419                         else
    420                             eval(['DataMean.' VarName '=DataMean.' VarName '+ Field.' VarName ';']); % update the sum
    421                         end
    422                     end
    423                 end
    424         end
    425     end %end averaging loop
     426                end
     427            end
     428        end
     429    end
     430    %end averaging loop
    426431    for ivar=1:length(Field.ListVarName)
    427432        VarName=Field.ListVarName{ivar};
     
    444449        DataMean.Time_end=time(end,num_i1{end}(end),num_j1{end}(end));
    445450    end
    446 
     451   
    447452    %writing the result file
    448    if testima   
    449         [filemean]=name_generator(filebase_out,num_i1{1}(1),num_j1{1}(1),'.png',NomTypeOut,1,num_i2{end}(end),num_j2{end}(end),subdir_result);
    450         if exist(filemean,'file')
    451             backupfile=filemean;
    452             testexist=2;
    453             while testexist==2
    454                 backupfile=[backupfile(1:end-4) '~.png'];
    455                 testexist=exist(backupfile,'file');
    456             end
    457             [success,message]=copyfile(filemean,backupfile);%make backup
    458             if ~isequal(success,1)
    459                 msgbox_uvmat('ERROR',['previous file result ' filemean ' already exists, problem in backup'])
    460                 return
    461             end
    462         end
     453    if testima %case of images
    463454        if isequal(Atype{1},'uint16')
    464             imwrite(uint16(DataMean.A),filemean,'BitDepth',16);
     455            imwrite(uint16(DataMean.A),fileresult{i_slice},'BitDepth',16); % case of 16 bit images
    465456        else
    466             imwrite(uint8(DataMean.A),filemean,'BitDepth',8);
    467         end
    468         display([filemean ' written']);
     457            imwrite(uint8(DataMean.A),fileresult{i_slice},'BitDepth',8); % case of 8 bit images
     458        end
     459        display([fileresult{i_slice} ' written']);
    469460    else %case of netcdf input file , determine global attributes
    470461        DataMean.ListGlobalAttribute=[DataMean.ListGlobalAttribute {Series.Action}];
     
    477468        if isfield(DataMean,'Time')
    478469            DataMean.ListGlobalAttribute=[DataMean.ListGlobalAttribute {'Time','Time_end'}];
    479         end 
    480         filemean=name_generator(filebase_out,num_i1{1}(1),num_j1{1}(1),'.nc',NomTypeOut,1,num_i2{end}(end),num_j2{end}(end),subdir_result);
    481         if exist(filemean,'file')
    482             backupfile=filemean;
    483             testexist=2;
    484             while testexist==2
    485                 backupfile=[backupfile(1:end-3) '~.nc'];
    486                 testexist=exist(backupfile,'file');
    487             end
    488             [success,message]=copyfile(filemean,backupfile);%make backup
    489             if ~isequal(success,1)
    490                 msgbox_uvmat('ERROR',['previous file result ' filemean ' already exists, problem in backup'])
    491                 display(['previous file result ' filemean ' already exists, problem in backup'])
    492                 return
    493             end
    494         end
    495         errormsg=struct2nc(filemean,DataMean); %save result file
     470        end
     471        errormsg=struct2nc(fileresult{i_slice},DataMean); %save result file
    496472        if isempty(errormsg)
    497             display([filemean ' written']);
     473            display([fileresult{i_slice} ' written']);
    498474        else
    499475            msgbox_uvmat('ERROR',['error in writting result file: ' errormsg])
    500476            display(errormsg)
    501477        end
    502    end
    503 end
    504 
     478    end  % end averaging  loop
     479end % end loop on slices
     480
     481%% reproduce ImaDoc/GeometryCalib for image series
     482if isfield(XmlData{1},'GeometryCalib') && ~isempty(XmlData{1}.GeometryCalib)
     483    [pp,RootFile]=fileparts(filebase_out);
     484    outputxml=fullfile(pathdir,[RootFile '.xml'])
     485    errormsg=update_imadoc(XmlData{1}.GeometryCalib,outputxml);% introduce the calibration data in the xml file
     486    if strcmp(errormsg,'')
     487        display(['GeometryCalib transferred to ' outputxml])
     488    else
     489        msgbox_uvmat('ERROR',errormsg);
     490    end
     491end
     492
     493%% open the result file with uvmat
    505494hget_field=findobj(allchild(0),'name','get_field');%find the get_field... GUI
    506495delete(hget_field)
    507 uvmat(filemean)
     496uvmat(fileresult{end})
Note: See TracChangeset for help on using the changeset viewer.