Ignore:
Timestamp:
Feb 22, 2015, 11:07:07 PM (9 years ago)
Author:
sommeria
Message:

deformation modified in civ_series and bug corrections

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/src/series/civ_series.m

    r873 r876  
    691691            end
    692692        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 guess
    696         %par_civ2.SearchBoxSize(2)=2*iby2+9;
    697693        if CheckInputFile % else Dt given by par_civ2
    698694            if strcmp(Param.ActionInput.ListCompareMode,'displacement')
     
    707703        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];
    708704        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);
    713709        end
    714710       
     
    910906
    911907%% 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);
     908ibx2=floor(par_civ.CorrBoxSize(1)/2);
     909iby2=floor(par_civ.CorrBoxSize(2)/2);
     910isx2=floor(par_civ.SearchBoxSize(1)/2);
     911isy2=floor(par_civ.SearchBoxSize(2)/2);
    916912shiftx=round(par_civ.SearchBoxShift(:,1));%use the input shift estimate, rounded to the next integer value
    917913shifty=-round(par_civ.SearchBoxShift(:,2));% sign minus because image j index increases when y decreases
     
    968964    end
    969965    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  area
    971     par_civ.ImageB(check_undefined)=0;% put image B to zero (i.e. the min image value) in the undefined  area
     966%     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
    972968end
    973969
     
    993989        image2_crop=MinA*ones(numel(subrange2_y),numel(subrange2_x));% default value=min of image A
    994990        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=min for mask
     991        mask2_crop=ones(numel(subrange2_y),numel(subrange2_x));% default value=1 for mask
    996992        check1_x=subrange1_x>=1 & subrange1_x<=par_civ.ImageWidth;% check which points in the subimage 1 are contained in the initial image 1
    997993        check1_y=subrange1_y>=1 & subrange1_y<=par_civ.ImageHeight;
     
    10061002            F(ivec)=3; %
    10071003        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;
    10081006            image1_mean=mean(mean(image1_crop))/(1-sizemask);
    10091007            image2_mean=mean(mean(image2_crop))/(1-sizemask);
     
    10261024                XIant=XI-par_civ.DUDX(ivec)*XI-par_civ.DUDY(ivec)*YI+ceil(size(image1_crop,2)/2);
    10271025                YIant=YI-par_civ.DVDX(ivec)*XI-par_civ.DVDY(ivec)*YI+ceil(size(image1_crop,1)/2);
     1026                        image1_crop_1=image1_crop;
    10281027                image1_crop=interp2(image1_crop,XIant,YIant);
     1028                        image1_crop_2=image1_crop;
    10291029                image1_crop(isnan(image1_crop))=0;
     1030                  image1_crop_3=image1_crop;
    10301031                xi=(1:mesh:size(image2_crop,2));
    10311032                yi=(1:mesh:size(image2_crop,1))';
    10321033                image2_crop=interp2(image2_crop,xi,yi);
    10331034                image2_crop(isnan(image2_crop))=0;
     1035                par_civ.CorrSmooth=3;
    10341036                %%
    10351037            end
     
    10411043            %Find the correlation max, at 255
    10421044            [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
    10431048            if ~isempty(y) && ~isempty(x)
    10441049                try
     
    10471052                    elseif par_civ.CorrSmooth==2
    10481053                        [vector,F(ivec)] = SUBPIX2DGAUSS (result_conv,x,y);
     1054                    else
     1055                        [vector,F(ivec)] = quadr_fit(result_conv,x,y);
    10491056                    end
    10501057                    utable(ivec)=vector(1)*mesh+shiftx(ivec);
     
    10521059                    xtable(ivec)=iref+utable(ivec)/2-0.5;% convec flow (velocity taken at the point middle from imgae 1 and 2)
    10531060                    ytable(ivec)=jref+vtable(ivec)/2-0.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 vector
     1061                    iref=round(xtable(ivec));% nearest image index for the middle of the vector
    10551062                    jref=round(ytable(ivec));
     1063                    % eliminate vectors located in the mask
    10561064                    if checkmask && par_civ.Mask(jref,iref)<200 && par_civ.Mask(jref,iref)>=100
    10571065                        utable(ivec)=0;
     
    11501158end
    11511159vector=[peakx-floor(npx/2)-1 peaky-floor(npy/2)-1];
     1160
     1161%------------------------------------------------------------------------
     1162% --- Find the maximum of the correlation function after quadratic interpolation
     1163function [vector,F] = quadr_fit(result_conv,x,y)
     1164[npy,npx]=size(result_conv);
     1165if x<4 || y<4 || npx-x<4 ||npy-y <4
     1166    F=-2;
     1167    vector=[x y];
     1168else
     1169    F=0;
     1170    x_ind=x-4:x+4;
     1171    y_ind=y-4:y+4;
     1172    x_vec=0.25*(x_ind-x);
     1173    y_vec=0.25*(y_ind-y);
     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=X-0.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
     1197end
    11521198
    11531199%'RUN_FIX': function for fixing velocity fields:
Note: See TracChangeset for help on using the changeset viewer.