Changeset 876 for trunk/src/series/civ_series.m
 Timestamp:
 Feb 22, 2015, 11:07:07 PM (9 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

trunk/src/series/civ_series.m
r873 r876 691 691 end 692 692 end 693 % ibx2=ceil(par_civ2.CorrBoxSize(1)/2);694 % iby2=ceil(par_civ2.CorrBoxSize(2)/2);695 % par_civ2.SearchBoxSize(1)=2*ibx2+9;% search ara +4 pixels around the guess696 %par_civ2.SearchBoxSize(2)=2*iby2+9;697 693 if CheckInputFile % else Dt given by par_civ2 698 694 if strcmp(Param.ActionInput.ListCompareMode,'displacement') … … 707 703 par_civ2.Grid=[par_civ2.Grid(nbval>=1,1)par_civ2.SearchBoxShift(:,1)/2 par_civ2.Grid(nbval>=1,2)par_civ2.SearchBoxShift(:,2)/2]; 708 704 if par_civ2.CheckDeformation 709 par_civ2.DUDX=DUDX ./nbval;710 par_civ2.DUDY=DUDY ./nbval;711 par_civ2.DVDX=DVDX ./nbval;712 par_civ2.DVDY=DVDY ./nbval;705 par_civ2.DUDX=DUDX(nbval>=1)./nbval(nbval>=1); 706 par_civ2.DUDY=DUDY(nbval>=1)./nbval(nbval>=1); 707 par_civ2.DVDX=DVDX(nbval>=1)./nbval(nbval>=1); 708 par_civ2.DVDY=DVDY(nbval>=1)./nbval(nbval>=1); 713 709 end 714 710 … … 910 906 911 907 %% prepare correlation and search boxes 912 ibx2= ceil(par_civ.CorrBoxSize(1)/2);913 iby2= ceil(par_civ.CorrBoxSize(2)/2);914 isx2= ceil(par_civ.SearchBoxSize(1)/2);915 isy2= ceil(par_civ.SearchBoxSize(2)/2);908 ibx2=floor(par_civ.CorrBoxSize(1)/2); 909 iby2=floor(par_civ.CorrBoxSize(2)/2); 910 isx2=floor(par_civ.SearchBoxSize(1)/2); 911 isy2=floor(par_civ.SearchBoxSize(2)/2); 916 912 shiftx=round(par_civ.SearchBoxShift(:,1));%use the input shift estimate, rounded to the next integer value 917 913 shifty=round(par_civ.SearchBoxShift(:,2));% sign minus because image j index increases when y decreases … … 968 964 end 969 965 check_undefined=(par_civ.Mask<200 & par_civ.Mask>=20 ); 970 par_civ.ImageA(check_undefined)=0;% put image A to zero (i.e. the min image value) in the undefined area971 par_civ.ImageB(check_undefined)=0;% put image B to zero (i.e. the min image value) in the undefined area966 % par_civ.ImageA(check_undefined)=0;% put image A to zero (i.e. the min image value) in the undefined area 967 % par_civ.ImageB(check_undefined)=0;% put image B to zero (i.e. the min image value) in the undefined area 972 968 end 973 969 … … 993 989 image2_crop=MinA*ones(numel(subrange2_y),numel(subrange2_x));% default value=min of image A 994 990 mask1_crop=ones(numel(subrange1_y),numel(subrange1_x));% default value=1 for mask 995 mask2_crop=ones(numel(subrange2_y),numel(subrange2_x));% default value= minfor mask991 mask2_crop=ones(numel(subrange2_y),numel(subrange2_x));% default value=1 for mask 996 992 check1_x=subrange1_x>=1 & subrange1_x<=par_civ.ImageWidth;% check which points in the subimage 1 are contained in the initial image 1 997 993 check1_y=subrange1_y>=1 & subrange1_y<=par_civ.ImageHeight; … … 1006 1002 F(ivec)=3; % 1007 1003 else 1004 image1_crop=image1_crop.*~mask1_crop;% put to zero the masked pixels (mask1_crop='true'=1) 1005 image2_crop=image2_crop.*~mask2_crop; 1008 1006 image1_mean=mean(mean(image1_crop))/(1sizemask); 1009 1007 image2_mean=mean(mean(image2_crop))/(1sizemask); … … 1026 1024 XIant=XIpar_civ.DUDX(ivec)*XIpar_civ.DUDY(ivec)*YI+ceil(size(image1_crop,2)/2); 1027 1025 YIant=YIpar_civ.DVDX(ivec)*XIpar_civ.DVDY(ivec)*YI+ceil(size(image1_crop,1)/2); 1026 image1_crop_1=image1_crop; 1028 1027 image1_crop=interp2(image1_crop,XIant,YIant); 1028 image1_crop_2=image1_crop; 1029 1029 image1_crop(isnan(image1_crop))=0; 1030 image1_crop_3=image1_crop; 1030 1031 xi=(1:mesh:size(image2_crop,2)); 1031 1032 yi=(1:mesh:size(image2_crop,1))'; 1032 1033 image2_crop=interp2(image2_crop,xi,yi); 1033 1034 image2_crop(isnan(image2_crop))=0; 1035 par_civ.CorrSmooth=3; 1034 1036 %% 1035 1037 end … … 1041 1043 %Find the correlation max, at 255 1042 1044 [y,x] = find(result_conv==255,1); 1045 subimage2_crop=image2_crop(y:y+2*iby2/mesh,x:x+2*ibx2/mesh);%subimage of image 2 corresponding to the optimum displacement of first image 1046 sum_square=sum_square*sum(sum(subimage2_crop.*subimage2_crop));% product of variances of image 1 and 2 1047 sum_square=sqrt(sum_square);% srt of the variance product to normalise correlation 1043 1048 if ~isempty(y) && ~isempty(x) 1044 1049 try … … 1047 1052 elseif par_civ.CorrSmooth==2 1048 1053 [vector,F(ivec)] = SUBPIX2DGAUSS (result_conv,x,y); 1054 else 1055 [vector,F(ivec)] = quadr_fit(result_conv,x,y); 1049 1056 end 1050 1057 utable(ivec)=vector(1)*mesh+shiftx(ivec); … … 1052 1059 xtable(ivec)=iref+utable(ivec)/20.5;% convec flow (velocity taken at the point middle from imgae 1 and 2) 1053 1060 ytable(ivec)=jref+vtable(ivec)/20.5;% and position of pixel 1=0.5 (convention for image coordinates=0 at the edge) 1054 iref=round(xtable(ivec));% image index for the middle of the vector1061 iref=round(xtable(ivec));% nearest image index for the middle of the vector 1055 1062 jref=round(ytable(ivec)); 1063 % eliminate vectors located in the mask 1056 1064 if checkmask && par_civ.Mask(jref,iref)<200 && par_civ.Mask(jref,iref)>=100 1057 1065 utable(ivec)=0; … … 1150 1158 end 1151 1159 vector=[peakxfloor(npx/2)1 peakyfloor(npy/2)1]; 1160 1161 % 1162 %  Find the maximum of the correlation function after quadratic interpolation 1163 function [vector,F] = quadr_fit(result_conv,x,y) 1164 [npy,npx]=size(result_conv); 1165 if x<4  y<4  npxx<4 npyy <4 1166 F=2; 1167 vector=[x y]; 1168 else 1169 F=0; 1170 x_ind=x4:x+4; 1171 y_ind=y4:y+4; 1172 x_vec=0.25*(x_indx); 1173 y_vec=0.25*(y_indy); 1174 [X,Y]=meshgrid(x_vec,y_vec); 1175 coord=[reshape(X,[],1) reshape(Y,[],1)]; 1176 result_conv=reshape(result_conv(y_ind,x_ind),[],1); 1177 1178 1179 % n=numel(X); 1180 % x=[X Y]; 1181 % X=X0.5; 1182 % Y=Y+0.5; 1183 % y = (X.*X+2*Y.*Y+X.*Y+6) + 0.1*rand(n,1); 1184 p = polyfitn(coord,result_conv,2); 1185 A(1,1)=2*p.Coefficients(1); 1186 A(1,2)=p.Coefficients(2); 1187 A(2,1)=p.Coefficients(2); 1188 A(2,2)=2*p.Coefficients(4); 1189 vector=[x y]'A\[p.Coefficients(3) p.Coefficients(5)]'; 1190 vector=vector'[floor(npx/2) floor(npy/2)]1 ; 1191 % zg = polyvaln(p,coord); 1192 % figure 1193 % surf(x_vec,y_vec,reshape(zg,9,9)) 1194 % hold on 1195 % plot3(X,Y,reshape(result_conv,9,9),'o') 1196 % hold off 1197 end 1152 1198 1153 1199 %'RUN_FIX': function for fixing velocity fields:
Note: See TracChangeset
for help on using the changeset viewer.