Changeset 1144 for trunk/src/series/civ2vel_3C.m
 Timestamp:
 May 13, 2024, 9:49:09 PM (6 months ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

trunk/src/series/civ2vel_3C.m
r1141 r1144 58 58 ParamOut.WholeIndexRange='off';% prescribes the file index ranges from min to max (options 'off'/'on', 'off' by default) 59 59 ParamOut.NbSlice='off'; %nbre of slices ('off' by default) 60 ParamOut.VelType='o ne';% menu for selecting the velocity type (options 'off'/'one'/'two', 'off' by default)60 ParamOut.VelType='off';% menu for selecting the velocity type (options 'off'/'one'/'two', 'off' by default) 61 61 ParamOut.FieldName='off';% menu for selecting the field (s) in the input file(options 'off'/'one'/'two', 'off' by default) 62 62 ParamOut.FieldTransform = 'off';%use the phys transform function without choice … … 108 108 hdisp=disp_uvmat('WAITING...','checking the file series',checkrun); 109 109 [filecell,i1_series,i2_series,j1_series,j2_series]=get_file_series(Param); 110 if ~isempty(hdisp),delete(hdisp),end ;110 if ~isempty(hdisp),delete(hdisp),end 111 111 %%%%%%%%%%%% 112 112 % The cell array filecell is the list of input file names, while … … 124 124 %% define the directory for result file (with path=RootPath{1}) 125 125 OutputDir=[Param.OutputSubDir Param.OutputDirExt];% subdirectory for output files 126 %127 % if ~isfield(Param,'InputFields')128 % Param.InputFields.FieldName='';129 % end130 126 131 127 %% calibration data and timing: read the ImaDoc files … … 175 171 CheckOverwrite=Param.CheckOverwrite; 176 172 end 173 177 174 for index=1:NbField 178 175 179 180 181 %% generating the name of the merged field 176 %% generating the name of the merged field 182 177 i1=i1_series{1}(index); 183 178 if ~isempty(i2_series{end}) … … 198 193 OutputFile=fullfile_uvmat(RootPath{1},OutputDir,RootFile{1},'.nc','_12',i1,i2,j1,j2); 199 194 200 %%201 202 203 195 if ~isempty(RUNHandle) && ~strcmp(get(RUNHandle,'BusyAction'),'queue') 204 196 disp('program stopped by user') … … 206 198 end 207 199 208 if (~CheckOverwrite && exist(OutputFile,'file'))209 210 211 end212 200 if (~CheckOverwrite && exist(OutputFile,'file')) 201 disp('existing output file already exists, skip to next field') 202 continue% skip iteration if the mode overwrite is desactivated and the result file already exists 203 end 204 213 205 %%%%%%%%%%%%%%%% loop on views (input lines) %%%%%%%%%%%%%%%% 214 206 Data=cell(1,NbView);%initiate the set Data … … 217 209 %get Xphys,Yphys,Zphys from 1 or 2 stereo folders. Positions are taken 218 210 %at the middle between to time step 219 clear ZItemp 220 ZItemp=zeros(size(XI,1),size(XI,2),2); 221 222 if index==1 211 clear ZItemp 212 ZItemp=zeros(size(XI,1),size(XI,2),2); 213 CheckZ=0; 214 215 if index==1 223 216 first_img=i1_series{1,1}(1,1); %id of the first image of the series 224 end 225 226 idtemp=0; 227 for indextemp=index:index+1; 228 idtemp=idtemp+1; 229 if NbView==3 % if there is only 1 stereo folder, extract directly Xphys,Yphys and Zphys 230 217 end 218 219 idtemp=0; 220 for indextemp=index:index+1 221 idtemp=idtemp+1; 222 if NbView==3 % if there is only 1 stereo folder, extract directly Xphys,Yphys and Zphys 223 224 [Data{3},tild,errormsg] = nc2struct([Param.InputTable{3,1},'/',Param.InputTable{3,2},'/',Param.InputTable{3,3},'_',int2str(first_img+indextemp1),'.nc']); 225 226 if exist('Data{3}.Civ3_FF','var') % FF is present, remove wrong vector 227 temp=find(Data{3}.Civ3_FF==0); 228 Zphys=Data{3}.Zphys(temp); 229 Yphys=Data{3}.Yphys(temp); 230 Xphys=Data{3}.Xphys(temp); 231 else 232 Zphys=Data{3}.Zphys; 233 Yphys=Data{3}.Yphys; 234 Xphys=Data{3}.Xphys; 235 end 236 237 elseif NbView==4 % is there is 2 stereo folders, get global U and V and compute Zphys 238 239 %test if the seconde camera is the same for both folder 240 for i=3:4 241 indpt(i)=strfind(Param.InputTable{i,2},'.'); % indice of the "." is the folder name 1 242 indline(i)=strfind(Param.InputTable{i,2},''); % indice of the "" is the folder name1 243 camname{i}=Param.InputTable{i,2}(indline(i)+1:indpt(i)1);% extract the second camera name 244 end 245 246 if strcmp(camname{3},camname{4})==0 247 disp_uvmat('ERROR','The 2 stereo folders should have the same camera for the second position',checkrun) 248 return 249 end 250 251 [Data{3},tild,errormsg] = nc2struct([Param.InputTable{3,1},'/',Param.InputTable{3,2},'/',Param.InputTable{3,3},'_',int2str(first_img+indextemp1),'.nc']); 252 253 if exist('Data{3}.Civ3_FF','var') % if FF is present, remove wrong vector 254 temp=find(Data{3}.Civ3_FF==0); 255 Xmid3=Data{3}.Xmid(temp); 256 Ymid3=Data{3}.Ymid(temp); 257 U3=Data{3}.Uphys(temp); 258 V3=Data{3}.Vphys(temp); 259 else 260 Xmid3=Data{3}.Xmid; 261 Ymid3=Data{3}.Ymid; 262 U3=Data{3}.Uphys; 263 V3=Data{3}.Vphys; 264 end 265 %temporary gridd of merging the 2 stereos datas 266 [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)); 267 268 %1st folder : interpolate the first camera (Dalsa1) points on the second (common) camera 269 %(Dalsa 3) 270 x3Q=griddata(Xmid3+(U3)/2,Ymid3+(V3)/2,Xmid3(U3)/2,xq,yq); 271 y3Q=griddata(Xmid3+(U3)/2,Ymid3+(V3)/2,Ymid3(V3)/2,xq,yq); 272 273 [Data{4},tild,errormsg] = nc2struct([Param.InputTable{4,1},'/',Param.InputTable{4,2},'/',Param.InputTable{4,3},'_',int2str(first_img+indextemp1),'.nc']); 274 if exist('Data{4}.Civ3_FF','var') % if FF is present, remove wrong vector 275 temp=find(Data{4}.Civ3_FF==0); 276 Xmid4=Data{4}.Xmid(temp); 277 Ymid4=Data{4}.Ymid(temp); 278 U4=Data{4}.Uphys(temp); 279 V4=Data{4}.Vphys(temp); 280 else 281 Xmid4=Data{4}.Xmid; 282 Ymid4=Data{4}.Ymid; 283 U4=Data{4}.Uphys; 284 V4=Data{4}.Vphys; 285 end 286 287 %2nd folder :interpolate the first camera (Dalsa2) points on the second (common) camera 288 %(Dalsa 3) 289 x4Q=griddata(Xmid4+(U4)/2,Ymid4+(V4)/2,Xmid4(U4)/2,xq,yq); 290 y4Q=griddata(Xmid4+(U4)/2,Ymid4+(V4)/2,Ymid4(V4)/2,xq,yq); 291 292 xmid=reshape((x4Q+x3Q)/2,length(xq(:,1)).*length(xq(1,:)),1); 293 ymid=reshape((y4Q+y3Q)/2,length(yq(:,1)).*length(yq(1,:)),1); 294 u=reshape(x4Qx3Q,length(xq(:,1)).*length(xq(1,:)),1); 295 v=reshape(y4Qy3Q,length(yq(:,1)).*length(yq(1,:)),1); 296 297 298 [Zphys,Xphys,Yphys,error]=shift2z(xmid, ymid, u, v,XmlData); %get Xphy,Yphy and Zphys 299 %remove NaN 300 tempNaN=isnan(Zphys);tempind=find(tempNaN==1); 301 Zphys(tempind)=[]; 302 Xphys(tempind)=[]; 303 Yphys(tempind)=[]; 304 end 231 305 232 233 [Data{3},tild,errormsg] = nc2struct([Param.InputTable{3,1},'/',Param.InputTable{3,2},'/',Param.InputTable{3,3},'_',int2str(first_img+indextemp1),'.nc']); 234 235 if exist('Data{3}.Civ3_FF','var') % FF is present, remove wrong vector 236 temp=find(Data{3}.Civ3_FF==0); 237 Zphys=Data{3}.Zphys(temp); 238 Yphys=Data{3}.Yphys(temp); 239 Xphys=Data{3}.Xphys(temp); 240 else 241 Zphys=Data{3}.Zphys; 242 Yphys=Data{3}.Yphys; 243 Xphys=Data{3}.Xphys; 244 end 245 246 247 248 elseif NbView==4 % is there is 2 stereo folders, get global U and V and compute Zphys 249 250 251 %test if the seconde camera is the same for both folder 252 for i=3:4 253 indpt(i)=strfind(Param.InputTable{i,2},'.'); % indice of the "." is the folder name 1 254 indline(i)=strfind(Param.InputTable{i,2},''); % indice of the "" is the folder name1 255 camname{i}=Param.InputTable{i,2}(indline(i)+1:indpt(i)1);% extract the second camera name 256 end 257 258 if strcmp(camname{3},camname{4})==0 259 disp_uvmat('ERROR','The 2 stereo folders should have the same camera for the second position',checkrun) 260 return 261 end 262 263 264 265 [Data{3},tild,errormsg] = nc2struct([Param.InputTable{3,1},'/',Param.InputTable{3,2},'/',Param.InputTable{3,3},'_',int2str(first_img+indextemp1),'.nc']); 266 267 if exist('Data{3}.Civ3_FF','var') % if FF is present, remove wrong vector 268 temp=find(Data{3}.Civ3_FF==0); 269 Xmid3=Data{3}.Xmid(temp); 270 Ymid3=Data{3}.Ymid(temp); 271 U3=Data{3}.Uphys(temp); 272 V3=Data{3}.Vphys(temp); 273 else 274 Xmid3=Data{3}.Xmid; 275 Ymid3=Data{3}.Ymid; 276 U3=Data{3}.Uphys; 277 V3=Data{3}.Vphys; 278 end 279 %temporary gridd of merging the 2 stereos datas 280 [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)); 281 282 %1st folder : interpolate the first camera (Dalsa1) points on the second (common) camera 283 %(Dalsa 3) 284 x3Q=griddata(Xmid3+(U3)/2,Ymid3+(V3)/2,Xmid3(U3)/2,xq,yq); 285 y3Q=griddata(Xmid3+(U3)/2,Ymid3+(V3)/2,Ymid3(V3)/2,xq,yq); 286 287 288 289 [Data{4},tild,errormsg] = nc2struct([Param.InputTable{4,1},'/',Param.InputTable{4,2},'/',Param.InputTable{4,3},'_',int2str(first_img+indextemp1),'.nc']); 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 303 %2nd folder :interpolate the first camera (Dalsa2) points on the second (common) camera 304 %(Dalsa 3) 305 x4Q=griddata(Xmid4+(U4)/2,Ymid4+(V4)/2,Xmid4(U4)/2,xq,yq); 306 y4Q=griddata(Xmid4+(U4)/2,Ymid4+(V4)/2,Ymid4(V4)/2,xq,yq); 307 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 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 324 if NbView>2 325 ZItemp(:,:,idtemp)=griddata(Xphys,Yphys,Zphys,XI,YI); %interpolation on the choosen gridd 326 end 327 328 end 329 ZI=mean(ZItemp,3); %mean between two the two time step 330 Vtest=ZItemp(:,:,2)ZItemp(:,:,1); 306 if NbView>2 307 ZItemp(:,:,idtemp)=griddata(Xphys,Yphys,Zphys,XI,YI); %interpolation on the choosen gridd 308 CheckZ=1; 309 end 310 end 311 ZI=mean(ZItemp,3); %mean between the two time step 331 312 332 313 [Xa,Ya]=px_XYZ(XmlData{1}.GeometryCalib,[],XI,YI,ZI);% set of image coordinates on view a 333 314 [Xb,Yb]=px_XYZ(XmlData{2}.GeometryCalib,[],XI,YI,ZI);% set of image coordinates on view b 334 315 335 316 336 317 for iview=1:2 337 318 %% reading input file(s) … … 342 323 end 343 324 % get the time defined in the current file if not already defined from the xml file 344 if isfield(Data{iview},'Time')&& isequal(Data{iview}.Time,Data{1}.Time)325 if isfield(Data{iview},'Time')&& (Data{iview}.TimeData{1}.Time)<0.0001 345 326 Time=Data{iview}.Time; 346 327 else … … 355 336 end 356 337 end 357 %remove wrong vector 338 %remove wrong vector 358 339 if isfield(Data{1},'FF') 359 340 temp=find(Data{1}.FF==0); … … 373 354 [A]=get_coeff(XmlData{1}.GeometryCalib,Xa,Ya,XI,YI,ZI); %get coef A~ 374 355 375 %remove wrong vector 356 %remove wrong vector 376 357 if isfield(Data{2},'FF') 377 358 temp=find(Data{2}.FF==0); … … 396 377 S=ones(size(XI,1),size(XI,2),3); 397 378 D=ones(size(XI,1),size(XI,2),3,3); 398 379 399 380 S(:,:,1)=A(:,:,1,1).*Ua+A(:,:,2,1).*Va+B(:,:,1,1).*Ub+B(:,:,2,1).*Vb; 400 381 S(:,:,2)=A(:,:,1,2).*Ua+A(:,:,2,2).*Va+B(:,:,1,2).*Ub+B(:,:,2,2).*Vb; … … 416 397 W(indj,indi)=dxyz(3); 417 398 end 418 end 399 end 419 400 Error=zeros(size(XI,1),size(XI,2),4); 420 401 Error(:,:,1)=A(:,:,1,1).*U+A(:,:,1,2).*V+A(:,:,1,3).*WUa; … … 424 405 425 406 426 427 428 429 430 407 %% recording the merged field 431 408 if index==1% initiate the structure at first index 432 MergeData.ListGlobalAttribute={'Conventions','Time','Dt' };409 MergeData.ListGlobalAttribute={'Conventions','Time','Dt','CoordUnit'}; 433 410 MergeData.Conventions='uvmat'; 434 MergeData.Time=Time; 435 MergeData.Dt=Dt; 436 MergeData.ListVarName={'coord_x','coord_y','Z','U','V','W','Error'}; 411 if isfield (XmlData{1}.GeometryCalib,'CoordUnit') && isfield (XmlData{2}.GeometryCalib,'CoordUnit') && strcmp(XmlData{1}.GeometryCalib.CoordUnit, XmlData{2}.GeometryCalib.CoordUnit) 412 MergeData.CoordUnit=XmlData{1}.GeometryCalib.CoordUnit; 413 else 414 disp_uvmat('ERROR','inconsistent coord units in the two input velocity series',checkrun) 415 return 416 end 417 MergeData.ListVarName={'coord_x','coord_y','U','V','W','Error'}; 437 418 MergeData.VarDimName={'coord_x','coord_y',{'coord_y','coord_x'},{'coord_y','coord_x'}... 438 {'coord_y','coord_x'},{'coord_y','coord_x'},{'coord_y','coord_x'}}; 419 {'coord_y','coord_x'},{'coord_y','coord_x'}}; 420 MergeData.VarAttribute{6}.unit='pixel'; %error estimate expressed in pixel 421 if CheckZ 422 MergeData.ListVarName=[MergeData.ListVarName {'Z'}]; 423 MergeData.VarDimName=[MergeData.ListVarName {'coord_y','coord_x'}]; 424 MergeData.Z=ZI; 425 end 439 426 MergeData.coord_x=xI; 440 427 MergeData.coord_y=yI; 441 428 end 429 MergeData.Time=Time; 430 MergeData.Dt=Dt; 442 431 MergeData.U=U/Dt; 443 432 MergeData.V=V/Dt; 444 433 MergeData.W=W/Dt; 445 MergeData.Z=ZI;446 447 %mfx=(XmlData{1}.GeometryCalib.fx_fy(1)+XmlData{2}.GeometryCalib.fx_fy(1))/2;448 %mfy=(XmlData{1}.GeometryCalib.fx_fy(2)+XmlData{2}.GeometryCalib.fx_fy(2))/2;449 MergeData.Error=0. 5*sqrt(sum(Error.^2,3));434 435 436 mfx=(XmlData{1}.GeometryCalib.fx_fy(1)+XmlData{2}.GeometryCalib.fx_fy(1))/2; 437 mfy=(XmlData{1}.GeometryCalib.fx_fy(2)+XmlData{2}.GeometryCalib.fx_fy(2))/2; 438 MergeData.Error=0.25*(mfx+mfy)*sqrt(sum(Error.^2,3)); 450 439 errormsg=struct2nc(OutputFile,MergeData);%save result file 451 440 if isempty(errormsg)
Note: See TracChangeset
for help on using the changeset viewer.