source: trunk/src/pivlab.m @ 239

Last change on this file since 239 was 238, checked in by sommeria, 14 years ago

various bug corrections

File size: 7.3 KB
RevLine 
[225]1% 'pivlab': function piv.m adapted from PIVlab http://pivlab.blogspot.com/
2%--------------------------------------------------------------------------
[227]3% function [xtable ytable utable vtable typevector] = pivlab (image1,image2,ibx,iby step, subpixfinder, mask, roi)
[225]4%
5% OUTPUT:
6% xtable: set of x coordinates
7% ytable: set of y coordiantes
8% utable: set of u displacements (along x)
9% vtable: set of v displacements (along y)
10% ctable: max image correlation for each vector
11% typevector: set of flags, =1 for good, =0 for NaN vectors
12%
13%INPUT:
14% image1:first image (matrix)
15% image2: second image (matrix)
[227]16% ibx,iby: size of the correlation box along x and y (in px)
[225]17% step: mesh of the measurement points (in px)
18% subpixfinder=1 or 2 controls the curve fitting of the image correlation
19% mask: =[] for no mask
20% roi: 4 element vector defining a region of interest: x position, y position, width, height, (in image indices), for the whole image, roi=[];
[236]21function [xtable ytable utable vtable ctable F result_conv errormsg] = pivlab (image1,image2,ibx2,iby2,isx2,isy2,shiftx,shifty, GridIndices, subpixfinder,mask)
[225]22%this funtion performs the DCC PIV analysis. Recent window-deformation
23%methods perform better and will maybe be implemented in the future.
[236]24nbvec=size(GridIndices,1);
25xtable=zeros(nbvec,1);
26ytable=xtable;
27utable=xtable;
28vtable=xtable;
29ctable=xtable;
30F=xtable;
31result_conv=[];
[231]32errormsg='';
[236]33%warning off %MATLAB:log:logOfZero
[231]34[npy_ima npx_ima]=size(image1);
[236]35if ~isequal(size(image2),[npy_ima npx_ima])
[231]36    errormsg='image pair with unequal size';
37    return
38end
[236]39
40%% mask
41testmask=0;
42if exist('mask','var') && ~isempty(mask)
43   testmask=1;
44   if ~isequal(size(mask),[npy_ima npx_ima])
45        errormsg='mask must be an image with the same size as the images';
46        return
47   end
48    % Convention for mask
49    % mask >200 : velocity calculated
50    %  200 >=mask>150;velocity not calculated, interpolation allowed (bad spots)
51    % 150>=mask >100: velocity not calculated, nor interpolated
52    %  100>=mask> 20: velocity not calculated, impermeable (no flux through mask boundaries)
53    %  20>=mask: velocity=0
54    test_noflux=(mask<=100) ;
55    test_undefined=(mask<=200 & mask>100 );
56    image1(test_undefined)=min(min(image1))*ones(size(image1));% put image to zero in the undefined  area
57    image2(test_undefined)=min(min(image1))*ones(size(image1));% put image to zero in the undefined  area
[225]58end
[236]59image1=double(image1);
60image2=double(image2);
[231]61
[236]62%% calculate correlations: MAINLOOP
63corrmax=0;
64sum_square=1;% default
[227]65for ivec=1:nbvec
[231]66    iref=GridIndices(ivec,1);
67    jref=GridIndices(ivec,2);
[236]68    testmask_ij=0;
69    test0=0;
70    if testmask
71        if mask(jref,iref)<=20
72           vector=[0 0];
73           test0=1;
74        else
75            mask_crop1=mask(jref-iby2:jref+iby2,iref-ibx2:iref+ibx2);
76            mask_crop2=mask(jref+shifty-isy2:jref+shifty+isy2,iref+shiftx-isx2:iref+shiftx+isx2);
77            if ~isempty(find(mask_crop1<=200 & mask_crop1>100,1)) || ~isempty(find(mask_crop2<=200 & mask_crop2>100,1));
78                testmask_ij=1;
79            end
80        end
81    end
82    if ~test0
83        image1_crop=image1(jref-iby2:jref+iby2,iref-ibx2:iref+ibx2);
84        image2_crop=image2(jref+shifty-isy2:jref+shifty+isy2,iref+shiftx-isx2:iref+shiftx+isx2);
85        image1_crop=image1_crop-mean(mean(image1_crop));
86        image2_crop=image2_crop-mean(mean(image2_crop));
87        %reference: Oliver Pust, PIV: Direct Cross-Correlation
88        result_conv= conv2(image2_crop,flipdim(flipdim(image1_crop,2),1),'valid');
89        corrmax= max(max(result_conv));
90        result_conv=(result_conv/corrmax)*255; %normalize, peak=always 255
91        %Find the correlation max, at 255
[238]92        [y,x] = find(result_conv==255,1);
93        if ~isnan(y) && ~isnan(x)
[236]94            try
95                if subpixfinder==1
[238]96                    [vector,F(ivec)] = SUBPIXGAUSS (result_conv,x,y);
[236]97                elseif subpixfinder==2
[238]98                    [vector,F(ivec)] = SUBPIX2DGAUSS (result_conv,x,y);
[225]99                end
[236]100                sum_square=sum(sum(image1_crop.*image1_crop));
101                ctable(ivec)=corrmax/sum_square;% correlation value
[238]102%                 if vector(1)>shiftx+isx2-ibx2+subpixfinder || vector(2)>shifty+isy2-iby2+subpixfinder
103%                     F(ivec)=-2;%vector reaches the border of the search zone
104%                 end
[236]105            catch ME
[231]106                vector=[0 0]; %if something goes wrong with cross correlation.....
[236]107                F(ivec)=3;
[225]108            end
[236]109        else
110            vector=[0 0]; %if something goes wrong with cross correlation.....
111            F(ivec)=3;
[225]112        end
[236]113        if testmask_ij
114            F(ivec)=3;
115        end
116    end
117   
118    %Create the vector matrix x, y, u, v
119    xtable(ivec)=iref+vector(1)/2;% convec flow (velocity taken at the point middle from imgae1 and 2)
120    ytable(ivec)=jref+vector(2)/2;
121    utable(ivec)=vector(1);
122    vtable(ivec)=vector(2);
[225]123end
[236]124result_conv=result_conv*corrmax/(255*sum_square);% keep the last correlation matrix for output
[225]125
[231]126
[238]127function [vector,F] = SUBPIXGAUSS (result_conv,x,y)
128vector=[0 0]; %default
129F=0;
130[npy,npx]=size(result_conv);
[231]131
[238]132% if (x <= (size(result_conv,1)-1)) && (y <= (size(result_conv,1)-1)) && (x >= 1) && (y >= 1)
[225]133    %the following 8 lines are copyright (c) 1998, Uri Shavit, Roi Gurka, Alex Liberzon, Technion – Israel Institute of Technology
134    %http://urapiv.wordpress.com
[238]135    peaky = y;
136    if y <= npy-1 && y >= 1
137        f0 = log(result_conv(y,x));
138        f1 = log(result_conv(y-1,x));
139        f2 = log(result_conv(y+1,x));
140        peaky = peaky+ (f1-f2)/(2*f1-4*f0+2*f2);
141    else
142        F=-2; % warning flag for vector truncated by the limited search box
143    end
144    peakx=x;
145    if x <= npx-1 && x >= 1
146        f0 = log(result_conv(y,x));
147        f1 = log(result_conv(y,x-1));
148        f2 = log(result_conv(y,x+1));
149        peakx = peakx+ (f1-f2)/(2*f1-4*f0+2*f2);
150    else
151        F=-2; % warning flag for vector truncated by the limited search box
152    end
[231]153    vector=[peakx-floor(npx/2)-1 peaky-floor(npy/2)-1];
[238]154% else
155%     vector=[NaN NaN];
156% end
[225]157
[238]158function [vector,F] = SUBPIX2DGAUSS (result_conv,x,y)
159vector=[0 0]; %default
160F=-2;
161peaky=y;
162peakx=x;
163[npy,npx]=size(result_conv);
164if (x <= npx-1) && (y <= npy-1) && (x >= 1) && (y >= 1)
165    F=0;
[225]166    for i=-1:1
167        for j=-1:1
168            %following 15 lines based on
169            %H. Nobach Æ M. Honkanen (2005)
170            %Two-dimensional Gaussian regression for sub-pixel displacement
171            %estimation in particle image velocimetry or particle position
172            %estimation in particle tracking velocimetry
173            %Experiments in Fluids (2005) 38: 511–515
174            c10(j+2,i+2)=i*log(result_conv(y+j, x+i));
175            c01(j+2,i+2)=j*log(result_conv(y+j, x+i));
176            c11(j+2,i+2)=i*j*log(result_conv(y+j, x+i));
177            c20(j+2,i+2)=(3*i^2-2)*log(result_conv(y+j, x+i));
178            c02(j+2,i+2)=(3*j^2-2)*log(result_conv(y+j, x+i));
179            %c00(j+2,i+2)=(5-3*i^2-3*j^2)*log(result_conv_norm(maxY+j, maxX+i));
180        end
181    end
182    c10=(1/6)*sum(sum(c10));
183    c01=(1/6)*sum(sum(c01));
184    c11=(1/4)*sum(sum(c11));
185    c20=(1/6)*sum(sum(c20));
[231]186    c02=(1/6)*sum(sum(c02));
[225]187    deltax=(c11*c01-2*c10*c02)/(4*c20*c02-c11^2);
188    deltay=(c11*c10-2*c01*c20)/(4*c20*c02-c11^2);
[238]189    if abs(deltax)<1
190        peakx=x+deltax;
191    end
192    if abs(deltay)<1
193        peaky=y+deltay;
194    end
195end
196vector=[peakx-floor(npx/2)-1 peaky-floor(npy/2)-1];
Note: See TracBrowser for help on using the repository browser.