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 | % ibx,iby: size of the correlation box along x and y (in px) |
---|
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=[]; |
---|
21 | function [xtable ytable utable vtable ctable F result_conv errormsg] = pivlab (image1,image2,ibx2,iby2,isx2,isy2,shiftx,shifty, GridIndices, subpixfinder,mask) |
---|
22 | %this funtion performs the DCC PIV analysis. Recent window-deformation |
---|
23 | %methods perform better and will maybe be implemented in the future. |
---|
24 | nbvec=size(GridIndices,1); |
---|
25 | xtable=zeros(nbvec,1); |
---|
26 | ytable=xtable; |
---|
27 | utable=xtable; |
---|
28 | vtable=xtable; |
---|
29 | ctable=xtable; |
---|
30 | F=xtable; |
---|
31 | result_conv=[]; |
---|
32 | errormsg=''; |
---|
33 | %warning off %MATLAB:log:logOfZero |
---|
34 | [npy_ima npx_ima]=size(image1); |
---|
35 | if ~isequal(size(image2),[npy_ima npx_ima]) |
---|
36 | errormsg='image pair with unequal size'; |
---|
37 | return |
---|
38 | end |
---|
39 | |
---|
40 | %% mask |
---|
41 | testmask=0; |
---|
42 | if 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 |
---|
58 | end |
---|
59 | image1=double(image1); |
---|
60 | image2=double(image2); |
---|
61 | |
---|
62 | %% calculate correlations: MAINLOOP |
---|
63 | corrmax=0; |
---|
64 | sum_square=1;% default |
---|
65 | for ivec=1:nbvec |
---|
66 | iref=GridIndices(ivec,1); |
---|
67 | jref=GridIndices(ivec,2); |
---|
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 |
---|
92 | [y,x] = find(result_conv==255,1); |
---|
93 | if ~isnan(y) && ~isnan(x) |
---|
94 | try |
---|
95 | if subpixfinder==1 |
---|
96 | [vector,F(ivec)] = SUBPIXGAUSS (result_conv,x,y); |
---|
97 | elseif subpixfinder==2 |
---|
98 | [vector,F(ivec)] = SUBPIX2DGAUSS (result_conv,x,y); |
---|
99 | end |
---|
100 | sum_square=sum(sum(image1_crop.*image1_crop)); |
---|
101 | ctable(ivec)=corrmax/sum_square;% correlation value |
---|
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 |
---|
105 | catch ME |
---|
106 | vector=[0 0]; %if something goes wrong with cross correlation..... |
---|
107 | F(ivec)=3; |
---|
108 | end |
---|
109 | else |
---|
110 | vector=[0 0]; %if something goes wrong with cross correlation..... |
---|
111 | F(ivec)=3; |
---|
112 | end |
---|
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); |
---|
123 | end |
---|
124 | result_conv=result_conv*corrmax/(255*sum_square);% keep the last correlation matrix for output |
---|
125 | |
---|
126 | |
---|
127 | function [vector,F] = SUBPIXGAUSS (result_conv,x,y) |
---|
128 | vector=[0 0]; %default |
---|
129 | F=0; |
---|
130 | [npy,npx]=size(result_conv); |
---|
131 | |
---|
132 | % if (x <= (size(result_conv,1)-1)) && (y <= (size(result_conv,1)-1)) && (x >= 1) && (y >= 1) |
---|
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 |
---|
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 |
---|
153 | vector=[peakx-floor(npx/2)-1 peaky-floor(npy/2)-1]; |
---|
154 | % else |
---|
155 | % vector=[NaN NaN]; |
---|
156 | % end |
---|
157 | |
---|
158 | function [vector,F] = SUBPIX2DGAUSS (result_conv,x,y) |
---|
159 | vector=[0 0]; %default |
---|
160 | F=-2; |
---|
161 | peaky=y; |
---|
162 | peakx=x; |
---|
163 | [npy,npx]=size(result_conv); |
---|
164 | if (x <= npx-1) && (y <= npy-1) && (x >= 1) && (y >= 1) |
---|
165 | F=0; |
---|
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: 511515 |
---|
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)); |
---|
186 | c02=(1/6)*sum(sum(c02)); |
---|
187 | deltax=(c11*c01-2*c10*c02)/(4*c20*c02-c11^2); |
---|
188 | deltay=(c11*c10-2*c01*c20)/(4*c20*c02-c11^2); |
---|
189 | if abs(deltax)<1 |
---|
190 | peakx=x+deltax; |
---|
191 | end |
---|
192 | if abs(deltay)<1 |
---|
193 | peaky=y+deltay; |
---|
194 | end |
---|
195 | end |
---|
196 | vector=[peakx-floor(npx/2)-1 peaky-floor(npy/2)-1]; |
---|