Changeset 1158
- Timestamp:
- Jul 12, 2024, 8:46:38 AM (3 months ago)
- Location:
- trunk/src
- Files:
-
- 1 added
- 1 deleted
- 6 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/src/mouse_motion.m
r1143 r1158 165 165 end 166 166 %display the field values 167 if isfield (CellInfo{icell},'XName') 168 XName=CellInfo{icell}.XName; 169 text_displ_2=[XName '=' num2str(Field.(XName)(ivec),4)]; 170 end 171 if isfield (CellInfo{icell},'YName') 172 YName=CellInfo{icell}.YName; 173 text_displ_2=[text_displ_2 ',' YName '=' num2str(Field.(YName)(ivec),4)]; 174 end 167 175 for ivar=1:numel(CellInfo{icell}.VarIndex) 168 176 VarName=Field.ListVarName{CellInfo{icell}.VarIndex(ivar)}; 169 177 VarVal=Field.(VarName)(ivec); 170 178 var_text=[VarName '=' num2str(VarVal,4) ',']; 179 171 180 if isequal(ivar,CellInfo{icell}.CoordIndex(end))||isequal(ivar,CellInfo{icell}.CoordIndex(end-1))||isequal(ivar,CellInfo{icell}.CoordIndex(1)) 172 181 text_displ_1=[text_displ_1 var_text]; -
trunk/src/plot_field.m
r1152 r1158 1331 1331 n=size(xc); 1332 1332 xN=NaN*ones(size(xc)); 1333 matx=[xc(:) -uc(:)/2 xc(:)+uc(:)/2xN(:)]';1333 matx=[xc(:) xc(:)+uc(:) xN(:)]'; 1334 1334 matx=reshape(matx,1,3*n(2)); 1335 maty=[yc(:) -vc(:)/2 yc(:)+vc(:)/2xN(:)]';1335 maty=[yc(:) yc(:)+vc(:) xN(:)]'; 1336 1336 maty=reshape(maty,1,3*n(2)); 1337 1337 … … 1339 1339 arrowplus=rot*[uc;vc]; 1340 1340 arrowmoins=rot'*[uc;vc]; 1341 x1=xc+uc/2-arrowplus(1,:); 1342 x2=xc+uc/2; 1343 x3=xc+uc/2-arrowmoins(1,:); 1344 y1=yc+vc/2-arrowplus(2,:); 1345 y2=yc+vc/2; 1346 y3=yc+vc/2-arrowmoins(2,:); 1341 % x1=xc+uc/2-arrowplus(1,:); 1342 % x2=xc+uc/2; 1343 % x3=xc+uc/2-arrowmoins(1,:); 1344 % y1=yc+vc/2-arrowplus(2,:); 1345 % y2=yc+vc/2; 1346 % y3=yc+vc/2-arrowmoins(2,:); 1347 x1=xc+uc-arrowplus(1,:); 1348 x2=xc+uc; 1349 x3=xc+uc-arrowmoins(1,:); 1350 y1=yc+vc-arrowplus(2,:); 1351 y2=yc+vc; 1352 y3=yc+vc-arrowmoins(2,:); 1347 1353 matxar=[x1(:) x2(:) x3(:) xN(:)]'; 1348 1354 matxar=reshape(matxar,1,4*n(2)); -
trunk/src/proj_field.m
r1135 r1158 1360 1360 if ~(check3D && ivar==CellInfo{icell}.CoordIndex(1)) 1361 1361 ListVarName=[ListVarName VarName]; 1362 VarDimName=[VarDimName DimCell];1362 VarDimName=[VarDimName {DimCell}]; 1363 1363 nbvar=nbvar+1; 1364 1364 if isfield(FieldData,'VarAttribute') && length(FieldData.VarAttribute) >=ivar -
trunk/src/series/civ_3D.m
r1157 r1158 180 180 Data.ListGlobalAttribute=[ListGlobalAttribute Civ1_param]; 181 181 % set the list of variables 182 Data.ListVarName={'Coord_z','C oord_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 record183 Data.VarDimName={'npz', 'npy','npx',{'npz','npy','npx'},{'npz','npy','npx'},{'npz','npy','npx'},{'npz','npy','npx'},...182 Data.ListVarName={'Coord_z','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',{'npz','npy','npx'},{'npz','npy','npx'},{'npz','npy','npx'},{'npz','npy','npx'},... 184 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'; 185 Data.VarAttribute{1}.Role='coord_z'; 186 Data.VarAttribute{2}.Role='coord_x'; 187 ata.VarAttribute{3}.Role='coord_y'; 188 Data.VarAttribute{4}.Role='vector_x'; 189 Data.VarAttribute{5}.Role='vector_y'; 190 Data.VarAttribute{6}.Role='vector_z'; 191 Data.VarAttribute{7}.Role='ancillary'; 192 Data.VarAttribute{8}.Role='errorflag'; 192 193 Data.Coord_z=j1_series_Civ1; 193 194 % path for output nc file … … 297 298 par_civ1.ImageB=zeros(2*SearchRange_z+1,npy,npx); 298 299 299 z_index=1;%first vertical block centered at image index z_index= par_civ1.Dz300 for iz=1: par_civ1.Dz+SearchRange_z300 z_index=1;%first vertical block centered at image index z_index=SearchRange_z+1 301 for iz=1:2*SearchRange_z+1 301 302 j_image_index=j1_series_Civ1(iz,1)% j index of the current image 302 303 ImageName_A=fullfile_uvmat(RootPath_A,SubDir_A,RootFile_A,FileExt_A,NomType_A,i1_series_Civ1(1,ifield),[],j_image_index);% … … 304 305 ImageName_B=fullfile_uvmat(RootPath_B,SubDir_B,RootFile_B,FileExt_B,NomType_B,i2_series_Civ1(1,ifield),[],j_image_index); 305 306 B= read_image(ImageName_B,FileType_B); 306 par_civ1.ImageA(iz +par_civ1.Dz-1,:,:) = A;307 par_civ1.ImageB(iz +par_civ1.Dz-1,:,:) = B;307 par_civ1.ImageA(iz,:,:) = A; 308 par_civ1.ImageB(iz,:,:) = B; 308 309 end 309 310 % calculate velocity data at the first position z, corresponding to j_index=par_civ1.Dz … … 338 339 Data.Civ1_C(Data.Civ1_C==Inf)=0; 339 340 [npz,npy,npx]=size(Data.Civ1_X); 340 Data.Coord_z= 1:npz;341 Data.Coord_z=SearchRange_z+1:par_civ1.Dz:NbSlice-1; 341 342 Data.Coord_y=1:npy; 342 343 Data.Coord_x=1:npx; … … 788 789 789 790 %% prepare measurement grid 790 791 minix=floor(par_civ.Dx/2)-0.5; 792 maxix=minix+par_civ.Dx*floor((par_civ.ImageWidth-1)/par_civ.Dx); 793 miniy=floor(par_civ.Dy/2)-0.5;% first automatic grid point at half the mesh Dy 794 maxiy=miniy+par_civ.Dy*floor((par_civ.ImageHeight-1)/par_civ.Dy); 795 [GridX,GridY]=meshgrid(minix:par_civ.Dx:maxix,miniy:par_civ.Dy:maxiy); 791 nbinterv_x=floor((par_civ.ImageWidth-1)/par_civ.Dx); 792 gridlength_x=nbinterv_x*par_civ.Dx; 793 minix=ceil((par_civ.ImageWidth-gridlength_x)/2); 794 nbinterv_y=floor((par_civ.ImageHeight-1)/par_civ.Dy); 795 gridlength_y=nbinterv_y*par_civ.Dy; 796 miniy=ceil((par_civ.ImageHeight-gridlength_y)/2); 797 [GridX,GridY]=meshgrid(minix:par_civ.Dx:par_civ.ImageWidth-1,miniy:par_civ.Dy:par_civ.ImageHeight-1); 798 % minix:par_civ.Dx:par_civ.ImageWidth-1 799 % miniy:par_civ.Dy:par_civ.ImageHeight-1 800 % minix=floor(par_civ.Dx/2)-0.5; 801 % maxix=minix+par_civ.Dx*floor((par_civ.ImageWidth-1)/par_civ.Dx); 802 % miniy=floor(par_civ.Dy/2)-0.5;% first automatic grid point at half the mesh Dy 803 % maxiy=miniy+par_civ.Dy*floor((par_civ.ImageHeight-1)/par_civ.Dy); 804 % [GridX,GridY]=meshgrid(minix:par_civ.Dx:maxix,miniy:par_civ.Dy:maxiy); 796 805 par_civ.Grid(:,:,1)=GridX; 797 806 par_civ.Grid(:,:,2)=GridY;% increases with array index, … … 931 940 end 932 941 933 934 935 for kz=1:par_civ.SearchBoxSize(3) 942 npxycorr=par_civ.SearchBoxSize(1:2)-par_civ.CorrBoxSize(1:2)+1; 943 result_conv=zeros([par_civ.SearchBoxSize(3) npxycorr]);%initialise the conv product 944 max_xy=zeros(par_civ.SearchBoxSize(3),1);%initialise the max correlation vs z 945 xk=ones(par_civ.SearchBoxSize(3),1);%initialise the x displacement corresponding to the max corr 946 yk=ones(par_civ.SearchBoxSize(3),1);%initialise the y displacement corresponding to the max corr 947 948 for kz=kref:kref%par_civ.SearchBoxSize(3) 936 949 subima2=squeeze(image2_crop(kz,:,:)); 937 950 subima1=squeeze(image1_crop(kref,:,:)); 938 951 correl_xy=conv2(subima2,flip(flip(subima1,2),1),'valid'); 939 result_conv(kz,:,:)=correl_xy;952 result_conv(kz,:,:)= 255*correl_xy; 940 953 max_xy(kz)=max(max(correl_xy)); 941 [xk(kz),yk(kz)]=find(correl_xy==max_xy(kz),1);942 954 % max_xy(kz)=max(max(correl_xy)); 955 [yk(kz),xk(kz)]=find(correl_xy==max_xy(kz),1); 943 956 end 944 [corrmax,zmax]=max(max_xy); 945 946 x=xk(zmax); 947 y=yk(zmax); 948 result_conv=(result_conv/corrmax)*255; %normalize, peak=always 255 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,:,:)))); 957 [corrmax,dz]=max(max_xy); 958 dx=xk(dz); 959 dy=yk(dz); 960 961 % 962 % result_conv=255*result_conv/corrmax; %normalize, peak=always 255 963 % for kz=1:par_civ.SearchBoxSize(3) 964 % [dy,dx] = find(result_conv(kz,:,:)==255,1) 965 % if ~isempty(dy) 966 % dz=kz; 967 % break 968 % end 969 % end 970 subimage2_crop=squeeze(image2_crop(dz,dy:dy+2*iby2/mesh,dx:dx+2*ibx2/mesh));%subimage of image 2 corresponding to the optimum displacement of first image 971 sum_square=sum(sum(squeeze(image1_crop(dz,:,:).*image1_crop(dz,:,:)))); 951 972 sum_square=sum_square*sum(sum(subimage2_crop.*subimage2_crop));% product of variances of image 1 and 2 952 973 sum_square=sqrt(sum_square);% srt of the variance product to normalise correlation 953 if ~isempty(y) && ~isempty(x) 954 955 if par_civ.CorrSmooth==1 956 [vector,FF(ivec_y,ivec_x)] = SUBPIXGAUSS (squeeze(result_conv(zmax,:,:)),x,y);%TODO: improve by max optimisation along z 957 elseif par_civ.CorrSmooth==2 958 [vector,FF(ivec_y,ivec_x)] = SUBPIX2DGAUSS (squeeze(result_conv(zmax,:,:)),x,y); 959 else 960 [vector,FF(ivec_y,ivec_x)] = quadr_fit(squeeze(result_conv(zmax,:,:)),x,y); 961 end 962 utable(ivec_y,ivec_x)=vector(1)*mesh+shiftx(ivec_y,ivec_x); 963 vtable(ivec_y,ivec_x)=vector(2)*mesh+shifty(ivec_y,ivec_x); 964 xtable(ivec_y,ivec_x)=iref+utable(ivec_y,ivec_x)/2-0.5;% convec flow (velocity taken at the point middle from imgae 1 and 2) 965 ytable(ivec_y,ivec_x)=jref+vtable(ivec_y,ivec_x)/2-0.5;% and position of pixel 1=0.5 (convention for image coordinates=0 at the edge) 966 iref=round(xtable(ivec_y,ivec_x)+0.5);% nearest image index for the middle of the vector 967 jref=round(ytable(ivec_y,ivec_x)+0.5); 968 wtable(ivec_y,ivec_x)=zmax-kref; 969 % eliminate vectors located in the mask 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)) 971 utable(ivec_y,ivec_x)=0; 972 vtable(ivec_y,ivec_x)=0; 973 FF(ivec_y,ivec_x)=1; 974 end 975 ctable(ivec_y,ivec_x)=corrmax/sum_square;% correlation value 976 974 if ~isempty(dz)&& ~isempty(dy) && ~isempty(dx) 975 if par_civ.CorrSmooth==1 976 [vector,FF(ivec_y,ivec_x)] = SUBPIXGAUSS (result_conv,dx,dy,dz); 977 elseif par_civ.CorrSmooth==2 978 [vector,FF(ivec_y,ivec_x)] = SUBPIX2DGAUSS (result_conv,dx,dy,dz); 979 else 980 [vector,FF(ivec_y,ivec_x)] = quadr_fit(result_conv,dx,dy,dz); 981 end 982 utable(ivec_y,ivec_x)=vector(1)+shiftx(ivec_y,ivec_x); 983 vtable(ivec_y,ivec_x)=vector(2)+shifty(ivec_y,ivec_x); 984 wtable(ivec_y,ivec_x)=vector(3); 985 xtable(ivec_y,ivec_x)=iref+utable(ivec_y,ivec_x)/2-0.5;% convec flow (velocity taken at the point middle from imgae 1 and 2) 986 ytable(ivec_y,ivec_x)=jref+vtable(ivec_y,ivec_x)/2-0.5;% and position of pixel 1=0.5 (convention for image coordinates=0 at the edge) 987 iref=round(xtable(ivec_y,ivec_x)+0.5);% nearest image index for the middle of the vector 988 jref=round(ytable(ivec_y,ivec_x)+0.5); 989 990 % eliminate vectors located in the mask 991 if checkmask && (iref<1 || jref<1 ||iref>npx_ima || jref>npy_ima ||( par_civ.Mask(jref,iref)<200 && par_civ.Mask(jref,iref)>=100)) 992 utable(ivec_y,ivec_x)=0; 993 vtable(ivec_y,ivec_x)=0; 994 FF(ivec_y,ivec_x)=1; 995 end 996 ctable(ivec_y,ivec_x)=corrmax/sum_square;% correlation value 977 997 else 978 998 FF(ivec_y,ivec_x)=true; … … 993 1013 % x,y: position of the maximum correlation at integer values 994 1014 995 function [vector,FF] = SUBPIXGAUSS (result_conv,x,y )1015 function [vector,FF] = SUBPIXGAUSS (result_conv,x,y,z) 996 1016 %------------------------------------------------------------------------ 997 vector=[0 0]; %default 998 [np y,npx]=size(result_conv);1017 1018 [npz,npy,npx]=size(result_conv); 999 1019 result_conv(result_conv<1)=1; %set to 1 correlation values smaller than 1 (=0 by discretisation, to avoid divergence in the log) 1000 1020 %the following 8 lines are copyright (c) 1998, Uri Shavit, Roi Gurka, Alex Liberzon, Technion ??? Israel Institute of Technology 1001 1021 %http://urapiv.wordpress.com 1002 1022 FF=false; 1003 peaky = y;1004 if y < npy && y> 11005 f0 = log(result_conv(y,x));1006 f1 = log(result_conv(y-1,x));1007 f2 = log(result_conv(y+1,x));1008 peaky = peaky+ (f1-f2)/(2*f1-4*f0+2*f2);1009 else 1010 FF=true; % error flag for vector truncated by the limited search box in y1011 end1012 peak x=x;1013 if x < npx-1 && x > 1 1014 f0 = log(result_conv(y,x));1015 f1 = log(result_conv(y,x-1));1016 f2 = log(result_conv(y,x+1));1017 peakx = peakx+ (f1-f2)/(2*f1-4*f0+2*f2); 1018 else1019 FF=true; % warning flag for vector truncated by the limited search box in x 1020 end 1021 vector=[peakx-floor(npx/2)-1 peaky-floor(npy/2)-1];1023 peakz=z;peaky=y;peakx=x; 1024 if z < npz && z > 1 && y < npy && y > 1 && x < npx && x > 1 1025 f0 = log(result_conv(z,y,x)); 1026 f1 = log(result_conv(z-1,y,x)); 1027 f2 = log(result_conv(z+1,y,x)); 1028 peakz = peakz+ (f1-f2)/(2*f1-4*f0+2*f2); 1029 1030 f1 = log(result_conv(z,y-1,x)); 1031 f2 = log(result_conv(z,y+1,x)); 1032 peaky = peaky+ (f1-f2)/(2*f1-4*f0+2*f2); 1033 1034 f1 = log(result_conv(z,y,x-1)); 1035 f2 = log(result_conv(z,y,x+1)); 1036 peakx = peakx+ (f1-f2)/(2*f1-4*f0+2*f2); 1037 else 1038 FF=true; 1039 end 1040 1041 vector=[peakx-floor(npx/2)-1 peaky-floor(npy/2)-1 peakz-floor(npz/2)-1]; 1022 1042 1023 1043 %------------------------------------------------------------------------ -
trunk/src/series/civ_series.m
r1157 r1158 962 962 % else par_civ.Grid is already an array, no action here 963 963 else% automatic grid in x, y image coordinates 964 minix=floor(par_civ.Dx/2)-0.5; 965 maxix=minix+par_civ.Dx*floor((par_civ.ImageWidth-1)/par_civ.Dx); 966 miniy=floor(par_civ.Dy/2)-0.5;% first automatic grid point at half the mesh Dy 967 maxiy=minix+par_civ.Dy*floor((par_civ.ImageHeight-1)/par_civ.Dy); 968 [GridX,GridY]=meshgrid(minix:par_civ.Dx:maxix,miniy:par_civ.Dy:maxiy); 964 nbinterv_x=floor((par_civ.ImageWidth-1)/par_civ.Dx); 965 gridlength_x=nbinterv_x*par_civ.Dx; 966 minix=ceil((par_civ.ImageWidth-gridlength_x)/2); 967 nbinterv_y=floor((par_civ.ImageHeight-1)/par_civ.Dy); 968 gridlength_y=nbinterv_y*par_civ.Dy; 969 miniy=ceil((par_civ.ImageHeight-gridlength_y)/2); 970 [GridX,GridY]=meshgrid(minix:par_civ.Dx:par_civ.ImageWidth-1,miniy:par_civ.Dy:par_civ.ImageHeight-1); 969 971 par_civ.Grid(:,1)=reshape(GridX,[],1); 970 972 par_civ.Grid(:,2)=reshape(GridY,[],1);% increases with array index 973 % 974 % 975 % minix=floor(par_civ.Dx/2)-0.5; 976 % maxix=minix+par_civ.Dx*floor((par_civ.ImageWidth-1)/par_civ.Dx); 977 % miniy=floor(par_civ.Dy/2)-0.5;% first automatic grid point at half the mesh Dy 978 % maxiy=minix+par_civ.Dy*floor((par_civ.ImageHeight-1)/par_civ.Dy); 979 % [GridX,GridY]=meshgrid(minix:par_civ.Dx:maxix,miniy:par_civ.Dy:maxiy); 980 % par_civ.Grid(:,1)=reshape(GridX,[],1); 981 % par_civ.Grid(:,2)=reshape(GridY,[],1);% increases with array index 971 982 end 972 983 nbvec=size(par_civ.Grid,1); … … 1007 1018 par_civ.ImageA=sum(double(par_civ.ImageA),3);%sum over rgb component for color images 1008 1019 par_civ.ImageB=sum(double(par_civ.ImageB),3); 1009 [npy_ima 1020 [npy_ima,npx_ima]=size(par_civ.ImageA); 1010 1021 if ~isequal(size(par_civ.ImageB),[npy_ima npx_ima]) 1011 1022 errormsg='image pair with unequal size'; … … 1116 1127 sum_square=sum(sum(image1_crop.*image1_crop)); 1117 1128 %reference: Oliver Pust, PIV: Direct Cross-Correlation 1118 result_conv= conv2(image2_crop,flip dim(flipdim(image1_crop,2),1),'valid');1129 result_conv= conv2(image2_crop,flip(flip(image1_crop,2),1),'valid'); 1119 1130 corrmax= max(max(result_conv)); 1120 1131 result_conv=(result_conv/corrmax)*255; %normalize, peak=always 255 … … 1176 1187 %the following 8 lines are copyright (c) 1998, Uri Shavit, Roi Gurka, Alex Liberzon, Technion ??? Israel Institute of Technology 1177 1188 %http://urapiv.wordpress.com 1178 peaky = y; 1179 if y < npy && y > 1 1189 peaky = y;peakx=x; 1190 if y < npy && y > 1 && x < npx-1 && x > 1 1180 1191 f0 = log(result_conv(y,x)); 1181 1192 f1 = log(result_conv(y-1,x)); 1182 1193 f2 = log(result_conv(y+1,x)); 1183 1194 peaky = peaky+ (f1-f2)/(2*f1-4*f0+2*f2); 1184 else1185 F=1; % warning flag for vector truncated by the limited search box1186 end1187 peakx=x;1188 if x < npx-1 && x > 11189 f0 = log(result_conv(y,x));1190 1195 f1 = log(result_conv(y,x-1)); 1191 1196 f2 = log(result_conv(y,x+1)); … … 1194 1199 F=1; % warning flag for vector truncated by the limited search box 1195 1200 end 1201 1196 1202 vector=[peakx-floor(npx/2)-1 peaky-floor(npy/2)-1]; 1197 1203 -
trunk/src/series/filter_time.m
r1157 r1158 195 195 % DataOut.Wfilter=squeeze(mean(Wblock(end-sumindex:end,:,:),1,'omitnan')); 196 196 %updating output the netcdf file 197 [i1,i2,j1,j2] = get_file_index(i1_series{1}(index),[],PairString {1});197 [i1,i2,j1,j2] = get_file_index(i1_series{1}(index),[],PairString); 198 198 OutputFile=fullfile_uvmat(OutputPath,OutputDir,Param.InputTable{1,3},'.nc',NomType{1},i1,i2,j1,j2); 199 199 errormsg=struct2nc(OutputFile, DataOut);
Note: See TracChangeset
for help on using the changeset viewer.