Instructor: Lichuan Gui lichuan-gui@uiowa

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

Lecture 26. Sub-pixel Displacement

Correlation Interrogation & FFT Acceleration g(i,j) (m*=10, n*=-5) (m*=-10, n*=5) Evaluation function Auto-correlation function • One double exposed evaluation sample, i.e. g(i,j)=g2(i,j)=g1(i,j) • Displacement determined by positions of the secondary maxima • Two possible velocity directions

Correlation Interrogation & FFT Acceleration (m*=3, n*=5) Evaluation function • Cross-correlation function • Two single exposed evaluation samples, i.e. g2(i,j) & g1(i,j) • Displacement determined by position of the maximum • Velocity direction clear g1(i,j) g2(i,j)

Sub-pixel Displacement Local surface fit of evaluation function Discrete correlation function around point of maximum (m*,n*)

Sub-pixel Displacement Local surface fit of evaluation function Discrete distribution fitted with continuous surface z(x,y)

Sub-pixel Displacement Local surface fit of evaluation function Sub-pixel displacement (x*,y*) determined with position of the maximum of z(x,y)

Sub-pixel Displacement Local surface fit of evaluation function 5-point fit most commonly used

Sub-pixel Displacement Fit function Linear equation system for determining constants Position of the maximum Local surface fit of evaluation function 5-point parabolic fit

Sub-pixel Displacement Fit function Linear equation system for determining constants Position of the maximum Local surface fit of evaluation function 5-point exponential fit (3-point Guassian curve fit)

Arbitrarily sized interrogation window Performance of FFT-based correlation - Effective correlation region A : Effective correlation region B,C,D: Periodically Padded regions

Arbitrarily sized interrogation window Reliability: Some test results with synthetic images  S S f: interrogation window size Performance of FFT-based correlation - Reliability & RMS error dependent on particle image displacement - Reliability & RMS error dependent on interrogation window size - Impossible to determine displacement components not less than half of window side length RMS error: S : measured value S’: real value

Arbitrarily sized interrogation window  Test results with synthetic images (f=32x32) 1  original 2  padded to f=64x64 S S Average padding to next power of 2

More about Correlation Interrogation Fast computation of evaluation function Acceleration with radix-2 based FFT algorithm - Computation time test

More about Correlation Interrogation Effect of linear transformation of evaluation samples on the correlation function - no influence

More about Correlation Interrogation Effect of linear transformation of evaluation samples on the sub-pixel displacement determination with Parabolic surface/curve fit - no influence

More about Correlation Interrogation Effect of linear transformation of evaluation samples on the sub-pixel displacement determination with exponential surface/curve fit - obvious influence - no zero and negative values allowed - normalization suggested before sub-pixel fit

Class project: example of Matlab functions for PIV recording evaluation File name: correlation.m fori=1:M2 % determine coordinates in correlation plane - begin for j=1:N2 m(i,j)=i-M2/2-1; n(i,j)=j-N2/2-1; end end % determine coordinates in correlation plane - end fori=1:M2 % periodical reconstruction - begin for j=1:N2/2 t=C(i,j); C(i,j)=C(i,N2-j+1); C(i,N2-j+1)=t; end end forj=1:N2 for i=1:M2/2 t=C(i,j); C(i,j)=C(M2-i+1,j); C(M2-i+1,j)=t; end end for i=1:M2/2 C1(i,1:N2)=C(i+M2/2,1:N2); end fori=M2/2+1:M2 C1(i,1:N2)=C(i-M2/2,1:N2); end for j=1:N2/2 c(1:M2,j)=C1(1:M2,j+N2/2); end for j=N2/2+1:N2 c(1:M2,j)=C1(1:M2,j-N2/2); end % periodical reconstruction of correlation function - end function[c m n]=correlation(g1,g2) % INPUT PARAMETERS % g1 - 1st evaluation sample, g2 - 2nd evaluation sample % OUTPUT PARAMETERS % c - correlation function, % m - horizontal coordinate, n - vertical coordinate [M N]=size(g1); % determine size of evaluation sample for k=3:12 % determine window size for FFT - begin if 2^k>=M & 2^k>=N break; end end M2=2^k; N2=2^k; % determine window size for FFT - end gm1=mean(mean(g1)); % average gray value padding - begin gm2=mean(mean(g2)); fori=1:M2 for j=1:N2 if i<=M & j<=N gg1(i,j)=g1(i,j); gg2(i,j)=g2(i,j); else gg1(i,j)=gm1; gg2(i,j)=gm2; end end end% average gray value padding - end gfft1=fft2(gg1); % compute correlation function - begin gfft2=conj(fft2(gg2)); C=real(ifft2(gfft1.*gfft2)); % compute correlation function - end

