Changeset 596 for trunk/src/civ_matlab.m


Ignore:
Timestamp:
Mar 30, 2013, 11:34:27 AM (12 years ago)
Author:
sommeria
Message:

corrections done in civ

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/src/civ_matlab.m

    r594 r596  
    55%   filter_tps: make interpolation-smoothing
    66%------------------------------------------------------------------------
    7 % function [Data,errormsg,result_conv]= civ_uvmat(Param,ncfile)
     7% function [Data,errormsg,result_conv]= civ_matlab(Param,ncfile)
    88%
    99%OUTPUT
    … …  
    584584    iref=round(par_civ.Grid(ivec,1)+0.5);% xindex on the image A for the middle of the correlation box
    585585    jref=round(par_civ.ImageHeight-par_civ.Grid(ivec,2)+0.5);% yindex on the image B for the middle of the correlation box
    586     if ~(checkmask && par_civ.Mask(jref,iref)<=20) %velocity not set to zero by the black mask
    587 %         if jref-iby2<1 || jref+iby2>par_civ.ImageHeight|| iref-ibx2<1 || iref+ibx2>par_civ.ImageWidth||...
    588 %               jref+shifty(ivec)-isy2<1||jref+shifty(ivec)+isy2>par_civ.ImageHeight|| iref+shiftx(ivec)-isx2<1 || iref+shiftx(ivec)+isx2>par_civ.ImageWidth  % we are outside the image
    589 %             F(ivec)=3;
    590 %         else
    591             F(ivec)=0;
    592             subrange1_x=iref-ibx2:iref+ibx2;% x indices defining the first subimage
    593             subrange1_y=jref-iby2:jref+iby2;% y indices defining the first subimage
    594             subrange2_x=iref+shiftx(ivec)-isx2:iref+shiftx(ivec)+isx2;%x indices defining the second subimage
    595             subrange2_y=jref+shifty(ivec)-isy2:jref+shifty(ivec)+isy2;%y indices defining the second subimage
    596             image1_crop=MinA*ones(numel(subrange1_y),numel(subrange1_x));% default value=min of image A
    597             image2_crop=MinA*ones(numel(subrange2_y),numel(subrange2_x));% default value=min of image A
    598             check1_x=subrange1_x>=1 & subrange1_x<=par_civ.ImageWidth;% check which points in the subimage 1 are contained in the initial image 1
    599             check1_y=subrange1_y>=1 & subrange1_y<=par_civ.ImageHeight;
    600             check2_x=subrange2_x>=1 & subrange2_x<=par_civ.ImageWidth;% check which points in the subimage 2 are contained in the initial image 2
    601             check2_y=subrange2_y>=1 & subrange2_y<=par_civ.ImageHeight;
    602             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
    603             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
    604             image1_mean=mean(mean(image1_crop));
    605             image2_mean=mean(mean(image2_crop));
    606             %threshold on image minimum
    607             if check_MinIma && (image1_mean < par_civ.MinIma || image2_mean < par_civ.MinIma)
     586    %if ~(checkmask && par_civ.Mask(jref,iref)<=20) %velocity not set to zero by the black mask
     587    %         if jref-iby2<1 || jref+iby2>par_civ.ImageHeight|| iref-ibx2<1 || iref+ibx2>par_civ.ImageWidth||...
     588    %               jref+shifty(ivec)-isy2<1||jref+shifty(ivec)+isy2>par_civ.ImageHeight|| iref+shiftx(ivec)-isx2<1 || iref+shiftx(ivec)+isx2>par_civ.ImageWidth  % we are outside the image
     589    %             F(ivec)=3;
     590    %         else
     591    F(ivec)=0;
     592    subrange1_x=iref-ibx2:iref+ibx2;% x indices defining the first subimage
     593    subrange1_y=jref-iby2:jref+iby2;% y indices defining the first subimage
     594    subrange2_x=iref+shiftx(ivec)-isx2:iref+shiftx(ivec)+isx2;%x indices defining the second subimage
     595    subrange2_y=jref+shifty(ivec)-isy2:jref+shifty(ivec)+isy2;%y indices defining the second subimage
     596    image1_crop=MinA*ones(numel(subrange1_y),numel(subrange1_x));% default value=min of image A
     597    image2_crop=MinA*ones(numel(subrange2_y),numel(subrange2_x));% default value=min of image A
     598    check1_x=subrange1_x>=1 & subrange1_x<=par_civ.ImageWidth;% check which points in the subimage 1 are contained in the initial image 1
     599    check1_y=subrange1_y>=1 & subrange1_y<=par_civ.ImageHeight;
     600    check2_x=subrange2_x>=1 & subrange2_x<=par_civ.ImageWidth;% check which points in the subimage 2 are contained in the initial image 2
     601    check2_y=subrange2_y>=1 & subrange2_y<=par_civ.ImageHeight;
     602   
     603    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
     604    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
     605    image1_mean=mean(mean(image1_crop));
     606    image2_mean=mean(mean(image2_crop));
     607    %threshold on image minimum
     608    if check_MinIma && (image1_mean < par_civ.MinIma || image2_mean < par_civ.MinIma)
     609        F(ivec)=3;
     610    end
     611    %threshold on image maximum
     612    if check_MaxIma && (image1_mean > par_civ.MaxIma || image2_mean > par_civ.MaxIma)
     613        F(ivec)=3;
     614    end
     615    %         end
     616    if F(ivec)~=3
     617        image1_crop=image1_crop-image1_mean;%substract the mean
     618        image2_crop=image2_crop-image2_mean;
     619        if CheckDecimal
     620            xi=(1:mesh:size(image1_crop,2));
     621            yi=(1:mesh:size(image1_crop,1))';
     622            if CheckDeformation
     623                [XI,YI]=meshgrid(xi-ceil(size(image1_crop,2)/2),yi-ceil(size(image1_crop,1)/2));
     624                XIant=XI-par_civ.DUDX(ivec)*XI-par_civ.DUDY(ivec)*YI+ceil(size(image1_crop,2)/2);
     625                YIant=YI-par_civ.DVDX(ivec)*XI-par_civ.DVDY(ivec)*YI+ceil(size(image1_crop,1)/2);
     626                image1_crop=interp2(image1_crop,XIant,YIant);
     627            else
     628                image1_crop=interp2(image1_crop,xi,yi);
     629            end
     630            xi=(1:mesh:size(image2_crop,2));
     631            yi=(1:mesh:size(image2_crop,1))';
     632            image2_crop=interp2(image2_crop,xi,yi);
     633        end
     634        sum_square=sum(sum(image1_crop.*image1_crop));
     635        %reference: Oliver Pust, PIV: Direct Cross-Correlation
     636        result_conv= conv2(image2_crop,flipdim(flipdim(image1_crop,2),1),'valid');
     637        corrmax= max(max(result_conv));
     638        result_conv=(result_conv/corrmax)*255; %normalize, peak=always 255
     639        %Find the correlation max, at 255
     640        [y,x] = find(result_conv==255,1);
     641        if ~isempty(y) && ~isempty(x)
     642            try
     643                if par_civ.CorrSmooth==1
     644                    [vector,F(ivec)] = SUBPIXGAUSS (result_conv,x,y);
     645                elseif par_civ.CorrSmooth==2
     646                    [vector,F(ivec)] = SUBPIX2DGAUSS (result_conv,x,y);
     647                end
     648                utable(ivec)=vector(1)*mesh+shiftx(ivec);
     649                vtable(ivec)=vector(2)*mesh+shifty(ivec);
     650                xtable(ivec)=iref+utable(ivec)/2-0.5;% convec flow (velocity taken at the point middle from imgae 1 and 2)
     651                ytable(ivec)=jref+vtable(ivec)/2-0.5;% and position of pixel 1=0.5 (convention for image coordinates=0 at the edge)
     652                iref=round(xtable(ivec));% image index for the middle of the vector
     653                jref=round(ytable(ivec));
     654                if checkmask && par_civ.Mask(jref,iref)<200 && par_civ.Mask(jref,iref)>=100
     655                    utable(ivec)=0;
     656                    vtable(ivec)=0;
     657                    F(ivec)=3;
     658                end
     659                ctable(ivec)=corrmax/sum_square;% correlation value
     660            catch ME
    608661                F(ivec)=3;
    609662            end
    610             %threshold on image maximum
    611             if check_MaxIma && (image1_mean > par_civ.MaxIma || image2_mean > par_civ.MaxIma)
    612                 F(ivec)=3;
    613             end
    614 %         end     
    615         if F(ivec)~=3
    616             image1_crop=image1_crop-image1_mean;%substract the mean
    617             image2_crop=image2_crop-image2_mean;
    618             if CheckDecimal
    619                 xi=(1:mesh:size(image1_crop,2));
    620                 yi=(1:mesh:size(image1_crop,1))';
    621                 if CheckDeformation
    622                     [XI,YI]=meshgrid(xi-ceil(size(image1_crop,2)/2),yi-ceil(size(image1_crop,1)/2));
    623                     XIant=XI-par_civ.DUDX(ivec)*XI-par_civ.DUDY(ivec)*YI+ceil(size(image1_crop,2)/2);
    624                     YIant=YI-par_civ.DVDX(ivec)*XI-par_civ.DVDY(ivec)*YI+ceil(size(image1_crop,1)/2);
    625                     image1_crop=interp2(image1_crop,XIant,YIant);
    626                 else
    627                     image1_crop=interp2(image1_crop,xi,yi);
    628                 end
    629                 xi=(1:mesh:size(image2_crop,2));
    630                 yi=(1:mesh:size(image2_crop,1))';
    631                 image2_crop=interp2(image2_crop,xi,yi);
    632             end
    633             sum_square=sum(sum(image1_crop.*image1_crop));
    634             %reference: Oliver Pust, PIV: Direct Cross-Correlation
    635             result_conv= conv2(image2_crop,flipdim(flipdim(image1_crop,2),1),'valid');
    636             corrmax= max(max(result_conv));
    637             result_conv=(result_conv/corrmax)*255; %normalize, peak=always 255
    638             %Find the correlation max, at 255
    639             [y,x] = find(result_conv==255,1);
    640             if ~isempty(y) && ~isempty(x)
    641                 try
    642                     if par_civ.CorrSmooth==1
    643                         [vector,F(ivec)] = SUBPIXGAUSS (result_conv,x,y);
    644                     elseif par_civ.CorrSmooth==2
    645                         [vector,F(ivec)] = SUBPIX2DGAUSS (result_conv,x,y);
    646                     end
    647                     utable(ivec)=vector(1)*mesh+shiftx(ivec);
    648                     vtable(ivec)=vector(2)*mesh+shifty(ivec);                 
    649                     xtable(ivec)=iref+utable(ivec)/2-0.5;% convec flow (velocity taken at the point middle from imgae 1 and 2)
    650                     ytable(ivec)=jref+vtable(ivec)/2-0.5;% and position of pixel 1=0.5 (convention for image coordinates=0 at the edge)
    651                     iref=round(xtable(ivec));% image index for the middle of the vector
    652                     jref=round(ytable(ivec));
    653                     if checkmask && par_civ.Mask(jref,iref)<200 && par_civ.Mask(jref,iref)>=100
    654                         utable(ivec)=0;
    655                         vtable(ivec)=0;
    656                         F(ivec)=3;
    657                     end
    658                     ctable(ivec)=corrmax/sum_square;% correlation value
    659                 catch ME
    660                     F(ivec)=3;
    661                 end
    662             else
    663                 F(ivec)=3;
    664             end
    665         end
    666     end
    667    
    668     %Create the vector matrix x, y, u, v
     663        else
     664            F(ivec)=3;
     665        end
     666    end
    669667end
    670668result_conv=result_conv*corrmax/(255*sum_square);% keep the last correlation matrix for output
Note: See TracChangeset for help on using the changeset viewer.