Changeset 959 for trunk/src/series/aver_stat.m
 Timestamp:
 Jun 24, 2016, 8:49:08 PM (5 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

trunk/src/series/aver_stat.m
r924 r959 111 111 end 112 112 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_YN','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 113 128 return 114 129 end … … 248 263 end 249 264 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 265 266 %% set volume scan 267 NbSlice_j=1; 268 index_series=1:size(filecell,2); 269 first_j_out=first_j; 270 last_j_out=last_j; 271 if 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); 276 end 277 278 %%% loop on slices (volume scan) 279 for islice=index_j 280 if NbSlice_j>1 281 first_j_out=islice; 282 last_j_out=islice; 257 283 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') 264 291 break 265 292 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 283 382 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 285 388 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_newMaxIndex)]; % append the new histo values 408 end 409 if MinIndex_new <= MinIndex1 410 DataOut.(VarName)=[VarMesh*(MinIndex_new:MinIndex1) DataOut.(VarName)];% insert the new variable values 411 DataOut.([VarName 'Histo'])=[zeros(1,MinIndexMinIndex_new) DataOut.([VarName 'Histo'])];% insert the new histo values 412 ind_start=1; 413 else 414 ind_start=MinIndex_newMinIndex+1; 415 end 416 DataOut.([VarName 'Histo'])(ind_start:ind_start+MaxIndex_newMinIndex_new)=... 417 DataOut.([VarName 'Histo'])(ind_start:ind_start+MaxIndex_newMinIndex_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 289 425 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 291 427 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 %%%%%%%%%%%%%%%% 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 404 440 if ~(isfield(Param,'ProjObject') && ismember(Param.ProjObject.ProjMode,{'inside','outside'})) 405 441 for ivar=1:length(Field.ListVarName) … … 426 462 427 463 %% writing the result file 428 OutputFile=fullfile_uvmat(RootPath{1},OutputDir,RootFile{1},FileExtOut,NomTypeOut,first_i,last_i,first_j ,last_j);464 OutputFile=fullfile_uvmat(RootPath{1},OutputDir,RootFile{1},FileExtOut,NomTypeOut,first_i,last_i,first_j_out,last_j_out); 429 465 if strcmp(FileExtOut,'.png') %case of images 430 466 if isequal(FileInfo{1}.BitDepth,16)(numel(FileInfo)==2 &&isequal(FileInfo{2}.BitDepth,16)) … … 444 480 end 445 481 end % end averaging loop 482 end 446 483 447 484 %% open the result file with uvmat (in RUN mode)
Note: See TracChangeset
for help on using the changeset viewer.