Changeset 879 for trunk/src/series
- Timestamp:
- Feb 25, 2015, 12:02:38 AM (10 years ago)
- Location:
- trunk/src/series
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/src/series/civ_series.m
r878 r879 958 958 MinA=min(min(par_civ.ImageA)); 959 959 %MinB=min(min(par_civ.ImageB)); 960 check_undefined=false(size(par_civ.ImageA));960 %check_undefined=false(size(par_civ.ImageA)); 961 961 if isfield(par_civ,'Mask') && ~isempty(par_civ.Mask) 962 962 checkmask=1; … … 990 990 image1_crop=MinA*ones(numel(subrange1_y),numel(subrange1_x));% default value=min of image A 991 991 image2_crop=MinA*ones(numel(subrange2_y),numel(subrange2_x));% default value=min of image A 992 mask1_crop=ones(numel(subrange1_y),numel(subrange1_x));% default value=1 for mask993 mask2_crop=ones(numel(subrange2_y),numel(subrange2_x));% default value=1 for mask994 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 995 993 check1_y=subrange1_y>=1 & subrange1_y<=par_civ.ImageHeight; … … 998 996 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 999 997 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 1000 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 1001 mask2_crop(check2_y,check2_x)=check_undefined(subrange2_y(check2_y),subrange2_x(check2_x));%extract a mask subimage (search box) from image B 1002 sizemask=sum(sum(mask1_crop))/(numel(subrange1_y)*numel(subrange1_x));%size of the masked part relative to the correlation sub-image 1003 if sizemask > 1/2% eliminate point if more than half of the correlation box is masked 1004 F(ivec)=3; % 998 if checkmask 999 mask1_crop=ones(numel(subrange1_y),numel(subrange1_x));% default value=1 for mask 1000 mask2_crop=ones(numel(subrange2_y),numel(subrange2_x));% default value=1 for mask 1001 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 1002 mask2_crop(check2_y,check2_x)=check_undefined(subrange2_y(check2_y),subrange2_x(check2_x));%extract a mask subimage (search box) from image B 1003 sizemask=sum(sum(mask1_crop))/(numel(subrange1_y)*numel(subrange1_x));%size of the masked part relative to the correlation sub-image 1004 if sizemask > 1/2% eliminate point if more than half of the correlation box is masked 1005 F(ivec)=3; % 1006 else 1007 image1_crop=image1_crop.*~mask1_crop;% put to zero the masked pixels (mask1_crop='true'=1) 1008 image2_crop=image2_crop.*~mask2_crop; 1009 image1_mean=mean(mean(image1_crop))/(1-sizemask); 1010 image2_mean=mean(mean(image2_crop))/(1-sizemask); 1011 end 1005 1012 else 1006 image1_crop=image1_crop.*~mask1_crop;% put to zero the masked pixels (mask1_crop='true'=1) 1007 image2_crop=image2_crop.*~mask2_crop; 1008 image1_mean=mean(mean(image1_crop))/(1-sizemask); 1009 image2_mean=mean(mean(image2_crop))/(1-sizemask); 1010 %threshold on image minimum 1011 if check_MinIma && (image1_mean < par_civ.MinIma || image2_mean < par_civ.MinIma) 1012 F(ivec)=3; 1013 end 1014 %threshold on image maximum 1015 if check_MaxIma && (image1_mean > par_civ.MaxIma || image2_mean > par_civ.MaxIma) 1016 F(ivec)=3; 1017 end 1018 end 1013 image1_mean=mean(mean(image1_crop)); 1014 image2_mean=mean(mean(image2_crop)); 1015 end 1016 %threshold on image minimum 1017 if check_MinIma && (image1_mean < par_civ.MinIma || image2_mean < par_civ.MinIma) 1018 F(ivec)=3; 1019 end 1020 %threshold on image maximum 1021 if check_MaxIma && (image1_mean > par_civ.MaxIma || image2_mean > par_civ.MaxIma) 1022 F(ivec)=3; 1023 end 1024 1019 1025 if F(ivec)~=3 1020 image1_crop=(image1_crop-image1_mean).*~mask1_crop;%substract the mean, put to zero the masked parts 1021 image2_crop=(image2_crop-image2_mean).*~mask2_crop; 1026 if checkmask 1027 image1_crop=(image1_crop-image1_mean).*~mask1_crop;%substract the mean, put to zero the masked parts 1028 image2_crop=(image2_crop-image2_mean).*~mask2_crop; 1029 else 1030 image1_crop=(image1_crop-image1_mean); 1031 image2_crop=(image2_crop-image2_mean); 1032 end 1022 1033 if CheckDeformation 1023 1034 xi=(1:mesh:size(image1_crop,2)); -
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 sub-image 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 sub-image 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))/(1-sizemask); 1011 1013 image2_mean=mean(mean(image2_crop))/(1-sizemask); 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_crop-image1_mean).*~mask1_crop;%substract the mean, put to zero the masked parts 1025 image2_crop=(image2_crop-image2_mean).*~mask2_crop; 1028 image1_crop=(image1_crop-image1_mean);%substract the mean, put to zero the masked parts 1029 image2_crop=(image2_crop-image2_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(xi-ceil(size(image1_crop,2)/2),yi-ceil(size(image1_crop,1)/2)); 1038 XIant=XI-par_civ.DUDX(ivec)*XI-par_civ.DUDY(ivec)*YI+ceil(size(image1_crop,2)/2); 1039 YIant=YI-par_civ.DVDX(ivec)*XI-par_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)/2-0.5;% convec flow (velocity taken at the point middle from imgae 1 and 2)
Note: See TracChangeset
for help on using the changeset viewer.