Changeset 596 for trunk/src/civ_matlab.m
 Timestamp:
 Mar 30, 2013, 11:34:27 AM (11 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

trunk/src/civ_matlab.m
r594 r596 5 5 % filter_tps: make interpolationsmoothing 6 6 % 7 % function [Data,errormsg,result_conv]= civ_ uvmat(Param,ncfile)7 % function [Data,errormsg,result_conv]= civ_matlab(Param,ncfile) 8 8 % 9 9 %OUTPUT … … 584 584 iref=round(par_civ.Grid(ivec,1)+0.5);% xindex on the image A for the middle of the correlation box 585 585 jref=round(par_civ.ImageHeightpar_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 jrefiby2<1  jref+iby2>par_civ.ImageHeight irefibx2<1  iref+ibx2>par_civ.ImageWidth... 588 % jref+shifty(ivec)isy2<1jref+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=irefibx2:iref+ibx2;% x indices defining the first subimage 593 subrange1_y=jrefiby2: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 jrefiby2<1  jref+iby2>par_civ.ImageHeight irefibx2<1  iref+ibx2>par_civ.ImageWidth... 588 % jref+shifty(ivec)isy2<1jref+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=irefibx2:iref+ibx2;% x indices defining the first subimage 593 subrange1_y=jrefiby2: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_cropimage1_mean;%substract the mean 618 image2_crop=image2_cropimage2_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(xiceil(size(image1_crop,2)/2),yiceil(size(image1_crop,1)/2)); 624 XIant=XIpar_civ.DUDX(ivec)*XIpar_civ.DUDY(ivec)*YI+ceil(size(image1_crop,2)/2); 625 YIant=YIpar_civ.DVDX(ivec)*XIpar_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 CrossCorrelation 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)/20.5;% convec flow (velocity taken at the point middle from imgae 1 and 2) 651 ytable(ivec)=jref+vtable(ivec)/20.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 608 661 F(ivec)=3; 609 662 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_cropimage1_mean;%substract the mean 617 image2_crop=image2_cropimage2_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(xiceil(size(image1_crop,2)/2),yiceil(size(image1_crop,1)/2)); 623 XIant=XIpar_civ.DUDX(ivec)*XIpar_civ.DUDY(ivec)*YI+ceil(size(image1_crop,2)/2); 624 YIant=YIpar_civ.DVDX(ivec)*XIpar_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 CrossCorrelation 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)/20.5;% convec flow (velocity taken at the point middle from imgae 1 and 2) 650 ytable(ivec)=jref+vtable(ivec)/20.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 669 667 end 670 668 result_conv=result_conv*corrmax/(255*sum_square);% keep the last correlation matrix for output
Note: See TracChangeset
for help on using the changeset viewer.