Changeset 1157 for trunk/src/series/civ_3D.m
 Timestamp:
 Jul 11, 2024, 4:13:03 PM (3 months ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

trunk/src/series/civ_3D.m
r1156 r1157 130 130 end 131 131 132 [ tild,i1_series,i2_series,j1_series,j2_series]=get_file_series(Param);132 [~,i1_series,~,j1_series,~]=get_file_series(Param); 133 133 iview_A=0;% series index (iview) for the first image series 134 134 iview_B=0;% series index (iview) for the second image series (only non zero for option 'shift' comparing two image series ) … … 156 156 PairCiv1=Param.ActionInput.PairIndices.ListPairCiv1; 157 157 158 [i1_series_Civ1,i2_series_Civ1,j1_series_Civ1,j2_series_Civ1, check_bounds,NomTypeNc]=...158 [i1_series_Civ1,i2_series_Civ1,j1_series_Civ1,j2_series_Civ1,~,NomTypeNc]=... 159 159 find_pair_indices(PairCiv1,i1_series{1},j1_series{1},MinIndex_i,MaxIndex_i,MinIndex_j,MaxIndex_j); 160 160 … … 168 168 %% prepare output Data 169 169 ListGlobalAttribute={'Conventions','Program','CivStage'}; 170 Data.Conventions='uvmat/civdata';% states the conventions used for the description of field variables and attributes 171 Data.Program='civ_series'; 170 Data.Conventions='uvmat/civdata_3D';% states the conventions used for the description of field variables and attributes 171 Data.Program='civ_3D'; 172 if isfield(Param,'UvmatRevision') 173 Data.Program=[Data.Program ', uvmat r' Param.UvmatRevision]; 174 end 172 175 Data.CivStage=0;%default 173 176 list_param=(fieldnames(Param.ActionInput.Civ1))'; 174 175 176 177 list_param(strcmp('TestCiv1',list_param))=[];% remove the parameter TestCiv1 from the list 178 Civ1_param=regexprep(list_param,'^.+','Civ1_$0');% insert 'Civ1_' before each string in list_param 179 Civ1_param=[{'Civ1_ImageA','Civ1_ImageB','Civ1_Time','Civ1_Dt'} Civ1_param]; %insert the names of the two input images 177 180 Data.ListGlobalAttribute=[ListGlobalAttribute Civ1_param]; 178 % set the list of variables 179 Data.ListVarName={'Civ1_X','Civ1_Y','Civ1_U','Civ1_V','Civ1_W','Civ1_C','Civ1_FF'};% cell array containing the names of the fields to record 180 Data.VarDimName={'nb_vec_1','nb_vec_1','nb_vec_1','nb_vec_1','nb_vec_1','nb_vec_1','nb_vec_1'}; 181 Data.VarAttribute{1}.Role='coord_x'; 182 Data.VarAttribute{2}.Role='coord_y'; 183 Data.VarAttribute{3}.Role='vector_x'; 184 Data.VarAttribute{4}.Role='vector_y'; 185 Data.VarAttribute{5}.Role='vector_z'; 186 Data.VarAttribute{6}.Role='ancillary'; 187 Data.VarAttribute{7}.Role='errorflag'; 181 % set the list of variables 182 Data.ListVarName={'Coord_z','Coord_y','Coord_x','Civ1_X','Civ1_Y','Civ1_U','Civ1_V','Civ1_W','Civ1_C','Civ1_FF'};% cell array containing the names of the fields to record 183 Data.VarDimName={'npz','npy','npx',{'npz','npy','npx'},{'npz','npy','npx'},{'npz','npy','npx'},{'npz','npy','npx'},... 184 {'npz','npy','npx'},{'npz','npy','npx'},{'npz','npy','npx'}}; 185 Data.VarAttribute{1}.Role='coord_x'; 186 Data.VarAttribute{2}.Role='coord_y'; 187 Data.VarAttribute{3}.Role='vector_x'; 188 Data.VarAttribute{4}.Role='vector_y'; 189 Data.VarAttribute{5}.Role='vector_z'; 190 Data.VarAttribute{6}.Role='ancillary'; 191 Data.VarAttribute{7}.Role='errorflag'; 192 Data.Coord_z=j1_series_Civ1; 193 % path for output nc file 194 OutputPath=fullfile(Param.OutputPath,num2str(Param.Experiment),num2str(Param.Device),[Param.OutputSubDir Param.OutputDirExt]); 188 195 189 196 %% get timing from the ImaDoc file or input video … … 227 234 par_civ1.ImageA=zeros(2*SearchRange_z+1,npy,npx); 228 235 par_civ1.ImageB=zeros(2*SearchRange_z+1,npy,npx); 236 237 229 238 230 239 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% … … 268 277 269 278 Data.CivStage=1; 270 % output nc file271 ncfile_out=fullfile_uvmat(OutputPath, Param.InputTable{1,2},Param.InputTable{1,3},Param.InputTable{1,5},...272 NomTypeNc,i1_series_Civ1(ifield),i2_series_Civ1(ifield),j1_series_Civ1(ifield),j2_series_Civ1(ifield));279 280 ncfile_out=fullfile_uvmat(OutputPath,'',Param.InputTable{1,3},'.nc',... 281 '_11',i1_series_Civ1(ifield),i2_series_Civ1(ifield)); %output file name 273 282 274 283 if ~CheckOverwrite && exist(ncfile_out,'file') … … 284 293 disp('civ1 started') 285 294 286 287 % read input images (except in mode Test where it is introduced directly in Param.ActionInput.Civ1.ImageNameA and B) 288 289 par_civ1.ImageA=zeros(2*SearchRange_z+1,npy,npx); 295 % read input images by vertical blocks with nbre of images equal to 2*SearchRange+1, 296 par_civ1.ImageA=zeros(2*SearchRange_z+1,npy,npx);%image block initiation 290 297 par_civ1.ImageB=zeros(2*SearchRange_z+1,npy,npx); 291 %first vertical block centered at image index islice=par_civ1.Dz 292 islice=par_civ1.Dz; 293 for iz=1:par_civ1.Dz+SearchRange_z 294 ImageName_A=fullfile_uvmat(RootPath_A,SubDir_A,RootFile_A,FileExt_A,NomType_A,i1_series_Civ1(1,ifield),[],j1_series_Civ1(iz,1));% 298 299 z_index=1;%first vertical block centered at image index z_index=par_civ1.Dz 300 for iz=1:par_civ1.Dz+SearchRange_z 301 j_image_index=j1_series_Civ1(iz,1)% j index of the current image 302 ImageName_A=fullfile_uvmat(RootPath_A,SubDir_A,RootFile_A,FileExt_A,NomType_A,i1_series_Civ1(1,ifield),[],j_image_index);% 295 303 A= read_image(ImageName_A,FileType_A); 296 ImageName_B=fullfile_uvmat(RootPath_B,SubDir_B,RootFile_B,FileExt_B,NomType_B,i2_series_Civ1(1,ifield),[],j 1_series_Civ1(iz,1));304 ImageName_B=fullfile_uvmat(RootPath_B,SubDir_B,RootFile_B,FileExt_B,NomType_B,i2_series_Civ1(1,ifield),[],j_image_index); 297 305 B= read_image(ImageName_B,FileType_B); 298 306 par_civ1.ImageA(iz+par_civ1.Dz1,:,:) = A; 299 307 par_civ1.ImageB(iz+par_civ1.Dz1,:,:) = B; 300 end 301 % caluclate velocity data (y and v in indices, reverse to y component) 302 [Data.Civ1_X(islice,:,:),Data.Civ1_Y(islice,:,:), utable, vtable,wtable, ctable, FF, result_conv, errormsg] = civ3D (par_civ1); 308 end 309 % calculate velocity data at the first position z, corresponding to j_index=par_civ1.Dz 310 [Data.Civ1_X(1,:,:),Data.Civ1_Y(1,:,:), Data.Civ1_U(1,:,:), Data.Civ1_V(1,:,:),Data.Civ1_W(1,:,:),... 311 Data.Civ1_C(1,:,:), Data.Civ1_FF(1,:,:), result_conv, errormsg] = civ3D (par_civ1); 312 if ~isempty(errormsg) 313 disp_uvmat('ERROR',errormsg,checkrun) 314 return 315 end 316 for z_index=2:floor((NbSlice1)/par_civ1.Dz)% loop on slices 317 par_civ1.ImageA=circshift(par_civ1.ImageA,par_civ1.Dz);%shift the indices in the image block upward by par_civ1.Dz 318 par_civ1.ImageB=circshift(par_civ1.ImageA,par_civ1.Dz); 319 for iz=1:par_civ1.Dz %read the new images at the end of the image block 320 j_image_index=j1_series_Civ1(z_index*par_civ1.Dz+SearchRange_zpar_civ1.Dz+iz,1) 321 ImageName_A=fullfile_uvmat(RootPath_A,SubDir_A,RootFile_A,FileExt_A,NomType_A,i1_series_Civ1(1,ifield),[],j_image_index);% 322 A= read_image(ImageName_A,FileType_A); 323 ImageName_B=fullfile_uvmat(RootPath_B,SubDir_B,RootFile_B,FileExt_B,NomType_B,i2_series_Civ1(1,ifield),[],j_image_index); 324 B= read_image(ImageName_B,FileType_B); 325 par_civ1.ImageA(2*SearchRange_z+1par_civ1.Dz+iz,:,:) = A; 326 par_civ1.ImageB(2*SearchRange_z+1par_civ1.Dz+iz,:,:) = B; 327 end 328 % caluclate velocity data (y and v in indices, reverse to y component) 329 [Data.Civ1_X(z_index,:,:),Data.Civ1_Y(z_index,:,:), Data.Civ1_U(z_index,:,:), Data.Civ1_V(z_index,:,:),Data.Civ1_W(z_index,:,:),... 330 Data.Civ1_C(z_index,:,:), Data.Civ1_FF(z_index,:,:), result_conv, errormsg] = civ3D (par_civ1); 331 303 332 if ~isempty(errormsg) 304 333 disp_uvmat('ERROR',errormsg,checkrun) 305 334 return 306 335 end 307 for islice=2*par_civ1.Dz:NbSlice% loop on slices for the first image in the pair 308 par_civ1.ImageA=circshift(par_civ1.ImageA,par_civ1.Dz);%shift the indces in the block upward by par_civ1.Dz 309 par_civ1.ImageB=circshift(par_civ1.ImageA,par_civ1.Dz); 310 for iz=1:par_civ1.Dz 311 ImageName_A=fullfile_uvmat(RootPath_A,SubDir_A,RootFile_A,FileExt_A,NomType_A,i1_series_Civ1(1,ifield),[],j1_series_Civ1(islice+SearchRange_zpar_civ1.Dz+iz,1));% 312 A= read_image(ImageName_A,FileType_A); 313 ImageName_B=fullfile_uvmat(RootPath_B,SubDir_B,RootFile_B,FileExt_B,NomType_B,i2_series_Civ1(1,ifield),[],j1_series_Civ1(islice+SearchRange_zpar_civ1.Dz+iz,1)); 314 B= read_image(ImageName_B,FileType_B); 315 par_civ1.ImageA(2*SearchRange_z+1par_civ1.Dz+iz,:,:) = A; 316 par_civ1.ImageB(2*SearchRange_z+1par_civ1.Dz+iz,:,:) = B; 317 end 318 % caluclate velocity data (y and v in indices, reverse to y component) 319 [Data.Civ1_X(islice,:,:),Data.Civ1_Y(islice,:,:), utable, vtable,wtable, ctable, FF, result_conv, errormsg] = civ3D (par_civ1); 320 if ~isempty(errormsg) 321 disp_uvmat('ERROR',errormsg,checkrun) 322 return 323 end 324 325 end 336 end 337 Data.Civ1_C=double(Data.Civ1_C); 338 Data.Civ1_C(Data.Civ1_C==Inf)=0; 339 [npz,npy,npx]=size(Data.Civ1_X); 340 Data.Coord_z=1:npz; 341 Data.Coord_y=1:npy; 342 Data.Coord_x=1:npx; 326 343 % case of mask TO ADAPT 327 344 if par_civ1.CheckMask&&~isempty(par_civ1.Mask) … … 668 685 end 669 686 Data.ListGlobalAttribute=[Data.ListGlobalAttribute Fix2_param]; 670 Data.Civ2_FF=d ouble(detect_false(Param.ActionInput.Fix2,Data.Civ2_C,Data.Civ2_U,Data.Civ2_V,Data.Civ2_FF));687 Data.Civ2_FF=detect_false(Param.ActionInput.Fix2,Data.Civ2_C,Data.Civ2_U,Data.Civ2_V,Data.Civ2_FF); 671 688 Data.CivStage=Data.CivStage+1; 672 689 end … … 728 745 time_total=toc(tstart); 729 746 disp(['ellapsed time ' num2str(time_total/60,2) ' minutes']) 730 disp(['time image reading ' num2str(time_input,2) ' s'])731 747 disp(['time civ1 ' num2str(time_civ1,2) ' s']) 732 748 disp(['time patch1 ' num2str(time_patch1,2) ' s']) 733 749 disp(['time civ2 ' num2str(time_civ2,2) ' s']) 734 750 disp(['time patch2 ' num2str(time_patch2,2) ' s']) 735 disp(['time other ' num2str((time_totaltime_inputtime_civ1time_patch1time_civ2time_patch2),2) ' s']) 736 % end 751 737 752 end 738 753 … … 860 875 for ivec_x=1:nbvec_x 861 876 for ivec_y=1:nbvec_y 862 ivec_y 863 iref=round(par_civ.Grid(ivec_y,ivec_x,1)+0.5)% xindex on the image A for the middle of the correlation box 864 jref=round(par_civ.ImageHeightpar_civ.Grid(ivec_y,ivec_x,2)+0.5)% j index for the middle of the correlation box in the image A 877 iref=round(par_civ.Grid(ivec_y,ivec_x,1)+0.5);% xindex on the image A for the middle of the correlation box 878 jref=round(par_civ.ImageHeightpar_civ.Grid(ivec_y,ivec_x,2)+0.5);% j index for the middle of the correlation box in the image A 865 879 subrange1_x=irefibx2:iref+ibx2;% x indices defining the first subimage 866 880 subrange1_y=jrefiby2:jref+iby2;% y indices defining the first subimage … … 918 932 919 933 920 %reference: Oliver Pust, PIV: Direct CrossCorrelation934 921 935 for kz=1:par_civ.SearchBoxSize(3) 922 936 subima2=squeeze(image2_crop(kz,:,:)); … … 928 942 929 943 end 930 [corrmax,z ]=max(max_xy);944 [corrmax,zmax]=max(max_xy); 931 945 932 x=xk(z );933 y=yk(z );946 x=xk(zmax); 947 y=yk(zmax); 934 948 result_conv=(result_conv/corrmax)*255; %normalize, peak=always 255 935 subimage2_crop=squeeze(image2_crop(z ,y:y+2*iby2/mesh,x:x+2*ibx2/mesh));%subimage of image 2 corresponding to the optimum displacement of first image936 sum_square=sum(sum(squeeze(image1_crop(z ,:,:).*image1_crop(z,:,:))));949 subimage2_crop=squeeze(image2_crop(zmax,y:y+2*iby2/mesh,x:x+2*ibx2/mesh));%subimage of image 2 corresponding to the optimum displacement of first image 950 sum_square=sum(sum(squeeze(image1_crop(zmax,:,:).*image1_crop(zmax,:,:)))); 937 951 sum_square=sum_square*sum(sum(subimage2_crop.*subimage2_crop));% product of variances of image 1 and 2 938 952 sum_square=sqrt(sum_square);% srt of the variance product to normalise correlation … … 940 954 941 955 if par_civ.CorrSmooth==1 942 [vector,FF(ivec_y,ivec_x)] = SUBPIXGAUSS ( result_conv(z,:,:),x,y);%TODO: improve by max optimisation along z956 [vector,FF(ivec_y,ivec_x)] = SUBPIXGAUSS (squeeze(result_conv(zmax,:,:)),x,y);%TODO: improve by max optimisation along z 943 957 elseif par_civ.CorrSmooth==2 944 [vector,FF(ivec_y,ivec_x)] = SUBPIX2DGAUSS ( result_conv(z,:,:),x,y);958 [vector,FF(ivec_y,ivec_x)] = SUBPIX2DGAUSS (squeeze(result_conv(zmax,:,:)),x,y); 945 959 else 946 [vector,FF(ivec_y,ivec_x)] = quadr_fit( result_conv(z,:,:),x,y);960 [vector,FF(ivec_y,ivec_x)] = quadr_fit(squeeze(result_conv(zmax,:,:)),x,y); 947 961 end 948 962 utable(ivec_y,ivec_x)=vector(1)*mesh+shiftx(ivec_y,ivec_x); … … 952 966 iref=round(xtable(ivec_y,ivec_x)+0.5);% nearest image index for the middle of the vector 953 967 jref=round(ytable(ivec_y,ivec_x)+0.5); 954 wtable(ivec_y,ivec_x)=z kref;968 wtable(ivec_y,ivec_x)=zmaxkref; 955 969 % eliminate vectors located in the mask 956 970 if checkmask && (iref<1  jref<1 iref>npx_ima  jref>npy_ima ( par_civ.Mask(jref,iref)<200 && par_civ.Mask(jref,iref)>=100)) … … 962 976 963 977 else 964 FF(ivec_y,ivec_x)= 1;978 FF(ivec_y,ivec_x)=true; 965 979 end 966 980 end … … 975 989 % OUPUT: 976 990 % vector = optimum displacement vector with subpixel correction 977 % F =flag: =0 OK 978 % =2 , warning: max too close to the edge of the search box (1 pixel margin) 991 % FF =flag: ='true' max too close to the edge of the search box (1 pixel margin) 979 992 % INPUT: 980 993 % x,y: position of the maximum correlation at integer values 981 994 982 function [vector,F ] = SUBPIXGAUSS (result_conv,x,y)995 function [vector,FF] = SUBPIXGAUSS (result_conv,x,y) 983 996 % 984 % vector=[0 0]; %default 985 F=0; 997 vector=[0 0]; %default 986 998 [npy,npx]=size(result_conv); 987 999 result_conv(result_conv<1)=1; %set to 1 correlation values smaller than 1 (=0 by discretisation, to avoid divergence in the log) 988 1000 %the following 8 lines are copyright (c) 1998, Uri Shavit, Roi Gurka, Alex Liberzon, Technion ??? Israel Institute of Technology 989 1001 %http://urapiv.wordpress.com 990 peaky = y; 991 if y < npy && y > 1 992 f0 = log(result_conv(y,x)); 993 f1 = log(result_conv(y1,x)); 994 f2 = log(result_conv(y+1,x)); 995 peaky = peaky+ (f1f2)/(2*f14*f0+2*f2); 996 else 997 F=1; % warning flag for vector truncated by the limited search box 998 end 999 peakx=x; 1000 if x < npx1 && x > 1 1001 f0 = log(result_conv(y,x)); 1002 f1 = log(result_conv(y,x1)); 1003 f2 = log(result_conv(y,x+1)); 1004 peakx = peakx+ (f1f2)/(2*f14*f0+2*f2); 1005 else 1006 F=1; % warning flag for vector truncated by the limited search box 1007 end 1008 vector=[peakxfloor(npx/2)1 peakyfloor(npy/2)1]; 1002 FF=false; 1003 peaky = y; 1004 if y < npy && y > 1 1005 f0 = log(result_conv(y,x)); 1006 f1 = log(result_conv(y1,x)); 1007 f2 = log(result_conv(y+1,x)); 1008 peaky = peaky+ (f1f2)/(2*f14*f0+2*f2); 1009 else 1010 FF=true; % error flag for vector truncated by the limited search box in y 1011 end 1012 peakx=x; 1013 if x < npx1 && x > 1 1014 f0 = log(result_conv(y,x)); 1015 f1 = log(result_conv(y,x1)); 1016 f2 = log(result_conv(y,x+1)); 1017 peakx = peakx+ (f1f2)/(2*f14*f0+2*f2); 1018 else 1019 FF=true; % warning flag for vector truncated by the limited search box in x 1020 end 1021 vector=[peakxfloor(npx/2)1 peakyfloor(npy/2)1]; 1009 1022 1010 1023 %
Note: See TracChangeset
for help on using the changeset viewer.