Instructor: Lichuan Gui lichuan-gui@uiowa

1. Measurements in Fluid Mechanics058:180 (ME:5180)Time & Location: 2:30P - 3:20PMWF 3315 SCOffice Hours: 4:00P – 5:00PMWF223B-5 HL Instructor: Lichuan Gui lichuan-gui@uiowa.edu Phone: 319-384-0594 (Lab), 319-400-5985 (Cell) http://lcgui.net

2. Lecture 35. Particle image identification & tracking

3. Particle Image Identification Purposes - To determined image specifics for particle tracking - To analyze particle size and distribution in flows - To separated phases according to image size - To mask unexpected objects in the flow field - To mask flow boundaries - Others

4. Particle Image Identification Step 1: Preprocess particle images - Increase contrast of the PIV recording - Increase brightness of dark images - Binarize particle images

5. Particle Image Identification Step 2: Identify particle images (method in EDPIV) - Number every image pixel and give a large number to background pixels - Reconstruct number field with minimum in 3×3 neighborhood of image pixels - Iterate until no changes can be made

6. F H  Particle Image Specifics - Particle image position (xs,ys) determined at center of mass - Particle size in pixels (F) - Average brightness of all pixels belong to one image (H) - Particle shape coefficient, e.g. roundness =F/(R2max)

7. Particle Image Specifics - Insect image identification example

8. 1. Original 2. Binary 3. Big particles 4. Expanded 5. Inverted Create Phase Mask

9. - Information from particle image identification Recording pair: 1st frame 2nd frame Particle image number: l1=1,2,···, L1 l2=1,2,···, L2 Size in pixels of particle images F1(l1) F2(l2) Average brightness H1(l1) H2(l2) Shape coefficient 1(l1) 2(l2) Position of particle image x1 (l1),y1 (l1) x2 (l2),y2 (l2) - Tracking function Tracking of Big Particle images Tow-frame tracking of big particles - Weighting coefficients: F +H +=1

10. - Pairing criterion: - Deformation limits: - Particle image displacement: Tracking of Big Particle images Tow-frame tracking of big particles

11. 1st frame 2nd frame Tracking results Tracking of Big Particle images Example

12. Tracking of Insect Motion Example: Tracking fire ant in successive image frames Fire ant trace Fire ant speed time history

13. B One pixel for one particle C Maxima search Erosion B B C C Combine Position of LID particle images

14. Evaluation with large window Evaluation at identified particles Evaluation of LID recordings • Individual particle image pattern tracking LID recordings LID recordings with small interrogation window

15. Simulated PIV recording pair of uniform flow Particle image diameter: 24 pixels Particle image displacement: Sx=2.3, Sy=1.5 pixels Particle tracking results RMS error=0.25 pixels Pattern tracking results 4x4-pixel pattern RMS error=0.06 pixels Tracking of small particle images - Nearest displacement criterion (distance between particles >> displacement) - Pairing with help of neighbor particle images (e.g. Okamoto et al. 1995) - Multi-frame tracking (e.g. Hassan & Canaan 1991 ) - Displacement determined according to particle image center positions - Individual particle image pattern tracking recommended for LID PIV

16. References • GuiL, Merzkirch W, Shu JZ (1997) Evaluation of low image density PIV recordings with the MQD method and application to the flow in a liquid bridge. Journal of Flow Visualization and Image Processing, vol. 4: 333-343  • Gui L and Seiner JM (2009) Application of an Image Tracking Algorithm in Fire Ant Motion Experiment. Algorithms 2, no. 2:735-749

17. Matlab function for 4-P CDIC 7 8 9 4 5 6 1 2 3 function[g]=sample4P(G,M,N,Xm,Ym,Sx,Sy,C) %INPUT PARAMETERS % G - gray value distribution of the PIV recording % M - interrogation sample width % N - interrogation sample height % Xm,Ym - interrogation sample location % Sx,Sy - displacements at 9 points % C=-1 for f1(i,j), C=1 for f2(i,j) % OUTPUT PARAMETERS % g - gray value distribution of the evaluation sample [nxny]=size(G); % image size Xws=Sx(5); % window shift Yws=Sy(5); Xdis=Sx-(Sx(1)+Sx(3)+Sx(7)+Sx(9))/4; % distortion function Ydis=Sy-(Sy(1)+Sy(3)+Sy(7)+Sy(9))/4; % at 9 points Xpix=C*(Xws+Xdis)/2; % pixel displacement Ypix=C*(Yws+Ydis)/2; % at 9 points - Particle image sisplacements at center and 4 corners (i.e. S1,S3,S5,S7,S9) determined according to a previus evaluation - Window shift determined with displacement in the window center, i.e. Sws=S5 - Image distortion at the 4 points determined as C=+1: C=-1:

