source: trunk/src/pivlab.m @ 297

Last change on this file since 297 was 273, checked in by sommeria, 13 years ago

GUI civ renovated with panels

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