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=[]; |
---|
23 | function [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. |
---|
26 | nbvec=size(GridIndices,1); |
---|
27 | xtable=zeros(nbvec,1); |
---|
28 | ytable=xtable; |
---|
29 | utable=xtable; |
---|
30 | vtable=xtable; |
---|
31 | ctable=xtable; |
---|
32 | F=xtable; |
---|
33 | result_conv=[]; |
---|
34 | errormsg=''; |
---|
35 | %warning off %MATLAB:log:logOfZero |
---|
36 | [npy_ima npx_ima]=size(image1); |
---|
37 | if ~isequal(size(image2),[npy_ima npx_ima]) |
---|
38 | errormsg='image pair with unequal size'; |
---|
39 | return |
---|
40 | end |
---|
41 | |
---|
42 | %% mask |
---|
43 | testmask=0; |
---|
44 | image1=double(image1); |
---|
45 | image2=double(image2); |
---|
46 | if 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 |
---|
62 | end |
---|
63 | |
---|
64 | %% calculate correlations: MAINLOOP on velocity vectors |
---|
65 | corrmax=0; |
---|
66 | sum_square=1;% default |
---|
67 | for 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; |
---|
125 | end |
---|
126 | result_conv=result_conv*corrmax/(255*sum_square);% keep the last correlation matrix for output |
---|
127 | |
---|
128 | |
---|
129 | function [vector,F] = SUBPIXGAUSS (result_conv,x,y) |
---|
130 | vector=[0 0]; %default |
---|
131 | F=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 | |
---|
160 | function [vector,F] = SUBPIX2DGAUSS (result_conv,x,y) |
---|
161 | vector=[0 0]; %default |
---|
162 | F=-2; |
---|
163 | peaky=y; |
---|
164 | peakx=x; |
---|
165 | [npy,npx]=size(result_conv); |
---|
166 | if (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: 511515 |
---|
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 |
---|
197 | end |
---|
198 | vector=[peakx-floor(npx/2)-1 peaky-floor(npy/2)-1]; |
---|