Changeset 954 for trunk/src/series/civ2vel_3C.m
 Timestamp:
 Jun 22, 2016, 1:36:50 PM (7 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

trunk/src/series/civ2vel_3C.m
r927 r954 1 %'civ2vel_3C': combine velocity fields from two cameras to get three velocity components 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 % 2 14 % 3 15 % function ParamOut=civ2vel_3C(Param) … … 7 19 % 8 20 %INPUT: 21 % 9 22 % In run mode, the input parameters are given as a Matlab structure Param copied from the GUI series. 10 23 % In batch mode, Param is the name of the corresponding xml file containing the same information … … 15 28 % see the current structure Param) 16 29 % .InputTable: cell of input file names, (several lines for multiple input) 17 % each line decomposed as {RootPath,SubDir,Rootfile,NomType,Extension} 30 % each line decomposed as {RootPath,SubDir,Rootfile,NomType,Extension}. 31 % 3 or 4 lines used as described above 18 32 % .OutputSubDir: name of the subdirectory for data outputs 19 33 % .OutputDirExt: directory extension for data outputs … … 22 36 % .ActionExt: fct extension ('.m', Matlab fct, '.sh', compiled Matlab fct 23 37 % .RUN =0 for GUI input, =1 for function activation 24 % .RunMode='local','background', 'cluster': type of function use 25 % 38 % .RunMode='local','background', 'cluster': type of function use 26 39 % .IndexRange: set the file or frame indices on which the action must be performed 27 40 % .InputFields: sub structure describing the input fields withfields … … 52 65 53 66 function ParamOut=civ2vel_3C(Param) 54 disp('test') 67 55 68 %% set the input elements needed on the GUI series when the function is selected in the menu ActionName or InputTable refreshed 56 69 if isstruct(Param) && isequal(Param.Action.RUN,0) … … 124 137 %% define the directory for result file (with path=RootPath{1}) 125 138 OutputDir=[Param.OutputSubDir Param.OutputDirExt];% subdirectory for output files 126 %127 % if ~isfield(Param,'InputFields')128 % Param.InputFields.FieldName='';129 % end130 139 131 140 %% calibration data and timing: read the ImaDoc files … … 160 169 return 161 170 end 162 ObjectData=Param.ProjObject; 171 ObjectData=Param.ProjObject;% Object attached to the GUI series 163 172 xI=ObjectData.RangeX(1):ObjectData.DX:ObjectData.RangeX(2); 164 173 yI=ObjectData.RangeY(1):ObjectData.DY:ObjectData.RangeY(2); … … 168 177 W=zeros(size(XI,1),size(XI,2)); 169 178 170 %% MAIN LOOP ON FIELDS179 %%%% MAIN LOOP ON FIELDS %%%%% 171 180 warning off 172 181 … … 175 184 CheckOverwrite=Param.CheckOverwrite; 176 185 end 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); 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); 186 192 if ~isempty(i2_series{end}) 187 i2=i2_series{end}( index);193 i2=i2_series{end}(field_index); 188 194 else 189 195 i2=i1; … … 192 198 j2=1; 193 199 if ~isempty(j1_series{1}) 194 j1=j1_series{1}( index);200 j1=j1_series{1}(field_index); 195 201 if ~isempty(j2_series{end}) 196 j2=j2_series{end}( index);202 j2=j2_series{end}(field_index); 197 203 else 198 204 j2=j1; … … 201 207 OutputFile=fullfile_uvmat(RootPath{1},OutputDir,RootFile{1},'.nc','_12',i1,i2,j1,j2); 202 208 203 %% 204 205 209 %% program stop if requested on the GUI 206 210 if ~isempty(RUNHandle) && ~strcmp(get(RUNHandle,'BusyAction'),'queue') 207 211 disp('program stopped by user') … … 209 213 end 210 214 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 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 216 221 %%%%%%%%%%%%%%%% loop on views (input lines) %%%%%%%%%%%%%%%% 217 222 Data=cell(1,NbView);%initiate the set Data 218 timeread=zeros(1,NbView);219 223 220 224 %get Xphys,Yphys,Zphys from 1 or 2 stereo folders. Positions are taken 221 225 %at the middle between to time step 222 clear ZItemp 223 ZItemp=zeros(size(XI,1),size(XI,2),2); 224 225 if index==1 226 ZItemp=zeros(size(XI,1),size(XI,2),2); 227 if field_index==1 226 228 first_img=i1_series{1,1}(1,1); %id of the first image of the series 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; 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 InputFile_3=fullfile(Param.InputTable{3,1},Param.InputTable{3,2},[Param.InputTable{3,3} '_' int2str(first_img+indextemp1) '.nc']); 236 if NbView==3 % if there is only 1 stereo folder (2 cameras only), extract directly Xphys,Yphys and Zphys 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; 250 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 247 322 end 248 323 324 ZItemp(:,:,idtemp)=griddata(Xphys,Yphys,Zphys,XI,YI); %interpolation on the choosen grid 249 325 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 328 329 330 ZItemp(:,:,idtemp)=griddata(Xphys,Yphys,Zphys,XI,YI); %interpolation on the choosen gridd 331 332 end 333 ZI=mean(ZItemp,3); %mean between two the two time step 326 end 327 ZI=mean(ZItemp,3); %mean between two the two times used for surface measurement 334 328 Vtest=ZItemp(:,:,2)ZItemp(:,:,1); 335 329 … … 337 331 [Xb,Yb]=px_XYZ(XmlData{2}.GeometryCalib,XI,YI,ZI);% set of image coordinates on view b 338 332 339 333 340 334 for iview=1:2 341 %% reading input file(s)342 [Data{iview},tild,errormsg]=read_civdata(filecell{iview, index},{'vec(U,V)'},'*');335 %% reading PIV input file(s) 336 [Data{iview},tild,errormsg]=read_civdata(filecell{iview,field_index},{'vec(U,V)'},'*'); 343 337 if ~isempty(errormsg) 344 338 disp_uvmat('ERROR',['ERROR in civ2vel_3C/read_field/' errormsg],checkrun) … … 359 353 end 360 354 end 361 %remove wrong vector355 %remove false vectors 362 356 temp=find(Data{1}.FF==0); 363 357 X1=Data{1}.X(temp); … … 369 363 Va=griddata(X1,Y1,V1,Xa,Ya); 370 364 371 % [Ua,Va,Xa,Ya]=Ud2U(XmlData{1}.GeometryCalib,Xa,Ya,Ua,Va); % convert Xd data to X 365 [Ua,Va,Xa,Ya]=Ud2U(XmlData{1}.GeometryCalib,Xa,Ya,Ua,Va); % convert Xd data to X 372 366 [A]=get_coeff(XmlData{1}.GeometryCalib,Xa,Ya,XI,YI,ZI); %get coef A~ 373 367 374 %remove wrong vector368 %remove false vectors 375 369 temp=find(Data{2}.FF==0); 376 370 X2=Data{2}.X(temp); … … 380 374 Ub=griddata(X2,Y2,U2,Xb,Yb); 381 375 Vb=griddata(X2,Y2,V2,Xb,Yb); 382 383 % [Ub,Vb,Xb,Yb]=Ud2U(XmlData{2}.GeometryCalib,Xb,Yb,Ub,Vb); % convert Xd data to X 376 377 [Ub,Vb,Xb,Yb]=Ud2U(XmlData{2}.GeometryCalib,Xb,Yb,Ub,Vb); % convert Xd data to X 384 378 [B]=get_coeff(XmlData{2}.GeometryCalib,Xb,Yb,XI,YI,ZI); %get coef B~ 385 379 386 380 387 381 % System to solve 388 382 S=ones(size(XI,1),size(XI,2),3); 389 383 D=ones(size(XI,1),size(XI,2),3,3); 390 384 391 385 S(:,:,1)=A(:,:,1,1).*Ua+A(:,:,2,1).*Va+B(:,:,1,1).*Ub+B(:,:,2,1).*Vb; 392 386 S(:,:,2)=A(:,:,1,2).*Ua+A(:,:,2,2).*Va+B(:,:,1,2).*Ub+B(:,:,2,2).*Vb; … … 408 402 W(indj,indi)=dxyz(3); 409 403 end 410 end 404 end 411 405 Error=zeros(size(XI,1),size(XI,2),4); 412 406 Error(:,:,1)=A(:,:,1,1).*U+A(:,:,1,2).*V+A(:,:,1,3).*WUa; … … 415 409 Error(:,:,4)=B(:,:,2,1).*U+B(:,:,2,2).*V+B(:,:,2,3).*WVb; 416 410 417 418 419 420 421 422 411 %% recording the merged field 423 if index==1% initiate the structure at first index412 if field_index==1% initiate the structure at first index 424 413 MergeData.ListGlobalAttribute={'Conventions','Time','Dt'}; 425 414 MergeData.Conventions='uvmat'; … … 428 417 MergeData.ListVarName={'coord_x','coord_y','Z','U','V','W','Error'}; 429 418 MergeData.VarDimName={'coord_x','coord_y',{'coord_y','coord_x'},{'coord_y','coord_x'}... 430 419 {'coord_y','coord_x'},{'coord_y','coord_x'},{'coord_y','coord_x'}}; 431 420 MergeData.coord_x=xI; 432 421 MergeData.coord_y=yI; … … 437 426 MergeData.Z=ZI; 438 427 439 % mfx=(XmlData{1}.GeometryCalib.fx_fy(1)+XmlData{2}.GeometryCalib.fx_fy(1))/2;440 % mfy=(XmlData{1}.GeometryCalib.fx_fy(2)+XmlData{2}.GeometryCalib.fx_fy(2))/2;428 % mfx=(XmlData{1}.GeometryCalib.fx_fy(1)+XmlData{2}.GeometryCalib.fx_fy(1))/2; 429 % mfy=(XmlData{1}.GeometryCalib.fx_fy(2)+XmlData{2}.GeometryCalib.fx_fy(2))/2; 441 430 MergeData.Error=0.5*sqrt(sum(Error.^2,3)); 442 431 errormsg=struct2nc(OutputFile,MergeData);%save result file … … 461 450 A(:,:,2,3)=(R(6)R(9)*Y)./T; 462 451 463 function [U,V,X,Y]=Ud2U(Calib,Xd,Yd,Ud,Vd) % convert Xd to X and Ud to U 464 452 function [U,V,X,Y]=Ud2U(Calib,Xd,Yd,Ud,Vd) 453 % convert image coordinates to view angles, after removal of quadratic distorsion 454 % input in pixel, output in radians 465 455 X1d=XdUd/2; 466 456 X2d=Xd+Ud/2; … … 483 473 z=0; 484 474 error=0; 485 486 475 487 476 %% first image
Note: See TracChangeset
for help on using the changeset viewer.