Changeset 939 for trunk/src/series/merge_proj_polar.m
 Timestamp:
 Apr 6, 2016, 10:35:30 PM (9 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

trunk/src/series/merge_proj_polar.m
r938 r939 61 61 %% set the input elements needed on the GUI series when the function is selected in the menu ActionName or InputTable refreshed 62 62 if isstruct(Param) && isequal(Param.Action.RUN,0) 63 ParamOut.AllowInputSort='o ff';% allow alphabetic sorting of the list of input file SubDir (options 'off'/'on', 'off' by default)63 ParamOut.AllowInputSort='on';% allow alphabetic sorting of the list of input file SubDir (options 'off'/'on', 'off' by default) 64 64 ParamOut.WholeIndexRange='off';% prescribes the file index ranges from min to max (options 'off'/'on', 'off' by default) 65 65 ParamOut.NbSlice='off'; %nbre of slices ('off' by default) … … 68 68 ParamOut.FieldTransform = 'on';%can use a transform function 69 69 ParamOut.TransformPath=fullfile(fileparts(which('uvmat')),'transform_field');% path to transform functions (needed for compilation only) 70 ParamOut.ProjObject='o n';%can use projection object(option 'off'/'on',71 ParamOut.Mask='o n';%can use mask option (option 'off'/'on', 'off' by default)72 ParamOut.OutputDirExt='. mproj';%set the output dir extension70 ParamOut.ProjObject='off';%can use projection object(option 'off'/'on', 71 ParamOut.Mask='off';%can use mask option (option 'off'/'on', 'off' by default) 72 ParamOut.OutputDirExt='.polar';%set the output dir extension 73 73 ParamOut.OutputFileMode='NbInput';% '=NbInput': 1 output file per input file index, '=NbInput_i': 1 file per input file index i, '=NbSlice': 1 file per slice 74 74 %check the input files … … 83 83 if ~exist(FirstFileName,'file') 84 84 msgbox_uvmat('WARNING',['the first input file ' FirstFileName ' does not exist']) 85 elseif isequal(size(Param.InputTable,1),1) && ~isfield(Param,'ProjObject')86 msgbox_uvmat('WARNING','You may need a projection object of type plane for merge_proj')87 85 end 88 86 return 89 87 end 88 89 %%%% specific input parameters 90 % calculate the positions on which to interpolate 91 radius_ref=450;% radius of the mountain top 92 radius_shifted=120:2:120;% radius shifted by the radius of the origin at the topography summit 93 radius=radius_ref+radius_shifted;%radius from centre of the tank 94 azimuth_arclength=(120:2:500);%azimuth in arc length at origin position 95 azimuth=pi/2azimuth_arclength/radius_ref;%azimuth in radian 96 [Radius,Azimuth]=meshgrid(radius,azimuth); 97 XI=Radius.*cos(Azimuth);% set of x axis of the points where interpolqtion needs to be done 98 YI=Radius.*sin(Azimuth)radius_ref;% set of y axis of the points where interpolqtion needs to be done 99 FieldNames={'vec(U,V)';'curl(U,V)';'div(U,V)'}; 100 HeadData.ListVarName= {'radius','azimuth'} ; 101 HeadData.VarDimName={'radius','azimuth'}; 102 HeadData.VarAttribute={'coord_y','coord_x'} ; 103 HeadData.radius=radius_shifted; 104 HeadData.azimuth=azimuth_arclength; 90 105 91 106 %%%%%%%%%%%% STANDARD PART (DO NOT EDIT) %%%%%%%%%%%% … … 105 120 106 121 %% define the directory for result file (with path=RootPath{1}) 107 OutputDir=[ 'three_cameras'Param.OutputDirExt];% subdirectory for output files122 OutputDir=[Param.OutputSubDir Param.OutputDirExt];% subdirectory for output files 108 123 109 124 if ~isfield(Param,'InputFields') … … 115 130 RootFile=Param.InputTable(:,3); 116 131 SubDir=Param.InputTable(:,2); 117 NomType=Param.InputTable(:,4);132 %NomType=Param.InputTable(:,4); 118 133 FileExt=Param.InputTable(:,5); 119 134 [filecell,i1_series,i2_series,j1_series,j2_series]=get_file_series(Param); … … 176 191 time=mean(time,1); %averaged time taken for the merged field 177 192 193 %% height z 194 % position of projection plane 195 196 ProjObjectCoord=XmlData{1}.GeometryCalib.SliceCoord; 197 CoordUnit=XmlData{1}.GeometryCalib.CoordUnit; 198 for iview =2:numel(XmlData) 199 if ~(isfield(XmlData{iview},'GeometryCalib')&& isequal(XmlData{iview}.GeometryCalib.SliceCoord,ProjObjectCoord))... 200 disp('error: geometric calibration missing') 201 return 202 end 203 end 204 205 178 206 %% coordinate transform or other user defined transform 179 207 transform_fct='';%default fct handle … … 194 222 %% check the validity of input file types 195 223 for iview=1:NbView 196 if ~isequal(Check Image{iview},1)&&~isequal(CheckNc{iview},1)197 disp_uvmat('ERROR','input set of input series: need either netcdf either image series',checkrun)224 if ~isequal(CheckNc{iview},1) 225 disp_uvmat('ERROR','input files needs to be in netcdf (extension .nc)',checkrun) 198 226 return 199 227 end … … 201 229 202 230 %% output file type 203 if min(cell2mat(CheckImage))==1 && (~Param.CheckObject  strcmp(Param.ProjObject.Type,'plane')) 204 FileExtOut='.png'; %image output (input and proj result = image) 205 for iview=1:NbView 206 BitDepth(iview)=FileInfo{iview}.BitDepth; 207 end 208 BitDepth=max(BitDepth); 209 else 210 FileExtOut='.nc'; %netcdf output 211 end 231 FileExtOut='.nc'; %netcdf output 212 232 if isempty(j1_series{1}) 213 233 NomTypeOut='_1'; … … 215 235 NomTypeOut='_1_1'; 216 236 end 217 %NomTypeOut=NomType;% output file index will indicate the first and last ref index in the series218 237 RootFileOut=RootFile{1}; 219 238 for iview=2:NbView … … 224 243 end 225 244 226 %% mask (TODO: case of multilevels)227 MaskData=cell(NbView,1);228 if Param.CheckMask229 if ischar(Param.MaskTable)% case of a single mask (char chain)230 Param.MaskTable={Param.MaskTable};231 end232 for iview=1:numel(Param.MaskTable)233 if exist(Param.MaskTable{iview},'file')234 [MaskData{iview},tild,errormsg] = read_field(Param.MaskTable{iview},'image');235 if ~isempty(transform_fct) && nargin(transform_fct)>=2236 MaskData{iview}=transform_fct(MaskData{iview},XmlData{iview});237 end238 end239 end240 end241 242 %% Set field names and velocity types243 %use Param.InputFields for all views244 245 245 246 %% MAIN LOOP ON FIELDS … … 250 251 % nbmissing=0; 251 252 252 253 %%%%%%%%%%%%%%%% loop on field indices %%%%%%%%%%%%%%%% 253 254 tstart=tic; %used to record the computing time 254 255 CheckOverwrite=1;%default … … 257 258 end 258 259 for index=1:NbField 259 260 update_waitbar(WaitbarHandle,index/NbField) 260 261 if ~isempty(RUNHandle) && ~strcmp(get(RUNHandle,'BusyAction'),'queue') 261 262 disp('program stopped by user') … … 263 264 end 264 265 265 266 %% generating the name of the merged field 266 267 i1=i1_series{1}(index); 267 268 if ~isempty(i2_series{end}) … … 282 283 OutputFile=fullfile_uvmat(RootPath{1},OutputDir,RootFileOut,FileExtOut,NomTypeOut,i1,i2,j1,j2); 283 284 if ~CheckOverwrite && exist(OutputFile,'file') 284 disp(['existing output file ' OutputFile ' already exists, skip to next field']) 285 continue% skip iteration if the mode overwrite is desactivated and the result file already exists 286 end 287 285 disp(['existing output file ' OutputFile ' already exists, skip to next field']) 286 continue% skip iteration if the mode overwrite is desactivated and the result file already exists 287 end 288 289 %% z position 290 ZIndex=mod(i1_series{1}(index)1,NbSlice_calib{1})+1;%Zindex for phys transform 291 ZPos=ProjObjectCoord(ZIndex,3); 292 % radius of the topography section at z pos 293 TopoRadius=0; 294 ind_mask=[]; 295 if ZPos<20 296 TopoRadius=40*sin(acos((20+ZPos)/40)); 297 ind_mask=(XI'.*XI'+YI'.*YI')<TopoRadius*TopoRadius;% indidces of data to mask 298 end 299 288 300 %%%%%%%%%%%%%%%% loop on views (input lines) %%%%%%%%%%%%%%%% 289 301 Data=cell(1,NbView);%initiate the set Data 290 302 timeread=zeros(1,NbView); 291 303 for iview=1:NbView 292 %% reading input file(s) 304 %% reading input file(s) 293 305 [Data{iview},tild,errormsg] = read_field(filecell{iview,index},FileType{iview},ParamIn{iview},frame_index{iview}(index)); 294 306 if ~isempty(errormsg) … … 318 330 end 319 331 320 %% calculate tps coefficients if needed 321 check_proj_tps= 1; 322 Data{iview}=tps_coeff_field(Data{iview},check_proj_tps); 332 %% calculate tps coefficients 333 Data{iview}=tps_coeff_field(Data{iview},1); 323 334 324 %% projection on object (gridded plane) 325 if Param.CheckObject 326 [Data{iview},errormsg]=proj_field(Data{iview},Param.ProjObject); 327 if ~isempty(errormsg) 328 disp_uvmat('ERROR',['ERROR in merge_proge/proj_field: ' errormsg],checkrun) 329 return 330 end 331 end 332 335 %% projection on the polar grid 336 [DataOut,VarAttribute,errormsg]=calc_field_tps(Data{iview}.Coord_tps,Data{iview}.NbCentre,Data{iview}.SubRange,... 337 cat(3,Data{iview}.U_tps,Data{iview}.V_tps),FieldNames,cat(3,XI,YI)); 338 ListVarName=(fieldnames(DataOut))'; 339 ProjData{iview}=HeadData; 340 ProjData{iview}.ListVarName= [ProjData{iview}.ListVarName ListVarName]; 341 ProjData{iview}.VarDimName={'radius','azimuth'}; 342 ProjData{iview}.VarAttribute=[{'coord_x'} {'coord_y'} VarAttribute]; 343 for ivar=1:numel(ListVarName) 344 ProjData{iview}.VarDimName{ivar+2}={'radius','azimuth'}; 345 ProjData{iview}.(ListVarName{ivar})=(DataOut.(ListVarName{ivar}))'; 346 end 333 347 %% mask 334 if Param.CheckMask && ~isempty(MaskData{iview})335 [Data{iview},errormsg]=mask_proj(Data{iview},MaskData{iview});336 end348 % if Param.CheckMask && ~isempty(MaskData{iview}) 349 % [Data{iview},errormsg]=mask_proj(Data{iview},MaskData{iview}); 350 % end 337 351 end 338 352 %%%%%%%%%%%%%%%% END LOOP ON VIEWS %%%%%%%%%%%%%%%% 339 353 340 354 %% merge the NbView fields 341 [MergeData,errormsg]=merge_field( Data);355 [MergeData,errormsg]=merge_field(ProjData); 342 356 if ~isempty(errormsg) 343 357 disp_uvmat('ERROR',errormsg,checkrun); 344 358 return 345 359 end 346 360 347 361 %% time of the merged field: take the average of the different views 348 362 if ~isempty(time) 349 timeread=time(index); 363 timeread=time(index); 350 364 elseif ~isempty(find(timeread))% time defined from ImaDoc 351 365 timeread=mean(timeread(timeread~=0));% take average over times form the files (when defined) 352 366 else 353 timeread=index;% take time=file index 354 end 355 356 367 timeread=index;% take time=file index 368 end 369 370 %% rotating the velocity vectors to the local axis of the polatr coordinates 371 Unew=MergeData.U.*sin(Azimuth')MergeData.V.*cos(Azimuth'); 372 MergeData.V=MergeData.U.*cos(Azimuth')+MergeData.V.*sin(Azimuth'); 373 MergeData.U=Unew; 374 if ~isempty(ind_mask) 375 MergeData.U(ind_mask)=NaN; 376 MergeData.V(ind_mask)=NaN; 377 MergeData.curl(ind_mask)=NaN; 378 MergeData.div(ind_mask)=NaN; 379 end 380 357 381 %% recording the merged field 358 if strcmp(FileExtOut,'.png') %output as image 359 if BitDepth==8 360 imwrite(uint8(MergeData.A),OutputFile,'BitDepth',8) 361 else 362 imwrite(uint16(MergeData.A),OutputFile,'BitDepth',16) 363 end 364 if index==1 365 %write xml calibration file, using the first file 366 siz=size(MergeData.A); 367 npy=siz(1); 368 npx=siz(2); 369 if isfield(MergeData,'Coord_x') && isfield(MergeData,'Coord_y') 370 Rangx=MergeData.Coord_x; 371 Rangy=MergeData.Coord_y; 372 elseif isfield(MergeData,'AX')&& isfield(MergeData,'AY') 373 Rangx=[MergeData.AX(1) MergeData.AX(end)]; 374 Rangy=[MergeData.AY(1) MergeData.AY(end)]; 375 else 376 Rangx=[0.5 npx0.5]; 377 Rangy=[npy0.5 0.5];%default 382 383 MergeData.ListGlobalAttribute={'Conventions','Project','InputFile_1','InputFile_end','NbCoord','NbDim','CoordUnit','ZPos'}; 384 MergeData.Conventions='uvmat'; 385 MergeData.NbCoord=2; 386 MergeData.NbDim=2; 387 MergeData.CoordUnit=CoordUnit; 388 MergeData.ZPos=ZPos; 389 % time interval of PIV 390 Dt=[]; 391 if isfield(Data{1},'Dt')&& isnumeric(Data{1}.Dt) 392 Dt=Data{1}.Dt; 393 end 394 for iview =2:numel(Data) 395 if ~(isfield(Data{iview},'Dt')&& isequal(Data{iview}.Dt,Dt)) 396 Dt=[];%dt not the same for all fields 397 end 398 end 399 if ~isempty(timeread) 400 MergeData.ListGlobalAttribute=[MergeData.ListGlobalAttribute {'Time'}]; 401 MergeData.Time=timeread; 402 end 403 404 % time unit 405 if isfield(Data{1},'TimeUnit') 406 TimeUnit=Data{1}.TimeUnit; 407 for iview =2:numel(Data) 408 if ~(isfield(Data{iview},'TimeUnit')&& isequal(Data{iview}.TimeUnit,TimeUnit)) 409 TimeUnit=[];%TimeUnit not the same for all fields 378 410 end 379 pxcmx=(npx1)/(Rangx(2)Rangx(1)); 380 pxcmy=(npy1)/(Rangy(1)Rangy(2)); 381 T_x=pxcmx*Rangx(1)+0.5; 382 T_y=pxcmy*Rangy(2)+0.5; 383 GeometryCal.CalibrationType='rescale'; 384 GeometryCal.CoordUnit=MergeData.CoordUnit; 385 GeometryCal.focal=1; 386 GeometryCal.R=[pxcmx,0,0;0,pxcmy,0;0,0,1]; 387 GeometryCal.Tx_Ty_Tz=[T_x T_y 1]; 388 ImaDoc.GeometryCalib=GeometryCal; 389 t=struct2xml(ImaDoc); 390 t=set(t,1,'name','ImaDoc'); 391 save(t,[fileparts(OutputFile) '.xml']) 392 end 393 394 else %output as netcdf files 395 MergeData.ListGlobalAttribute={'Conventions','Project','InputFile_1','InputFile_end','NbCoord','NbDim'}; 396 MergeData.Conventions='uvmat'; 397 MergeData.NbCoord=2; 398 MergeData.NbDim=2; 399 % time interval of PIV 400 Dt=[]; 401 if isfield(Data{1},'Dt')&& isnumeric(Data{1}.Dt) 402 Dt=Data{1}.Dt; 403 end 404 for iview =2:numel(Data) 405 if ~(isfield(Data{iview},'Dt')&& isequal(Data{iview}.Dt,Dt)) 406 Dt=[];%dt not the same for all fields 407 end 408 end 409 if ~isempty(timeread) 410 MergeData.ListGlobalAttribute=[MergeData.ListGlobalAttribute {'Time'}]; 411 MergeData.Time=timeread; 412 end 413 % position of projection plane 414 if isfield(Data{1},'ProjObjectCoord')&& isfield(Data{1},'ProjObjectAngle') 415 'test' 416 ProjObjectCoord=Data{1}.ProjObjectCoord; 417 ProjObjectAngle=Data{1}.ProjObjectAngle; 418 for iview =2:numel(Data) 419 if ~(isfield(Data{iview},'ProjObjectCoord')&& isequal(Data{iview}.ProjObjectCoord,ProjObjectCoord))... 420 ~(isfield(Data{iview},'ProjObjectAngle')&& isequal(Data{iview}.ProjObjectAngle,ProjObjectAngle)) 421 ProjObjectCoord=[];%dt not the same for all fields 422 end 423 end 424 if ~isempty(ProjObjectCoord) 425 MergeData.ListGlobalAttribute=[MergeData.ListGlobalAttribute {'ProjObjectCoord'} {'ProjObjectAngle'}]; 426 MergeData.ProjObjectCoord=ProjObjectCoord; 427 MergeData.ProjObjectAngle=ProjObjectAngle; 428 end 429 end 430 % coord unit 431 if isfield(Data{1},'CoordUnit') 432 CoordUnit=Data{1}.CoordUnit; 433 for iview =2:numel(Data) 434 if ~(isfield(Data{iview},'CoordUnit')&& isequal(Data{iview}.CoordUnit,CoordUnit)) 435 CoordUnit=[];%CoordUnit not the same for all fields 436 end 437 end 438 if ~isempty(CoordUnit) 439 MergeData.ListGlobalAttribute=[MergeData.ListGlobalAttribute {'CoordUnit'}]; 440 MergeData.CoordUnit=CoordUnit; 441 end 442 end 443 % time unit 444 if isfield(Data{1},'TimeUnit') 445 TimeUnit=Data{1}.TimeUnit; 446 for iview =2:numel(Data) 447 if ~(isfield(Data{iview},'TimeUnit')&& isequal(Data{iview}.TimeUnit,TimeUnit)) 448 TimeUnit=[];%TimeUnit not the same for all fields 449 end 450 end 451 if ~isempty(TimeUnit) 452 MergeData.ListGlobalAttribute=[MergeData.ListGlobalAttribute {'TimeUnit'}]; 453 MergeData.TimeUnit=TimeUnit; 454 end 455 end 456 error=struct2nc(OutputFile,MergeData);%save result file 457 if isempty(error) 458 disp(['output file ' OutputFile ' written']) 459 else 460 disp(error) 461 end 462 end 463 end 411 end 412 if ~isempty(TimeUnit) 413 MergeData.ListGlobalAttribute=[MergeData.ListGlobalAttribute {'TimeUnit'}]; 414 MergeData.TimeUnit=TimeUnit; 415 end 416 end 417 error=struct2nc(OutputFile,MergeData);%save result file 418 if isempty(error) 419 disp(['output file ' OutputFile ' written']) 420 else 421 disp(error) 422 end 423 end 424 464 425 ellapsed_time=toc(tstart); 465 426 disp(['total ellapsed time ' num2str(ellapsed_time/60,2) ' minutes'])
Note: See TracChangeset
for help on using the changeset viewer.