[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] | 21 | function [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] | 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=[]; |
---|
[231] | 32 | errormsg=''; |
---|
[236] | 33 | %warning off %MATLAB:log:logOfZero |
---|
[231] | 34 | [npy_ima npx_ima]=size(image1); |
---|
[236] | 35 | if ~isequal(size(image2),[npy_ima npx_ima]) |
---|
[231] | 36 | errormsg='image pair with unequal size'; |
---|
| 37 | return |
---|
| 38 | end |
---|
[236] | 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 |
---|
[225] | 58 | end |
---|
[236] | 59 | image1=double(image1); |
---|
| 60 | image2=double(image2); |
---|
[231] | 61 | |
---|
[236] | 62 | %% calculate correlations: MAINLOOP |
---|
| 63 | corrmax=0; |
---|
| 64 | sum_square=1;% default |
---|
[227] | 65 | for 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] | 123 | end |
---|
[236] | 124 | result_conv=result_conv*corrmax/(255*sum_square);% keep the last correlation matrix for output |
---|
[225] | 125 | |
---|
[231] | 126 | |
---|
[238] | 127 | function [vector,F] = SUBPIXGAUSS (result_conv,x,y) |
---|
| 128 | vector=[0 0]; %default |
---|
| 129 | F=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] | 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; |
---|
[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: 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)); |
---|
[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 |
---|
| 195 | end |
---|
| 196 | vector=[peakx-floor(npx/2)-1 peaky-floor(npy/2)-1]; |
---|