Changeset 985 for trunk/src/series/civ2vel_3C.m
 Timestamp:
 Jan 20, 2017, 5:17:25 PM (4 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

trunk/src/series/civ2vel_3C.m
r984 r985 1 %'civ2vel_3C': combine velocity fields from two cameras to get the three velocity components 2 % used with the GUI 'series': 3 % first input line =raw PIV camera 1 (image coordinates) 4 % second input line=raw PIV camera 2 (image coordinates) 5 % Three modes: 6 % 1) no additional input: measurements assumed in the reference plane (laser sheet) 7 % 2) measurement surface obtained by stereoscopic comparison of the images of the two cameras. 8 % third input line =surface z(x,y) given by correlation between camera 1 and 2 (expressed in phys apparent coordinates) 9 % 3) surface z(x,y) given by adding the displacements of each camera with a third intermediate camera (#3) used to reduce the 10 % to reduce the angle for stereoscopic view. 11 % third input line =correlation between camera 1 and 3 (expressed in phys apparent coordinates) 12 % fourth input line =corelation between camera 2 and 3 (expressed in phys apparent coordinates) 13 % 1 %'civ2vel_3C': combine velocity fields from two cameras to get three velocity components 14 2 % 15 3 % function ParamOut=civ2vel_3C(Param) … … 19 7 % 20 8 %INPUT: 21 %22 9 % In run mode, the input parameters are given as a Matlab structure Param copied from the GUI series. 23 10 % In batch mode, Param is the name of the corresponding xml file containing the same information … … 28 15 % see the current structure Param) 29 16 % .InputTable: cell of input file names, (several lines for multiple input) 30 % each line decomposed as {RootPath,SubDir,Rootfile,NomType,Extension}. 31 % 3 or 4 lines used as described above 17 % each line decomposed as {RootPath,SubDir,Rootfile,NomType,Extension} 32 18 % .OutputSubDir: name of the subdirectory for data outputs 33 19 % .OutputDirExt: directory extension for data outputs … … 36 22 % .ActionExt: fct extension ('.m', Matlab fct, '.sh', compiled Matlab fct 37 23 % .RUN =0 for GUI input, =1 for function activation 38 % .RunMode='local','background', 'cluster': type of function use 24 % .RunMode='local','background', 'cluster': type of function use 25 % 39 26 % .IndexRange: set the file or frame indices on which the action must be performed 40 27 % .InputFields: sub structure describing the input fields withfields … … 47 34 48 35 %======================================================================= 49 % Copyright 2008201 7, LEGI UMR 5519 / CNRS UGAGINP, Grenoble, France36 % Copyright 20082015, LEGI UMR 5519 / CNRS UJF GINP, Grenoble, France 50 37 % http://www.legi.grenobleinp.fr 51 38 % Joel.Sommeria  Joel.Sommeria (A) legi.cnrs.fr … … 65 52 66 53 function ParamOut=civ2vel_3C(Param) 67 54 disp('test') 68 55 %% set the input elements needed on the GUI series when the function is selected in the menu ActionName or InputTable refreshed 69 56 if isstruct(Param) && isequal(Param.Action.RUN,0) … … 137 124 %% define the directory for result file (with path=RootPath{1}) 138 125 OutputDir=[Param.OutputSubDir Param.OutputDirExt];% subdirectory for output files 126 % 127 % if ~isfield(Param,'InputFields') 128 % Param.InputFields.FieldName=''; 129 % end 139 130 140 131 %% calibration data and timing: read the ImaDoc files … … 169 160 return 170 161 end 171 ObjectData=Param.ProjObject; % Object attached to the GUI series162 ObjectData=Param.ProjObject; 172 163 xI=ObjectData.RangeX(1):ObjectData.DX:ObjectData.RangeX(2); 173 164 yI=ObjectData.RangeY(1):ObjectData.DY:ObjectData.RangeY(2); … … 177 168 W=zeros(size(XI,1),size(XI,2)); 178 169 179 %% %% MAIN LOOP ON FIELDS %%%%%170 %% MAIN LOOP ON FIELDS 180 171 warning off 181 172 … … 184 175 CheckOverwrite=Param.CheckOverwrite; 185 176 end 186 for field_index=1:NbField 187 188 update_waitbar(WaitbarHandle,field_index/NbField)% waitbar to visualise progress (in RUN mode) 189 190 %% generating the name of the output file for the merged field 191 i1=i1_series{1}(field_index); 177 for index=1:NbField 178 179 update_waitbar(WaitbarHandle,index/NbField) 180 181 182 183 184 %% generating the name of the merged field 185 i1=i1_series{1}(index); 192 186 if ~isempty(i2_series{end}) 193 i2=i2_series{end}( field_index);187 i2=i2_series{end}(index); 194 188 else 195 189 i2=i1; … … 198 192 j2=1; 199 193 if ~isempty(j1_series{1}) 200 j1=j1_series{1}( field_index);194 j1=j1_series{1}(index); 201 195 if ~isempty(j2_series{end}) 202 j2=j2_series{end}( field_index);196 j2=j2_series{end}(index); 203 197 else 204 198 j2=j1; … … 207 201 OutputFile=fullfile_uvmat(RootPath{1},OutputDir,RootFile{1},'.nc','_12',i1,i2,j1,j2); 208 202 209 %% program stop if requested on the GUI 203 %% 204 205 210 206 if ~isempty(RUNHandle) && ~strcmp(get(RUNHandle,'BusyAction'),'queue') 211 207 disp('program stopped by user') … … 213 209 end 214 210 215 %% check existence of the output file 216 if (~CheckOverwrite && exist(OutputFile,'file')) 217 disp('existing output file already exists, skip to next field') 218 continue% skip iteration if the mode overwrite is desactivated and the result file already exists 219 end 220 211 if (~CheckOverwrite && exist(OutputFile,'file')) 212 disp('existing output file already exists, skip to next field') 213 continue% skip iteration if the mode overwrite is desactivated and the result file already exists 214 end 215 221 216 %%%%%%%%%%%%%%%% loop on views (input lines) %%%%%%%%%%%%%%%% 222 217 Data=cell(1,NbView);%initiate the set Data 218 timeread=zeros(1,NbView); 223 219 224 220 %get Xphys,Yphys,Zphys from 1 or 2 stereo folders. Positions are taken 225 221 %at the middle between to time step 226 ZItemp=zeros(size(XI,1),size(XI,2),2); 227 if field_index==1 222 clear ZItemp 223 ZItemp=zeros(size(XI,1),size(XI,2),2); 224 225 if index==1 228 226 first_img=i1_series{1,1}(1,1); %id of the first image of the series 229 end 230 231 idtemp=0; 232 % get the surface shape corresponding to the PIV measurements 233 for indextemp=field_index:field_index+1;%TODO: generalise to field index intervals>1 for PIV 234 idtemp=idtemp+1; 235 if NbView==3 % if there is only 1 stereo folder (2 cameras only), extract directly Xphys,Yphys and Zphys 236 InputFile_3=fullfile(Param.InputTable{3,1},Param.InputTable{3,2},[Param.InputTable{3,3} '_' int2str(first_img+indextemp1) '.nc']); 237 % Data{1}: =raw PIV camera 1 only 238 % Data{2}: =raw PIV camera 2 only 239 % Data{3}: =correlation between camera 1 and 2 240 [Data{3},tild,errormsg] = nc2struct(InputFile_3);%read input file 241 if exist('Data{3}.Civ3_FF','var') % FF is present, remove wrong vector 242 temp=find(Data{3}.Civ3_FF==0); 243 Zphys=Data{3}.Zphys(temp); 244 Yphys=Data{3}.Yphys(temp); 245 Xphys=Data{3}.Xphys(temp); 246 else 247 Zphys=Data{3}.Zphys; 248 Yphys=Data{3}.Yphys; 249 Xphys=Data{3}.Xphys; 227 end 228 229 idtemp=0; 230 for indextemp=index:index+1; 231 idtemp=idtemp+1; 232 if NbView==3 % if there is only 1 stereo folder, extract directly Xphys,Yphys and Zphys 233 234 235 236 [Data{3},tild,errormsg] = nc2struct([Param.InputTable{3,1},'/',Param.InputTable{3,2},'/',Param.InputTable{3,3},'_',int2str(first_img+indextemp1),'.nc']); 237 238 if exist('Data{3}.Civ3_FF','var') % FF is present, remove wrong vector 239 temp=find(Data{3}.Civ3_FF==0); 240 Zphys=Data{3}.Zphys(temp); 241 Yphys=Data{3}.Yphys(temp); 242 Xphys=Data{3}.Xphys(temp); 243 else 244 Zphys=Data{3}.Zphys; 245 Yphys=Data{3}.Yphys; 246 Xphys=Data{3}.Xphys; 247 end 248 249 250 251 elseif NbView==4 % is there is 2 stereo folders, get global U and V and compute Zphys 252 253 254 %test if the seconde camera is the same for both folder 255 for i=3:4 256 indpt(i)=strfind(Param.InputTable{i,2},'.'); % indice of the "." is the folder name 1 257 indline(i)=strfind(Param.InputTable{i,2},''); % indice of the "" is the folder name1 258 camname{i}=Param.InputTable{i,2}(indline(i)+1:indpt(i)1);% extract the second camera name 259 end 260 261 if strcmp(camname{3},camname{4})==0 262 disp_uvmat('ERROR','The 2 stereo folders should have the same camera for the second position',checkrun) 263 return 264 end 265 266 267 268 [Data{3},tild,errormsg] = nc2struct([Param.InputTable{3,1},'/',Param.InputTable{3,2},'/',Param.InputTable{3,3},'_',int2str(first_img+indextemp1),'.nc']); 269 270 if exist('Data{3}.Civ3_FF','var') % if FF is present, remove wrong vector 271 temp=find(Data{3}.Civ3_FF==0); 272 Xmid3=Data{3}.Xmid(temp); 273 Ymid3=Data{3}.Ymid(temp); 274 U3=Data{3}.Uphys(temp); 275 V3=Data{3}.Vphys(temp); 276 else 277 Xmid3=Data{3}.Xmid; 278 Ymid3=Data{3}.Ymid; 279 U3=Data{3}.Uphys; 280 V3=Data{3}.Vphys; 281 end 282 %temporary gridd of merging the 2 stereos datas 283 [xq,yq] = meshgrid(min(Xmid3+(U3)/2):(max(Xmid3+(U3)/2)min(Xmid3+(U3)/2))/128:max(Xmid3+(U3)/2),min(Ymid3+(V3)/2):(max(Ymid3+(V3)/2)min(Ymid3+(V3)/2))/128:max(Ymid3+(V3)/2)); 284 285 %1st folder : interpolate the first camera (Dalsa1) points on the second (common) camera 286 %(Dalsa 3) 287 x3Q=griddata(Xmid3+(U3)/2,Ymid3+(V3)/2,Xmid3(U3)/2,xq,yq); 288 y3Q=griddata(Xmid3+(U3)/2,Ymid3+(V3)/2,Ymid3(V3)/2,xq,yq); 289 290 291 292 [Data{4},tild,errormsg] = nc2struct([Param.InputTable{4,1},'/',Param.InputTable{4,2},'/',Param.InputTable{4,3},'_',int2str(first_img+indextemp1),'.nc']); 293 if exist('Data{4}.Civ3_FF','var') % if FF is present, remove wrong vector 294 temp=find(Data{4}.Civ3_FF==0); 295 Xmid4=Data{4}.Xmid(temp); 296 Ymid4=Data{4}.Ymid(temp); 297 U4=Data{4}.Uphys(temp); 298 V4=Data{4}.Vphys(temp); 299 else 300 Xmid4=Data{4}.Xmid; 301 Ymid4=Data{4}.Ymid; 302 U4=Data{4}.Uphys; 303 V4=Data{4}.Vphys; 304 end 305 306 %2nd folder :interpolate the first camera (Dalsa2) points on the second (common) camera 307 %(Dalsa 3) 308 x4Q=griddata(Xmid4+(U4)/2,Ymid4+(V4)/2,Xmid4(U4)/2,xq,yq); 309 y4Q=griddata(Xmid4+(U4)/2,Ymid4+(V4)/2,Ymid4(V4)/2,xq,yq); 310 311 xmid=reshape((x4Q+x3Q)/2,length(xq(:,1)).*length(xq(1,:)),1); 312 ymid=reshape((y4Q+y3Q)/2,length(yq(:,1)).*length(yq(1,:)),1); 313 u=reshape(x4Qx3Q,length(xq(:,1)).*length(xq(1,:)),1); 314 v=reshape(y4Qy3Q,length(yq(:,1)).*length(yq(1,:)),1); 315 316 317 [Zphys,Xphys,Yphys,error]=shift2z(xmid, ymid, u, v,XmlData); %get Xphy,Yphy and Zphys 318 %remove NaN 319 tempNaN=isnan(Zphys);tempind=find(tempNaN==1); 320 Zphys(tempind)=[]; 321 Xphys(tempind)=[]; 322 Yphys(tempind)=[]; 323 error(tempind)=[]; 324 325 end 326 327 if NbView>2 328 ZItemp(:,:,idtemp)=griddata(Xphys,Yphys,Zphys,XI,YI); %interpolation on the choosen gridd 250 329 end 251 252 elseif NbView==4 % is there is 2 stereo folders, get global U and V and compute Zphys 253 % Data{1}: =raw PIV camera 1 only (left) 254 % Data{2}: =raw PIV camera 2 only (right) (no PIV done with middle camera) 255 % Data{3}: =corelation between camera 1 and 3 (left and middle) 256 % Data{4}: =corelation between camera 2 and 3 (right and middle) 257 % test if the second camera (3) is the same for both folders 258 for i=3:4 259 indpt(i)=strfind(Param.InputTable{i,2},'.'); % indice of the "." is the folder name 1 260 indline(i)=strfind(Param.InputTable{i,2},''); % indice of the "" is the folder name1 261 camname{i}=Param.InputTable{i,2}(indline(i)+1:indpt(i)1);% extract the second camera name 262 end 263 if strcmp(camname{3},camname{4})==0 264 disp_uvmat('ERROR','The 2 stereo folders should have the same camera for the second position',checkrun) 265 return 266 end 267 [Data{3},tild,errormsg] = nc2struct(InputFile_3); 268 if exist('Data{3}.Civ3_FF','var') % if FF is present, remove wrong vector 269 temp=find(Data{3}.Civ3_FF==0); 270 Xmid3=Data{3}.Xmid(temp); 271 Ymid3=Data{3}.Ymid(temp); 272 U3=Data{3}.Uphys(temp); 273 V3=Data{3}.Vphys(temp); 274 else 275 Xmid3=Data{3}.Xmid; 276 Ymid3=Data{3}.Ymid; 277 U3=Data{3}.Uphys; 278 V3=Data{3}.Vphys; 279 end 280 %temporary grid of merging the 2 stereos data 281 [xq,yq] = meshgrid(min(Xmid3+(U3)/2):(max(Xmid3+(U3)/2)min(Xmid3+(U3)/2))/128:max(Xmid3+(U3)/2),min(Ymid3+(V3)/2):(max(Ymid3+(V3)/2)min(Ymid3+(V3)/2))/128:max(Ymid3+(V3)/2)); 282 283 %1st folder : interpolate the first camera points on the second (common) camera 284 %(Dalsa 3) 285 x3Q=griddata(Xmid3+(U3)/2,Ymid3+(V3)/2,Xmid3(U3)/2,xq,yq); 286 y3Q=griddata(Xmid3+(U3)/2,Ymid3+(V3)/2,Ymid3(V3)/2,xq,yq); 287 288 InputFile_4=fullfile(Param.InputTable{4,1},Param.InputTable{4,2},[Param.InputTable{4,3} '_' int2str(first_img+indextemp1) '.nc']); 289 [Data{4},tild,errormsg] = nc2struct(InputFile_4); 290 if exist('Data{4}.Civ3_FF','var') % if FF is present, remove wrong vector 291 temp=find(Data{4}.Civ3_FF==0); 292 Xmid4=Data{4}.Xmid(temp); 293 Ymid4=Data{4}.Ymid(temp); 294 U4=Data{4}.Uphys(temp); 295 V4=Data{4}.Vphys(temp); 296 else 297 Xmid4=Data{4}.Xmid; 298 Ymid4=Data{4}.Ymid; 299 U4=Data{4}.Uphys; 300 V4=Data{4}.Vphys; 301 end 302 %2nd folder :interpolate the first camera (Dalsa2) points on the second (common) camera 303 %(Dalsa 3) 304 x4Q=griddata(Xmid4+(U4)/2,Ymid4+(V4)/2,Xmid4(U4)/2,xq,yq); 305 y4Q=griddata(Xmid4+(U4)/2,Ymid4+(V4)/2,Ymid4(V4)/2,xq,yq); 306 307 %add the displacements of the two camera pairs 308 xmid=reshape((x4Q+x3Q)/2,length(xq(:,1)).*length(xq(1,:)),1); 309 ymid=reshape((y4Q+y3Q)/2,length(yq(:,1)).*length(yq(1,:)),1); 310 u=reshape(x4Qx3Q,length(xq(:,1)).*length(xq(1,:)),1); 311 v=reshape(y4Qy3Q,length(yq(:,1)).*length(yq(1,:)),1); 312 313 % get the surface z(x,y) from the combined displacement 314 [Zphys,Xphys,Yphys,error]=shift2z(xmid, ymid, u, v,XmlData); %get Xphy,Yphy and Zphys 315 %remove NaN 316 tempNaN=isnan(Zphys);tempind=find(tempNaN==1); 317 Zphys(tempind)=[]; 318 Xphys(tempind)=[]; 319 Yphys(tempind)=[]; 320 error(tempind)=[]; 321 322 end 323 if NbView>2 324 ZItemp(:,:,idtemp)=griddata(Xphys,Yphys,Zphys,XI,YI); %interpolation on the choosen grid 325 end 326 327 end 328 ZI=mean(ZItemp,3); %mean between two the two times used for surface measurement 330 331 end 332 ZI=mean(ZItemp,3); %mean between two the two time step 329 333 Vtest=ZItemp(:,:,2)ZItemp(:,:,1); 330 334 … … 332 336 [Xb,Yb]=px_XYZ(XmlData{2}.GeometryCalib,XI,YI,ZI);% set of image coordinates on view b 333 337 334 338 335 339 for iview=1:2 336 %% reading PIVinput file(s)337 [Data{iview},tild,errormsg]=read_civdata(filecell{iview, field_index},{'vec(U,V)'},Param.InputFields.VelType);340 %% reading input file(s) 341 [Data{iview},tild,errormsg]=read_civdata(filecell{iview,index},{'vec(U,V)'},'*'); 338 342 if ~isempty(errormsg) 339 343 disp_uvmat('ERROR',['ERROR in civ2vel_3C/read_field/' errormsg],checkrun) … … 354 358 end 355 359 end 356 %remove false vectors360 %remove wrong vector 357 361 temp=find(Data{1}.FF==0); 358 362 X1=Data{1}.X(temp); … … 364 368 Va=griddata(X1,Y1,V1,Xa,Ya); 365 369 366 [Ua,Va,Xa,Ya]=Ud2U(XmlData{1}.GeometryCalib,Xa,Ya,Ua,Va); % convert Xd data to X 370 [Ua,Va,Xa,Ya]=Ud2U(XmlData{1}.GeometryCalib,Xa,Ya,Ua,Va); % convert Xd data to X 367 371 [A]=get_coeff(XmlData{1}.GeometryCalib,Xa,Ya,XI,YI,ZI); %get coef A~ 368 372 369 %remove false vectors373 %remove wrong vector 370 374 temp=find(Data{2}.FF==0); 371 375 X2=Data{2}.X(temp); … … 375 379 Ub=griddata(X2,Y2,U2,Xb,Yb); 376 380 Vb=griddata(X2,Y2,V2,Xb,Yb); 377 378 [Ub,Vb,Xb,Yb]=Ud2U(XmlData{2}.GeometryCalib,Xb,Yb,Ub,Vb); % convert Xd data to X 381 382 [Ub,Vb,Xb,Yb]=Ud2U(XmlData{2}.GeometryCalib,Xb,Yb,Ub,Vb); % convert Xd data to X 379 383 [B]=get_coeff(XmlData{2}.GeometryCalib,Xb,Yb,XI,YI,ZI); %get coef B~ 380 384 381 385 382 386 % System to solve 383 387 S=ones(size(XI,1),size(XI,2),3); 384 388 D=ones(size(XI,1),size(XI,2),3,3); 385 389 386 390 S(:,:,1)=A(:,:,1,1).*Ua+A(:,:,2,1).*Va+B(:,:,1,1).*Ub+B(:,:,2,1).*Vb; 387 391 S(:,:,2)=A(:,:,1,2).*Ua+A(:,:,2,2).*Va+B(:,:,1,2).*Ub+B(:,:,2,2).*Vb; … … 403 407 W(indj,indi)=dxyz(3); 404 408 end 405 end 409 end 406 410 Error=zeros(size(XI,1),size(XI,2),4); 407 411 Error(:,:,1)=A(:,:,1,1).*U+A(:,:,1,2).*V+A(:,:,1,3).*WUa; … … 410 414 Error(:,:,4)=B(:,:,2,1).*U+B(:,:,2,2).*V+B(:,:,2,3).*WVb; 411 415 416 417 418 419 420 412 421 %% recording the merged field 413 if field_index==1% initiate the structure at first index422 if index==1% initiate the structure at first index 414 423 MergeData.ListGlobalAttribute={'Conventions','Time','Dt'}; 415 424 MergeData.Conventions='uvmat'; … … 418 427 MergeData.ListVarName={'coord_x','coord_y','Z','U','V','W','Error'}; 419 428 MergeData.VarDimName={'coord_x','coord_y',{'coord_y','coord_x'},{'coord_y','coord_x'}... 420 {'coord_y','coord_x'},{'coord_y','coord_x'},{'coord_y','coord_x'}};429 {'coord_y','coord_x'},{'coord_y','coord_x'},{'coord_y','coord_x'}}; 421 430 MergeData.coord_x=xI; 422 431 MergeData.coord_y=yI; … … 427 436 MergeData.Z=ZI; 428 437 429 430 438 % mfx=(XmlData{1}.GeometryCalib.fx_fy(1)+XmlData{2}.GeometryCalib.fx_fy(1))/2; 439 % mfy=(XmlData{1}.GeometryCalib.fx_fy(2)+XmlData{2}.GeometryCalib.fx_fy(2))/2; 431 440 MergeData.Error=0.5*sqrt(sum(Error.^2,3)); 432 441 errormsg=struct2nc(OutputFile,MergeData);%save result file … … 451 460 A(:,:,2,3)=(R(6)R(9)*Y)./T; 452 461 453 function [U,V,X,Y]=Ud2U(Calib,Xd,Yd,Ud,Vd) 454 % convert image coordinates to view angles, after removal of quadratic distorsion 455 % input in pixel, output in radians 462 function [U,V,X,Y]=Ud2U(Calib,Xd,Yd,Ud,Vd) % convert Xd to X and Ud to U 463 456 464 X1d=XdUd/2; 457 465 X2d=Xd+Ud/2; … … 474 482 z=0; 475 483 error=0; 484 476 485 477 486 %% first image
Note: See TracChangeset
for help on using the changeset viewer.