[225] | 1 | % 'pivlab': function piv.m adapted from PIVlab http://pivlab.blogspot.com/ |
---|
| 2 | %-------------------------------------------------------------------------- |
---|
| 3 | % function [xtable ytable utable vtable typevector] = pivlab (image1,image2,interrogationarea, 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 | % interrogationarea: size of the correlation box (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] = pivlab (image1,image2,interrogationarea, step, subpixfinder, mask, roi) |
---|
| 22 | %this funtion performs the DCC PIV analysis. Recent window-deformation |
---|
| 23 | %methods perform better and will maybe be implemented in the future. |
---|
| 24 | warning off %MATLAB:log:logOfZero |
---|
| 25 | if numel(roi)>0 |
---|
| 26 | xroi=roi(1); |
---|
| 27 | yroi=roi(2); |
---|
| 28 | widthroi=roi(3); |
---|
| 29 | heightroi=roi(4); |
---|
| 30 | image1_roi=double(image1(yroi:yroi+heightroi,xroi:xroi+widthroi)); |
---|
| 31 | image2_roi=double(image2(yroi:yroi+heightroi,xroi:xroi+widthroi)); |
---|
| 32 | else |
---|
| 33 | xroi=0; |
---|
| 34 | yroi=0; |
---|
| 35 | image1_roi=double(image1); |
---|
| 36 | image2_roi=double(image2); |
---|
| 37 | end |
---|
| 38 | if numel(mask)>0 |
---|
| 39 | cellmask=mask; |
---|
| 40 | mask=zeros(size(image1_roi)); |
---|
| 41 | for i=1:size(cellmask,1); |
---|
| 42 | masklayerx=cellmask{i,1}; |
---|
| 43 | masklayery=cellmask{i,2}; |
---|
| 44 | mask = mask + poly2mask(masklayerx-xroi,masklayery-yroi,size(image1_roi,1),size(image1_roi,2)); %kleineres eingangsbild und maske geshiftet |
---|
| 45 | end |
---|
| 46 | else |
---|
| 47 | mask=zeros(size(image1_roi)); |
---|
| 48 | end |
---|
| 49 | mask(mask>1)=1; |
---|
| 50 | |
---|
| 51 | miniy=1+(ceil(interrogationarea/2)); |
---|
| 52 | minix=1+(ceil(interrogationarea/2)); |
---|
| 53 | maxiy=step*(floor(size(image1_roi,1)/step))-(interrogationarea-1)+(ceil(interrogationarea/2)); %statt size deltax von ROI nehmen |
---|
| 54 | maxix=step*(floor(size(image1_roi,2)/step))-(interrogationarea-1)+(ceil(interrogationarea/2)); |
---|
| 55 | |
---|
| 56 | numelementsy=floor((maxiy-miniy)/step+1); |
---|
| 57 | numelementsx=floor((maxix-minix)/step+1); |
---|
| 58 | |
---|
| 59 | LAy=miniy; |
---|
| 60 | LAx=minix; |
---|
| 61 | LUy=size(image1_roi,1)-maxiy; |
---|
| 62 | LUx=size(image1_roi,2)-maxix; |
---|
| 63 | shift4centery=round((LUy-LAy)/2); |
---|
| 64 | shift4centerx=round((LUx-LAx)/2); |
---|
| 65 | 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.... |
---|
| 66 | shift4centery=0; |
---|
| 67 | end |
---|
| 68 | 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.... |
---|
| 69 | shift4centerx=0; |
---|
| 70 | end |
---|
| 71 | miniy=miniy+shift4centery; |
---|
| 72 | minix=minix+shift4centerx; |
---|
| 73 | maxix=maxix+shift4centerx; |
---|
| 74 | maxiy=maxiy+shift4centery; |
---|
| 75 | |
---|
| 76 | image1_roi=padarray(image1_roi,[ceil(interrogationarea/2) ceil(interrogationarea/2)], min(min(image1_roi))); |
---|
| 77 | image2_roi=padarray(image2_roi,[ceil(interrogationarea/2) ceil(interrogationarea/2)], min(min(image1_roi))); |
---|
| 78 | mask=padarray(mask,[ceil(interrogationarea/2) ceil(interrogationarea/2)],0); |
---|
| 79 | if (rem(interrogationarea,2) == 0) %for the subpixel displacement measurement |
---|
| 80 | SubPixOffset=1; |
---|
| 81 | else |
---|
| 82 | SubPixOffset=0.5; |
---|
| 83 | end |
---|
| 84 | xtable=zeros(numelementsy,numelementsx); |
---|
| 85 | ytable=xtable; |
---|
| 86 | utable=xtable; |
---|
| 87 | vtable=xtable; |
---|
| 88 | u2table=xtable; |
---|
| 89 | v2table=xtable; |
---|
| 90 | s2n=xtable; |
---|
| 91 | typevector=ones(numelementsy,numelementsx);; |
---|
| 92 | |
---|
| 93 | nrx=0; |
---|
| 94 | nrxreal=0; |
---|
| 95 | nry=0; |
---|
| 96 | increments=0; |
---|
| 97 | |
---|
| 98 | %% MAINLOOP |
---|
| 99 | %handles=guihandles(getappdata(0,'hgui')); |
---|
| 100 | for j = miniy:step:maxiy %vertical loop |
---|
| 101 | nry=nry+1; |
---|
| 102 | % if increments<6 %reduced display refreshing rate |
---|
| 103 | % increments=increments+1; |
---|
| 104 | % else |
---|
| 105 | % increments=1; |
---|
| 106 | % %set(handles.progress, 'string' , ['Frame progress: ' int2str(j/maxiy*100) '%']);drawnow; |
---|
| 107 | % end |
---|
| 108 | n=round((j-miniy)/maxiy*100); |
---|
| 109 | for i = minix:step:maxix % horizontal loop |
---|
| 110 | nrx=nrx+1;%used to determine the pos of the vector in resulting matrix |
---|
| 111 | if nrxreal < numelementsx |
---|
| 112 | nrxreal=nrxreal+1; |
---|
| 113 | else |
---|
| 114 | nrxreal=1; |
---|
| 115 | end |
---|
| 116 | startpoint=[i j]; |
---|
| 117 | image1_crop=image1_roi(j:j+interrogationarea-1, i:i+interrogationarea-1); |
---|
| 118 | image2_crop=image2_roi(ceil(j-interrogationarea/2):ceil(j+1.5*interrogationarea-1), ceil(i-interrogationarea/2):ceil(i+1.5*interrogationarea-1)); |
---|
| 119 | corrmax=0; |
---|
| 120 | if mask(round(j+interrogationarea/2),round(i+interrogationarea/2))==0 |
---|
| 121 | %reference: Oliver Pust, PIV: Direct Cross-Correlation |
---|
| 122 | % image2_crop: sub image with the size of the search area in image 2 |
---|
| 123 | % image1_crop: sub image of the correlation box in image 1 |
---|
| 124 | % %image2_crop is bigger than image1_crop. Zeropading is therefore not |
---|
| 125 | result_conv= conv2(image2_crop,rot90(image1_crop,2),'valid'); |
---|
| 126 | %necessary. 'Valid' makes sure that no zero padded content is returned. |
---|
| 127 | corrmax= max(max(result_conv)); |
---|
| 128 | result_conv=(result_conv/corrmax)*255; %normalize, peak=always 255 |
---|
| 129 | %Find the 255 peak |
---|
| 130 | [y,x] = find(result_conv==255); |
---|
| 131 | if isnan(y)==0 & isnan(x)==0 |
---|
| 132 | try |
---|
| 133 | if subpixfinder==1 |
---|
| 134 | [vector] = SUBPIXGAUSS (result_conv,interrogationarea,x,y,SubPixOffset); |
---|
| 135 | elseif subpixfinder==2 |
---|
| 136 | [vector] = SUBPIX2DGAUSS (result_conv,interrogationarea,x,y,SubPixOffset); |
---|
| 137 | end |
---|
| 138 | catch |
---|
| 139 | vector=[NaN NaN]; %if something goes wrong with cross correlation..... |
---|
| 140 | end |
---|
| 141 | else |
---|
| 142 | vector=[NaN NaN]; %if something goes wrong with cross correlation..... |
---|
| 143 | end |
---|
| 144 | else %if mask was not 0 then |
---|
| 145 | vector=[NaN NaN]; |
---|
| 146 | typevector(nry,nrxreal)=0; |
---|
| 147 | end |
---|
| 148 | |
---|
| 149 | %Create the vector matrix x, y, u, v |
---|
| 150 | xtable(nry,nrxreal)=startpoint(1)+interrogationarea/2; |
---|
| 151 | ytable(nry,:)=startpoint(1,2)+interrogationarea/2; |
---|
| 152 | utable(nry,nrxreal)=vector(1); |
---|
| 153 | vtable(nry,nrxreal)=vector(2); |
---|
| 154 | ctable(nry,nrxreal)=corrmax/sum(sum(image1_crop.*image1_crop)); |
---|
| 155 | end |
---|
| 156 | |
---|
| 157 | end |
---|
| 158 | xtable=xtable-ceil(interrogationarea/2); |
---|
| 159 | ytable=ytable-ceil(interrogationarea/2); |
---|
| 160 | |
---|
| 161 | xtable=xtable+xroi; |
---|
| 162 | ytable=ytable+yroi; |
---|
| 163 | |
---|
| 164 | utable(utable>interrogationarea/1.5)=NaN; |
---|
| 165 | vtable(utable>interrogationarea/1.5)=NaN; |
---|
| 166 | vtable(vtable>interrogationarea/1.5)=NaN; |
---|
| 167 | utable(vtable>interrogationarea/1.5)=NaN; |
---|
| 168 | |
---|
| 169 | function [vector] = SUBPIXGAUSS (result_conv,interrogationarea,x,y,SubPixOffset) |
---|
| 170 | if size(x,1)>1 %if there are more than 1 peaks just take the first |
---|
| 171 | x=x(1:1); |
---|
| 172 | end |
---|
| 173 | if size(y,1)>1 %if there are more than 1 peaks just take the first |
---|
| 174 | y=y(1:1); |
---|
| 175 | end |
---|
| 176 | if (x <= (size(result_conv,1)-1)) && (y <= (size(result_conv,1)-1)) && (x >= 1) && (y >= 1) |
---|
| 177 | %the following 8 lines are copyright (c) 1998, Uri Shavit, Roi Gurka, Alex Liberzon, Technion Israel Institute of Technology |
---|
| 178 | %http://urapiv.wordpress.com |
---|
| 179 | f0 = log(result_conv(y,x)); |
---|
| 180 | f1 = log(result_conv(y-1,x)); |
---|
| 181 | f2 = log(result_conv(y+1,x)); |
---|
| 182 | peaky = y+ (f1-f2)/(2*f1-4*f0+2*f2); |
---|
| 183 | f0 = log(result_conv(y,x)); |
---|
| 184 | f1 = log(result_conv(y,x-1)); |
---|
| 185 | f2 = log(result_conv(y,x+1)); |
---|
| 186 | peakx = x+ (f1-f2)/(2*f1-4*f0+2*f2); |
---|
| 187 | % |
---|
| 188 | SubpixelX=peakx-(interrogationarea/2)-SubPixOffset; |
---|
| 189 | SubpixelY=peaky-(interrogationarea/2)-SubPixOffset; |
---|
| 190 | vector=[SubpixelX, SubpixelY]; |
---|
| 191 | else |
---|
| 192 | vector=[NaN NaN]; |
---|
| 193 | end |
---|
| 194 | |
---|
| 195 | function [vector] = SUBPIX2DGAUSS (result_conv,interrogationarea,x,y,SubPixOffset) |
---|
| 196 | if size(x,1)>1 %if there are more than 1 peaks just take the first |
---|
| 197 | x=x(1:1); |
---|
| 198 | end |
---|
| 199 | if size(y,1)>1 %if there are more than 1 peaks just take the first |
---|
| 200 | y=y(1:1); |
---|
| 201 | end |
---|
| 202 | if (x <= (size(result_conv,1)-1)) && (y <= (size(result_conv,1)-1)) && (x >= 1) && (y >= 1) |
---|
| 203 | for i=-1:1 |
---|
| 204 | for j=-1:1 |
---|
| 205 | %following 15 lines based on |
---|
| 206 | %H. Nobach Æ M. Honkanen (2005) |
---|
| 207 | %Two-dimensional Gaussian regression for sub-pixel displacement |
---|
| 208 | %estimation in particle image velocimetry or particle position |
---|
| 209 | %estimation in particle tracking velocimetry |
---|
| 210 | %Experiments in Fluids (2005) 38: 511515 |
---|
| 211 | c10(j+2,i+2)=i*log(result_conv(y+j, x+i)); |
---|
| 212 | c01(j+2,i+2)=j*log(result_conv(y+j, x+i)); |
---|
| 213 | c11(j+2,i+2)=i*j*log(result_conv(y+j, x+i)); |
---|
| 214 | c20(j+2,i+2)=(3*i^2-2)*log(result_conv(y+j, x+i)); |
---|
| 215 | c02(j+2,i+2)=(3*j^2-2)*log(result_conv(y+j, x+i)); |
---|
| 216 | %c00(j+2,i+2)=(5-3*i^2-3*j^2)*log(result_conv_norm(maxY+j, maxX+i)); |
---|
| 217 | end |
---|
| 218 | end |
---|
| 219 | c10=(1/6)*sum(sum(c10)); |
---|
| 220 | c01=(1/6)*sum(sum(c01)); |
---|
| 221 | c11=(1/4)*sum(sum(c11)); |
---|
| 222 | c20=(1/6)*sum(sum(c20)); |
---|
| 223 | c02=(1/6)*sum(sum(c02)); |
---|
| 224 | %c00=(1/9)*sum(sum(c00)); |
---|
| 225 | |
---|
| 226 | deltax=(c11*c01-2*c10*c02)/(4*c20*c02-c11^2); |
---|
| 227 | deltay=(c11*c10-2*c01*c20)/(4*c20*c02-c11^2); |
---|
| 228 | peakx=x+deltax; |
---|
| 229 | peaky=y+deltay; |
---|
| 230 | |
---|
| 231 | SubpixelX=peakx-(interrogationarea/2)-SubPixOffset; |
---|
| 232 | SubpixelY=peaky-(interrogationarea/2)-SubPixOffset; |
---|
| 233 | vector=[SubpixelX, SubpixelY]; |
---|
| 234 | else |
---|
| 235 | vector=[NaN NaN]; |
---|
| 236 | end |
---|