18. Matlab function for 4-P CDIC gm=0; % initial average gray value nr=0; % initial number of effective pixels for i=1:M % column loop start for j=1:N % row loop start A=(M-i)*(N-j)/double((M-1)*(N-1)); % weighting coefficient for point 1 B=(i-1)*(N-j)/double((M-1)*(N-1)); % weighting coefficient for point 3 C=(M-i)*(j-1)/double((M-1)*(N-1)); % weighting coefficient for point 7 D=(i-1)*(j-1)/double((M-1)*(N-1)); % weighting coefficient for point 9 x_pix=Xpix(1)*A+Xpix(3)*B+Xpix(7)*C+Xpix(9)*D; % pixel displacement at current pixel y_pix=Ypix(1)*A+Ypix(3)*B+Ypix(7)*C+Ypix(9)*D; % pixel displacement at current pixel X=Xm+x_pix-M/2+i; % corresponding x position of current pixel in the PIV recording Y=Ym+y_pix-N/2+j; % corresponding y position of current pixel in the PIV recording I=int16(X); % integer portion of x-position J=int16(Y); % integer portion of y-position x=double(X)-double(I); % decimal portion of x-position y=double(Y)-double(J); % decimal portion of y-position if x<0 % adjust values so that x≥0, y≥0 I=I-1; x=x+1; end if y<0 J=J-1; y=y+1; end j=N B A D C j=1 i=1 i=M 1 3 7 9

19. Matlab function for 4-P CDIC if I>=1 & I<nx & J>=1 & J<ny% limited in the image frame Ga=double(G(I,J)); % gray value at integer pixels Gb=double(G(I+1,J)); Gc=double(G(I,J+1)); Gd=double(G(I+1,J+1)); A=(1-x)*(1-y); % weighting coefficients for interpolation B=x*(1-y); C=(1-x)*y; D=x*y; g(i,j)=A*Ga+B*Gb+C*Gc+D*Gd; % bilinear interpolation gm=gm+g(i,j); % sum of gray values for averaging nr=nr+1; % count number of effective pixels else g(i,j)=-1; % temporary value for pixel out of image frame end end % row loop end end % column loop end gm=gm/double(nr); % average gray value of effective pixels for i=1:M for j=1:N if g(i,j)<0 g(i,j)=gm; % fill with average value for pixel out of image frame end end end J+1 B A D C J I I+1 Gd Gc Ga Gb

20. Matlab program for 4-P CDIC clear; % clear variables A1=imread('A001_1.bmp'); % input 1st image in the recording pair A2=imread('A001_2.bmp'); % input 2nd image file G1=img2xy(A1); % convert image to gray value distribution G2=img2xy(A2); % convert image to gray value distribution Mg=16; % interrogation grid width Ng=16; % interrogation grid height M=2*Mg; % interrogation window width w. 50% overlap N=2*Ng; % interrogation window height w. 50% overlap sr1=12; % initial search radius sr2=6; % final search radius NN=6; % iteration number dU=[-12 12 3]; % parameters for error detection dV=[-12 12 3]; % parameters for error detection [nxny]=size(G1); % determine size of the image col=400/Mg; % number of grid rows in limited area of 400-pixel in height row=400/Ng; % number of grid columns in limited area of 400-pixel in width

21. Matlab program for 4-P CDIC for i=1:col for j=1:row X(i,j)=double((i-1)*Mg+400); % x-position of interrogation point Y(i,j)=double((j-1)*Ng+300); % y-position of interrogation point U(i,j)=double(0); % initial particle image displacement in x-direction V(i,j)=double(0);% initial particle image displacement in y-direction end end for nn=1:NN % iteration begin sr=int16((nn-1)*(sr2-sr1)/(NN-1)+sr1); % determine search radius if nn>1 [U V valid]=interpolation(U,V, valid); % interpolation for at wrong vectors [U V valid]=interpolation(U,V, valid); % second pass of interpolation end % iteration may be necessary in complicated case

22. Matlab program for 4-P CDIC for i=1:col % column loop start for j=1:row % row loop start if nn==1 wsx=0;% set window shift to 0 in the first run wsy=0; else if valid(i,j)>0 wsx=U(i,j); % window shift determined with previous evaluation wsy=V(i,j); end end nr=0; % determining particle image displacement at 9 points in the window begin for q=-1:1 for p=-1:1 nr=nr+1; % number of grid point in the window if i>1 & i<col & j>1 & j<row & nn>1 % after the first run & when all the 9 pints have valid vectors sx(nr)=U(i+p,j+q); % determine displacements at 9 points in the window sy(nr)=V(i+p,j+q); % with results of previous evaluation else sx(nr)=wsx; % ignore image distortion sy(nr)=wsy; end end end % determining particle image displacement at 9 points in the window end

23. Matlab program for 4-P CDIC x=X(i,j); % determine horizontal coordinate of interrogation point y=Y(i,j); % determine vertical coordinate of interrogation point g1=sample4P(G1,M,N,x,y, sx, sy, -1); % evaluation sample with backward image correction g2=sample4P(G2,M,N,x,y, sx, sy, 1); % evaluation sample with forward image correction [C m n]=correlation(g1,g2); % calculating correlation function [cm vxvy]=peaksearch(C,m,n,sr,0,0); % determine particle image displacement U(i,j)=vx+wsx; % adjust particle image displacement with window shift V(i,j)=vy+wsy; % adjust particle image displacement with window shift end % row loop end end % column loop end valid=errordetection(U,V,dU,dV); % detect evaluation errors end % iteration end quiver(X,Y,U,V); % plot vector map