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 typevector 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 | errormsg=''; |
---|
25 | warning off %MATLAB:log:logOfZero |
---|
26 | [npy_ima npx_ima]=size(image1); |
---|
27 | if ~isequal(size(image1),size(image2)) |
---|
28 | errormsg='image pair with unequal size'; |
---|
29 | return |
---|
30 | end |
---|
31 | xroi=0; |
---|
32 | yroi=0; |
---|
33 | image1_roi=double(image1); |
---|
34 | image2_roi=double(image2); |
---|
35 | if numel(mask)>0 |
---|
36 | cellmask=mask; |
---|
37 | mask=zeros(size(image1)); |
---|
38 | for i=1:size(cellmask,1); |
---|
39 | masklayerx=cellmask{i,1}; |
---|
40 | masklayery=cellmask{i,2}; |
---|
41 | mask = mask + poly2mask(masklayerx-xroi,masklayery-yroi,npy_ima,npx_ima); %kleineres eingangsbild und maske geshiftet |
---|
42 | end |
---|
43 | else |
---|
44 | mask=zeros(size(image1)); |
---|
45 | end |
---|
46 | mask(mask>1)=1; |
---|
47 | |
---|
48 | % ibx=2*ibx2-1;%ibx and iby odd, reduced by 1 if even |
---|
49 | % iby=2*iby2-1; |
---|
50 | % miniy=1+iby2 |
---|
51 | % minix=1+ibx2 |
---|
52 | % maxiy=step*(floor(size(image1_roi,1)/step))-(iby-1)+iby2 %statt size deltax von ROI nehmen |
---|
53 | % maxix=step*(floor(size(image1_roi,2)/step))-(ibx-1)+ibx2 |
---|
54 | % numelementsy=floor((maxiy-miniy)/step+1); |
---|
55 | % numelementsx=floor((maxix-minix)/step+1); |
---|
56 | % |
---|
57 | % LAy=miniy; |
---|
58 | % LAx=minix; |
---|
59 | % LUy=size(image1_roi,1)-maxiy; |
---|
60 | % LUx=size(image1_roi,2)-maxix; |
---|
61 | % shift4centery=round((LUy-LAy)/2); |
---|
62 | % shift4centerx=round((LUx-LAx)/2); |
---|
63 | % if shift4centery<0 %shift4center will be negative if in the unshifted case the left border is bigger than the right border. the vectormatrix is hence not centered on the image. the matrix cannot be shifted more towards the left border because then image2_crop would have a negative index. The only way to center the matrix would be to remove a column of vectors on the right side. but then we weould have less data.... |
---|
64 | % shift4centery=0; |
---|
65 | % end |
---|
66 | % if shift4centerx<0 %shift4center will be negative if in the unshifted case the left border is bigger than the right border. the vectormatrix is hence not centered on the image. the matrix cannot be shifted more towards the left border because then image2_crop would have a negative index. The only way to center the matrix would be to remove a column of vectors on the right side. but then we weould have less data.... |
---|
67 | % shift4centerx=0; |
---|
68 | % end |
---|
69 | % miniy=miniy+shift4centery; |
---|
70 | % minix=minix+shift4centerx; |
---|
71 | % maxix=maxix+shift4centerx; |
---|
72 | % maxiy=maxiy+shift4centery; |
---|
73 | |
---|
74 | image1_roi=padarray(image1_roi,[iby2 ibx2], min(min(image1_roi)));%add a border around the image with minimum image value |
---|
75 | image2_roi=padarray(image2_roi,[iby2 ibx2], min(min(image1_roi))); |
---|
76 | mask=padarray(mask,[iby2 ibx2],0); |
---|
77 | %SubPixOffset=0.5;%odd values chosen for ibx and iby |
---|
78 | |
---|
79 | nbvec=size(GridIndices,1); |
---|
80 | xtable=zeros(nbvec,1); |
---|
81 | ytable=xtable; |
---|
82 | utable=xtable; |
---|
83 | vtable=xtable; |
---|
84 | u2table=xtable; |
---|
85 | v2table=xtable; |
---|
86 | s2n=xtable; |
---|
87 | typevector=ones(size(xtable)); |
---|
88 | |
---|
89 | nrx=0; |
---|
90 | nrxreal=0; |
---|
91 | nry=0; |
---|
92 | increments=0; |
---|
93 | |
---|
94 | %% MAINLOOP |
---|
95 | for ivec=1:nbvec |
---|
96 | iref=GridIndices(ivec,1); |
---|
97 | jref=GridIndices(ivec,2); |
---|
98 | %jref=npy_ima-PointCoord(ivec,2)+1; |
---|
99 | image1_crop=image1_roi(jref-iby2:jref+iby2,iref-ibx2:iref+ibx2); |
---|
100 | image2_crop=image2_roi(jref+shifty-isy2:jref+shifty+isy2,iref+shiftx-isx2:iref+shiftx+isx2); |
---|
101 | if mask(jref,iref)==0 |
---|
102 | %reference: Oliver Pust, PIV: Direct Cross-Correlation |
---|
103 | % image2_crop: sub image with the size of the search area in image 2 |
---|
104 | % image1_crop: sub image of the correlation box in image 1 |
---|
105 | % %image2_crop is bigger than image1_crop. Zeropading is therefore not |
---|
106 | result_conv= conv2(image2_crop,rot90(image1_crop,2),'valid'); |
---|
107 | %necessary. 'Valid' makes sure that no zero padded content is returned. |
---|
108 | corrmax= max(max(result_conv)); |
---|
109 | result_conv=(result_conv/corrmax)*255; %normalize, peak=always 255 |
---|
110 | % result_conv=flipdim(result_conv,2);%reverse x direction |
---|
111 | %Find the 255 peak |
---|
112 | [y,x] = find(result_conv==255); |
---|
113 | if isnan(y)==0 & isnan(x)==0 |
---|
114 | try |
---|
115 | if subpixfinder==1 |
---|
116 | [vector] = SUBPIXGAUSS (result_conv,x,y); |
---|
117 | elseif subpixfinder==2 |
---|
118 | [vector] = SUBPIX2DGAUSS (result_conv,x,y); |
---|
119 | end |
---|
120 | catch ME |
---|
121 | errormsg=ME.message |
---|
122 | vector=[0 0]; %if something goes wrong with cross correlation..... |
---|
123 | end |
---|
124 | else |
---|
125 | vector=[0 0]; %if something goes wrong with cross correlation..... |
---|
126 | end |
---|
127 | else %if mask was not 0 then |
---|
128 | vector=[0 0]; |
---|
129 | typevector(ivec)=0; |
---|
130 | end |
---|
131 | |
---|
132 | %Create the vector matrix x, y, u, v |
---|
133 | xtable(ivec)=GridIndices(ivec,1); |
---|
134 | ytable(ivec)=GridIndices(ivec,2); |
---|
135 | utable(ivec)=vector(1); |
---|
136 | vtable(ivec)=vector(2); |
---|
137 | sum_square=sum(sum(image1_crop.*image1_crop)); |
---|
138 | ctable(ivec)=corrmax/sum_square; |
---|
139 | end |
---|
140 | result_conv=result_conv/255; |
---|
141 | |
---|
142 | |
---|
143 | function [vector] = SUBPIXGAUSS (result_conv,x,y) |
---|
144 | |
---|
145 | if size(x,1)>1 %if there are more than 1 peaks just take the first |
---|
146 | x=x(1:1); |
---|
147 | end |
---|
148 | if size(y,1)>1 %if there are more than 1 peaks just take the first |
---|
149 | y=y(1:1); |
---|
150 | end |
---|
151 | if (x <= (size(result_conv,1)-1)) && (y <= (size(result_conv,1)-1)) && (x >= 1) && (y >= 1) |
---|
152 | %the following 8 lines are copyright (c) 1998, Uri Shavit, Roi Gurka, Alex Liberzon, Technion Israel Institute of Technology |
---|
153 | %http://urapiv.wordpress.com |
---|
154 | f0 = log(result_conv(y,x)); |
---|
155 | f1 = log(result_conv(y-1,x)); |
---|
156 | f2 = log(result_conv(y+1,x)); |
---|
157 | peaky = y+ (f1-f2)/(2*f1-4*f0+2*f2); |
---|
158 | f0 = log(result_conv(y,x)); |
---|
159 | f1 = log(result_conv(y,x-1)); |
---|
160 | f2 = log(result_conv(y,x+1)); |
---|
161 | peakx = x+ (f1-f2)/(2*f1-4*f0+2*f2); |
---|
162 | [npy,npx]=size(result_conv); |
---|
163 | vector=[peakx-floor(npx/2)-1 peaky-floor(npy/2)-1]; |
---|
164 | else |
---|
165 | vector=[NaN NaN]; |
---|
166 | end |
---|
167 | |
---|
168 | function [vector] = SUBPIX2DGAUSS (result_conv,x,y) |
---|
169 | if size(x,1)>1 %if there are more than 1 peaks just take the first |
---|
170 | x=x(1:1); |
---|
171 | end |
---|
172 | if size(y,1)>1 %if there are more than 1 peaks just take the first |
---|
173 | y=y(1:1); |
---|
174 | end |
---|
175 | if (x <= (size(result_conv,1)-1)) && (y <= (size(result_conv,1)-1)) && (x >= 1) && (y >= 1) |
---|
176 | for i=-1:1 |
---|
177 | for j=-1:1 |
---|
178 | %following 15 lines based on |
---|
179 | %H. Nobach Æ M. Honkanen (2005) |
---|
180 | %Two-dimensional Gaussian regression for sub-pixel displacement |
---|
181 | %estimation in particle image velocimetry or particle position |
---|
182 | %estimation in particle tracking velocimetry |
---|
183 | %Experiments in Fluids (2005) 38: 511515 |
---|
184 | c10(j+2,i+2)=i*log(result_conv(y+j, x+i)); |
---|
185 | c01(j+2,i+2)=j*log(result_conv(y+j, x+i)); |
---|
186 | c11(j+2,i+2)=i*j*log(result_conv(y+j, x+i)); |
---|
187 | c20(j+2,i+2)=(3*i^2-2)*log(result_conv(y+j, x+i)); |
---|
188 | c02(j+2,i+2)=(3*j^2-2)*log(result_conv(y+j, x+i)); |
---|
189 | %c00(j+2,i+2)=(5-3*i^2-3*j^2)*log(result_conv_norm(maxY+j, maxX+i)); |
---|
190 | end |
---|
191 | end |
---|
192 | c10=(1/6)*sum(sum(c10)); |
---|
193 | c01=(1/6)*sum(sum(c01)); |
---|
194 | c11=(1/4)*sum(sum(c11)); |
---|
195 | c20=(1/6)*sum(sum(c20)); |
---|
196 | c02=(1/6)*sum(sum(c02)); |
---|
197 | deltax=(c11*c01-2*c10*c02)/(4*c20*c02-c11^2); |
---|
198 | deltay=(c11*c10-2*c01*c20)/(4*c20*c02-c11^2); |
---|
199 | peakx=x+deltax; |
---|
200 | peaky=y+deltay; |
---|
201 | |
---|
202 | [npy,npx]=size(result_conv); |
---|
203 | vector=[peakx-floor(npx/2)-1 peaky-floor(npy/2)-1]; |
---|
204 | % SubpixelX=peakx-(ibx/2)-SubPixOffset; |
---|
205 | % SubpixelY=peaky-(iby/2)-SubPixOffset; |
---|
206 | % vector=[SubpixelX, SubpixelY]; |
---|
207 | else |
---|
208 | vector=[NaN NaN]; |
---|
209 | end |
---|