Class project: example of Matlab functions for PIV recording evaluation File name: peaksearch.m function [cm vxvy]=peaksearch(C,m,n,sr,mode,direction) % INPUT PARAMETERS % C - correlation function % m - horizontal coordinate % n - vertical coordinate % sr - search radius % mode - (=0) cross-correlation, (=1) auto-correlation % direction - 1. to right, 2. to left, 3. upwards, 4. downwards % OUTPUT PARAMETERS % cm - high correlation peak value % (vx,vy) – particle image displacement [M N]=size(C); % high peak search - begin ix=0; jy=0; cm=0; fori=M/2+1-sr:M/2+1+sr forj=N/2+1-sr:N/2+1+sr if mode==1 & direction==1 & i<=M/2+1 continue; end if mode==1 & direction==2 & i>=M/2+1 continue; end if mode==1 & direction==3 & j>=N/2+1 continue; end if mode==1 & direction==4 & j<=N/2+1 continue; end max=0; for i1=i-1:i+1 for j1=j-1:j+1 if max<C(i1,j1) max=C(i1,j1); end end end if C(i,j)>=max & C(i,j)>cm ix=i; jy=j; cm=C(i,j); end end end mm=m(ix,jy); nn=n(ix,jy); % high peak search -end % SUB-PIXEL FIT dx=(log(C(ix+1,jy))-log(C(ix-1,jy)))/(4*log(C(ix,jy))-2*log(C(ix+1,jy))-2*log(C(ix-1,jy))); dy=(log(C(ix,jy+1))-log(C(ix,jy-1)))/(4*log(C(ix,jy))-2*log(C(ix,jy+1))-2*log(C(ix,jy-1))); %Displacement vx=mm+dx; vy=nn+dy;

Class project: example of Matlab functions for PIV recording evaluation File name: sample01.m gm=0; % initial value of mean gray value nr=0; % initial value of number of effective pixels fori=1:M for j=1:N ii=i+I-int16(M/2); jj=j+J-int16(N/2); if ii>=1 & ii<nx & jj>=1 & jj<ny%pixels in the image frame Ga=double(G(ii,jj)); Gb=double(G(ii+1,jj)); Gc=double(G(ii,jj+1)); Gd=double(G(ii+1,jj+1)); A=(1-x)*(1-y); 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); nr=nr+1; % count effective pixels else g(i,j)=-1; % pixels out of image frame end end end gm=gm/double(nr); % determine mean gray value fori=1:M for j=1:N if g(i,j)<0 g(i,j)=gm; % fill with mean value end end end function[g]=sample01(G,M,N,X,Y) %INPUT PARAMETERS % G - gray value distribution of the PIV recording % M - interrogation sample width % N - interrogation sample height % X, Y - interrogation samplecoordinates % ---------------------------------------------------- % OUTPUT PARAMETERS % g - gray value distribution of the evaluation sample % ---------------------------------------------------- % I, J – integer part of evaluation sample coordinates % x,y – decimal part of evaluation sample coordinates [nxny]=size(G); I=int16(X); %integer part of coordinates J=int16(Y); x=double(X)-double(I); %decimal part of coordinates y=double(Y)-double(J); if x<0 I=I-1; x=x+1; % x should be positive end ify<0 J=J-1; y=y+1; % y should be positive end

Class project: example of Matlab functions for PIV recording evaluation File name: img2xy.m function [G]=img2xy(A) % A - image % G - gray value distribution [nynx]=size(A); for x=1:nx for y=1:ny G(x,y)=A(ny-y+1,x); end end File name: xy2img.m function [A]=xy2img(G) % A - image % G - gray value distribution [nxny]=size(G); Gmax=max(max(G)); Gmin=min(min(G)); for x=1:nx for y=1:ny A(ny-y+1,x)=int8(double(G(x,y)-Gmin)/double(Gmax-Gmin)*245+5); end end

Class project: example of Matlab functions for PIV recording evaluation File name: main.m clear; A=imread('D001_1.bmp'); % input image file G=img2xy(A); % conver image to gray value distribution % evaluation point x=100; y=200; % interrogation window size M=64; N=64; % evaluation samples g1=sample(G,M,N,x,y); g2=g1; % correlation function [C m n]=correlation(g1,g2); sr=20; % search radius mode=1; % PIV recording mode direction=1; % main flow direction [cm vxvy]=peaksearch(C,m,n,sr,mode,direction) % particle image displacement D=xy2img(C); % convert correlation function to image imshow(D); % display correlation function

Instructor: Lichuan Gui lichuan-gui@uiowa