Changeset 879


Ignore:
Timestamp:
Feb 25, 2015, 12:02:38 AM (10 years ago)
Author:
sommeria
Message:

bugs corrected in stereo_civ

Location:
trunk/src
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • trunk/src/get_field.m

    r876 r879  
    788788        set(handles.Coord_x,'Value',coord_val(1))
    789789        set(handles.Coord_y,'Value',coord_val(2))
    790         if numel(ListCoord)>=3
     790        if numel(coord_val)>=3
    791791            set(handles.Coord_z,'Value',coord_val(3))
    792792        end

  • trunk/src/series/civ_series.m

    r878 r879  
    958958MinA=min(min(par_civ.ImageA));
    959959%MinB=min(min(par_civ.ImageB));
    960 check_undefined=false(size(par_civ.ImageA));
     960%check_undefined=false(size(par_civ.ImageA));
    961961if isfield(par_civ,'Mask') && ~isempty(par_civ.Mask)
    962962    checkmask=1;
    … …  
    990990        image1_crop=MinA*ones(numel(subrange1_y),numel(subrange1_x));% default value=min of image A
    991991        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 mask
    993         mask2_crop=ones(numel(subrange2_y),numel(subrange2_x));% default value=1 for mask
    994992        check1_x=subrange1_x>=1 & subrange1_x<=par_civ.ImageWidth;% check which points in the subimage 1 are contained in the initial image 1
    995993        check1_y=subrange1_y>=1 & subrange1_y<=par_civ.ImageHeight;
    … …  
    998996        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
    999997        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
    10051012        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       
    10191025        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
    10221033            if CheckDeformation
    10231034                xi=(1:mesh:size(image1_crop,2));

  • trunk/src/series/stereo_civ.m

    r878 r879  
    903903nbvec=size(par_civ.Grid,1);
    904904
    905 %% prepare correlation and search boxes 
     905%% prepare correlation and search boxes
    906906ibx2=floor(par_civ.CorrBoxSize(1)/2);
    907907iby2=floor(par_civ.CorrBoxSize(2)/2);
    … …  
    945945
    946946%% Apply mask
    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
     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
    953953checkmask=0;
    954954MinA=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));
    956957if 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])
    959960        errormsg='mask must be an image with the same size as the images';
    960961        return
    961    end
    962   %  check_noflux=(par_civ.Mask<100) ;%TODO: to implement
     962    end
     963    %  check_noflux=(par_civ.Mask<100) ;%TODO: to implement
    963964    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  area
    965 %     par_civ.ImageB(check_undefined)=MinB;% put image B to zero (i.e. the min image value) in the undefined  area
     965    %     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
    966967end
    967968
    … …  
    997998    check1_y=subrange1_y>=1 & subrange1_y<=par_civ.ImageHeight;
    998999    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;
    10001001    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
    10011002    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
    10051007        if sizemask > 1/2% eliminate point if more than half of the correlation box is masked
    10061008            F(ivec)=3; %
    … …  
    10101012            image1_mean=mean(mean(image1_crop))/(1-sizemask);
    10111013            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
    10141019    %threshold on image minimum
    10151020    if check_MinIma && (image1_mean < par_civ.MinIma || image2_mean < par_civ.MinIma)
    … …  
    10201025        F(ivec)=3;
    10211026    end
    1022     end
    10231027    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
    10261034        if CheckDeformation
    1027               xi=(1:mesh:size(image1_crop,2));
    1028                 yi=(1:mesh:size(image1_crop,1))';
    1029                 [XI,YI]=meshgrid(xi-ceil(size(image1_crop,2)/2),yi-ceil(size(image1_crop,1)/2));
    1030                 XIant=XI-par_civ.DUDX(ivec)*XI-par_civ.DUDY(ivec)*YI+ceil(size(image1_crop,2)/2);
    1031                 YIant=YI-par_civ.DVDX(ivec)*XI-par_civ.DVDY(ivec)*YI+ceil(size(image1_crop,1)/2);
    1032                 image1_crop=interp2(image1_crop,XIant,YIant);
    1033                 image1_crop(isnan(image1_crop))=0;
    1034                 xi=(1:mesh:size(image2_crop,2));
    1035                 yi=(1:mesh:size(image2_crop,1))';
    1036                 image2_crop=interp2(image2_crop,xi,yi,'*spline');
    1037                 image2_crop(isnan(image2_crop))=0;
     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;
    10381046        end
    10391047        sum_square=(sum(sum(image1_crop.*image1_crop)));%+sum(sum(image2_crop.*image2_crop)))/2;
    … …  
    10451053        [y,x] = find(result_conv==255,1);
    10461054        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             sum_square=sum_square*sum(sum(subimage2_crop.*subimage2_crop));% product of variances of image 1 and 2
    1048             sum_square=sqrt(sum_square);% sqrt of the variance product to normalise correlation
     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
    10491057        if ~isempty(y) && ~isempty(x)
    10501058            try
    … …  
    10551063                end
    10561064               
    1057 
     1065               
    10581066                utable(ivec)=vector(1)*mesh+shiftx(ivec);
    10591067                vtable(ivec)=vector(2)*mesh+shifty(ivec);
    1060 
     1068               
    10611069               
    10621070                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.