Changeset 879 for trunk/src/series/stereo_civ.m
 Timestamp:
 Feb 25, 2015, 12:02:38 AM (9 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

trunk/src/series/stereo_civ.m
r878 r879 903 903 nbvec=size(par_civ.Grid,1); 904 904 905 %% prepare correlation and search boxes 905 %% prepare correlation and search boxes 906 906 ibx2=floor(par_civ.CorrBoxSize(1)/2); 907 907 iby2=floor(par_civ.CorrBoxSize(2)/2); … … 945 945 946 946 %% Apply mask 947 948 949 950 951 % 100>=mask> 20: velocity not calculated, impermeable (no flux through mask boundaries) 952 947 % Convention for mask IDEAS TO IMPLEMENT ? 948 % mask >200 : velocity calculated 949 % 200 >=mask>150;velocity not calculated, interpolation allowed (bad spots) 950 % 150>=mask >100: velocity not calculated, nor interpolated 951 % 100>=mask> 20: velocity not calculated, impermeable (no flux through mask boundaries) 952 % 20>=mask: velocity=0 953 953 checkmask=0; 954 954 MinA=min(min(par_civ.ImageA)); 955 MinB=min(min(par_civ.ImageB)); 955 %MinB=min(min(par_civ.ImageB)); 956 %check_undefined=false(size(par_civ.ImageA)); 956 957 if isfield(par_civ,'Mask') && ~isempty(par_civ.Mask) 957 checkmask=1;958 if ~isequal(size(par_civ.Mask),[npy_ima npx_ima])958 checkmask=1; 959 if ~isequal(size(par_civ.Mask),[npy_ima npx_ima]) 959 960 errormsg='mask must be an image with the same size as the images'; 960 961 return 961 end962 % check_noflux=(par_civ.Mask<100) ;%TODO: to implement962 end 963 % check_noflux=(par_civ.Mask<100) ;%TODO: to implement 963 964 check_undefined=(par_civ.Mask<200 & par_civ.Mask>=20 ); 964 % par_civ.ImageA(check_undefined)=MinA;% put image A to zero (i.e. the min image value) in the undefined area965 % par_civ.ImageB(check_undefined)=MinB;% put image B to zero (i.e. the min image value) in the undefined area965 % par_civ.ImageA(check_undefined)=MinA;% put image A to zero (i.e. the min image value) in the undefined area 966 % par_civ.ImageB(check_undefined)=MinB;% put image B to zero (i.e. the min image value) in the undefined area 966 967 end 967 968 … … 997 998 check1_y=subrange1_y>=1 & subrange1_y<=par_civ.ImageHeight; 998 999 check2_x=subrange2_x>=1 & subrange2_x<=par_civ.ImageWidth;% check which points in the subimage 2 are contained in the initial image 2 999 check2_y=subrange2_y>=1 & subrange2_y<=par_civ.ImageHeight; 1000 check2_y=subrange2_y>=1 & subrange2_y<=par_civ.ImageHeight; 1000 1001 image1_crop(check1_y,check1_x)=par_civ.ImageA(subrange1_y(check1_y),subrange1_x(check1_x));%extract a subimage (correlation box) from image A 1001 1002 image2_crop(check2_y,check2_x)=par_civ.ImageB(subrange2_y(check2_y),subrange2_x(check2_x));%extract a larger subimage (search box) from image B 1002 mask1_crop(check1_y,check1_x)=check_undefined(subrange1_y(check1_y),subrange1_x(check1_x));%extract a mask subimage (correlation box) from image A 1003 mask2_crop(check2_y,check2_x)=check_undefined(subrange2_y(check2_y),subrange2_x(check2_x));%extract a mask subimage (search box) from imag 1004 sizemask=sum(sum(mask1_crop))/(numel(subrange1_y)*numel(subrange1_x));%size of the masked part relative to the correlation subimage 1003 if checkmask 1004 mask1_crop(check1_y,check1_x)=check_undefined(subrange1_y(check1_y),subrange1_x(check1_x));%extract a mask subimage (correlation box) from image A 1005 mask2_crop(check2_y,check2_x)=check_undefined(subrange2_y(check2_y),subrange2_x(check2_x));%extract a mask subimage (search box) from imag 1006 sizemask=sum(sum(mask1_crop))/(numel(subrange1_y)*numel(subrange1_x));%size of the masked part relative to the correlation subimage 1005 1007 if sizemask > 1/2% eliminate point if more than half of the correlation box is masked 1006 1008 F(ivec)=3; % … … 1010 1012 image1_mean=mean(mean(image1_crop))/(1sizemask); 1011 1013 image2_mean=mean(mean(image2_crop))/(1sizemask); 1012 % image1_mean=mean(mean(image1_crop)); 1013 % image2_mean=mean(mean(image2_crop)); 1014 end 1015 else 1016 image1_mean=mean(mean(image1_crop)); 1017 image2_mean=mean(mean(image2_crop)); 1018 end 1014 1019 %threshold on image minimum 1015 1020 if check_MinIma && (image1_mean < par_civ.MinIma  image2_mean < par_civ.MinIma) … … 1020 1025 F(ivec)=3; 1021 1026 end 1022 end1023 1027 if F(ivec)~=3 1024 image1_crop=(image1_cropimage1_mean).*~mask1_crop;%substract the mean, put to zero the masked parts 1025 image2_crop=(image2_cropimage2_mean).*~mask2_crop; 1028 image1_crop=(image1_cropimage1_mean);%substract the mean, put to zero the masked parts 1029 image2_crop=(image2_cropimage2_mean); 1030 if checkmask 1031 image1_crop=image1_crop.*~mask1_crop;% put to zero the masked parts 1032 image2_crop=image2_crop.*~mask2_crop; 1033 end 1026 1034 if CheckDeformation 1027 1028 1029 1030 1031 1032 1033 1034 1035 1036 1037 1035 xi=(1:mesh:size(image1_crop,2)); 1036 yi=(1:mesh:size(image1_crop,1))'; 1037 [XI,YI]=meshgrid(xiceil(size(image1_crop,2)/2),yiceil(size(image1_crop,1)/2)); 1038 XIant=XIpar_civ.DUDX(ivec)*XIpar_civ.DUDY(ivec)*YI+ceil(size(image1_crop,2)/2); 1039 YIant=YIpar_civ.DVDX(ivec)*XIpar_civ.DVDY(ivec)*YI+ceil(size(image1_crop,1)/2); 1040 image1_crop=interp2(image1_crop,XIant,YIant); 1041 image1_crop(isnan(image1_crop))=0; 1042 xi=(1:mesh:size(image2_crop,2)); 1043 yi=(1:mesh:size(image2_crop,1))'; 1044 image2_crop=interp2(image2_crop,xi,yi,'*spline'); 1045 image2_crop(isnan(image2_crop))=0; 1038 1046 end 1039 1047 sum_square=(sum(sum(image1_crop.*image1_crop)));%+sum(sum(image2_crop.*image2_crop)))/2; … … 1045 1053 [y,x] = find(result_conv==255,1); 1046 1054 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 1047 1048 1055 sum_square=sum_square*sum(sum(subimage2_crop.*subimage2_crop));% product of variances of image 1 and 2 1056 sum_square=sqrt(sum_square);% sqrt of the variance product to normalise correlation 1049 1057 if ~isempty(y) && ~isempty(x) 1050 1058 try … … 1055 1063 end 1056 1064 1057 1065 1058 1066 utable(ivec)=vector(1)*mesh+shiftx(ivec); 1059 1067 vtable(ivec)=vector(2)*mesh+shifty(ivec); 1060 1068 1061 1069 1062 1070 xtable(ivec)=iref+utable(ivec)/20.5;% convec flow (velocity taken at the point middle from imgae 1 and 2)
Note: See TracChangeset
for help on using the changeset viewer.