Experiments in Fluids [Suppl.] S3±S12 Ó Springer-Verlag 2000

Theoretical analysis of the measurement precision in particle image velocimetry J. Westerweel

Abstract An analytical expression for the error in the subpixel displacement is derived. The expression relates the error to the second- and fourth-order moments of the intensity probability density function, which can both be expressed in terms of the source density. Scaling relations for the error as a function of the particle-image diameter are determined by means of the CrameÁr±Rao lower bound, for both constant image density and constant source density. The analytical expression is evaluated for the error as a function of the sub-pixel displacement, the particleimage diameter, and the image density, and the results are validated against results from Monte-Carlo simulations. It is demonstrated that results obtained from synthetic PIV images that do not represent the correct image intensity probability density function systematically underestimate the displacement measurement error.

1 Introduction Over the past decade PIV has undergone substantial improvement in measurement precision and reliability. Each new interrogation algorithm challenges the numerical precision and reliability of preceding ones. Initially, it was thought that the measurement uncertainty would be given by the pixel dimension, i.e., 0.5 pixel units (Hesselink 1988). As it turned out, this estimate for the precision was too conservative. With the introduction of sub-pixel interpolation, the measurement precision was improved substantially, and it became possible to obtain measurements with an uncertainty level of about 0.1 pixel units (Willert and Gharib 1991). The measurement reliability and precision were further improved by means of using windows with different size (Keane and Adrian 1992), and by using offset interrogation windows, by which the measurement error could be reduced to 0.04 pixel units (under controlled circumstances), as was demonstrated in an actual experiment (Westerweel et al. 1997). Even newer algorithms have been proposed that further challenge the numerical precision and reliability of the preceding ones. This raises the question of to what level the measurement precision can ultimately be reduced. The precision of a given algorithm is usually determined by means of Monte-Carlo simulations. PIV images with a known displacement are generated synthetically, J. Westerweel Laboratory for Aero and Hydrodynamics Delft University of Technology Rotterdamseweg 145, 2628 AL Delft The Netherlands

and the measurement precision and reliability are determined by comparing the analysis results against the actual displacement. Extensive parameter studies were carried out by Keane and Adrian (1992) and Willert (1996). For such simulation studies, it is very important that the synthetic PIV images are realistic representations of actual PIV images. However, the simulation studies only describe but do not explain the behavior of the measurement precision in relation to the experimental parameters. Another approach of investigating the measurement precision is by means of a theoretical investigation, where the aim is to obtain an analytical expression for the measurement error. The groundwork for this theory was laid down by Adrian (1988), and was later adopted for digital PIV images (Westerweel 1993). This paper describes recent improvements and additions to this theoretical description. For example, the initial theory for digital images did not take into account the proper nature of the image statistics for PIV images, so that empirical constants were required to ®nd agreement with observations. Also, the theory initially only described the measurement error for exact integer displacements, but has now been extended to include also any fractional displacement. The structure of this paper is as follows. In Sect. 2 an expression is derived for the measurement error of the displacement at sub-pixel level. The expression includes parameters that depend on the image statistics, which will be discussed in Sect. 3. Before the expression is evaluated (Sect. 5), the scaling for the error is discussed in Sect. 4. The paper concludes with a discussion in Sect. 6 and a summary of the main conclusions in Sect. 7.

2 PIV Interrogation Consider a single-exposure frame-pair. It is assumed that both frames were recorded with identical illumination conditions. The pixel gray value / at coordinates [m, n] is given by: /‰m; nŠ ˆ_ cdr2 I…Xm ; Yn † ‡ mm;n …1† where Xm ˆ mdr, Yn ˆ ndr, I(X, Y) is the exposure per unit area on the image sensor, c is a conversion parameter between the image exposure and the pixel gray value, dr is the pixel dimension, and mm,n is a random noise term that accounts for sensor noise and discretization noise. (The pixel area dr2 could have been absorbed in c, but this particular form was chosen so that the dimension of c is gray value per unit exposure.) The expression in Eq. (1) assumes a linear sensor response, which is generally satis®ed for a charge coupled device (CCD), and it would be

S3

valid when the pixel size is small with respect to the particle-image diameter; to account for under-resolved sampling I(X,Y)would need to be replaced by ~I…X; Y†, which is the convolution of I(X,Y) with the spatial sensitivity of the pixel (Westerweel 1993, 1997). The image frames are interrogated by the computation of the discrete spatial cross-correlation in an M ´ N-pixel interrogation region

R‰s; tŠ ˆ S4

M X N 1 X W1 ‰m; nŠ/1 ‰m; nŠ MN mˆ1 nˆ1

PIV), and the location of this peak corresponds to the displacement of the particle images. The ¯uctuation of RD is the origin of a random error between the measured and actual value of the displacement. In addition to the random error, the measured displacement usually includes also a bias error. Sub-pixel interpolation estimators, such as the peak centroid and the Gaussian peak-®t in Eq. (5), have a bias error towards an integer pixel value of the displacement. This leads to an effect called pixel locking, which typically occurs when the particle images become small with respect to the pixel dimension, e.g., see Prasad et al. (1992). Peak locking effects are supposed to occur only when the particle-image diameter is less than two pixel units (Westerweel 2000). A more subtle bias error occurs in relation to the ®nite dimensions of W1 and W2, which is the cause of a slight skew of the displacementcorrelation peak. The result is a small negative bias of the measured displacement magnitude; to remove this bias the computed correlation values R±1, R0 and R+1 must be divided by the local value of FI, i.e., the in-plane loss-ofcorrelation (Westerweel 1997). The bias vanishes implicitly when the correlation is computed with offset interrogation windows or with a window W2 that is large enough to contain all particle images that occur in W1, or vice versa. In the subsequent analysis, it is assumed that m0 ˆ n0 ˆ 0. This can always be accomplished by a proper shift of /1 with respect to /2 (Westerweel et al. 1997). Also, it is assumed that the ¯uctuations in ^eX and ^eY are orthogonal, i.e., cov…^eX ; ^eY † ˆ 0, so the estimation of the sub-pixel displacement can be reduced to a one-dimensional problem. This assumption is valid for a (nearly) uniform displacement ®eld. Finally, to keep the mathematics compact, only square interrogation windows (i.e., M ˆ N) of equal size are considered. The variance of ^e as a result of the ¯uctuations in RD is given by the error propagation equation

 W2 ‰m ‡ s; n ‡ tŠ/2 ‰m ‡ s; n ‡ tŠ …2† where /1 and /2 are the ®rst and second image frames respectively, and W1 and W2 the interrogation windows in /1 and /2 respectively. When the image intensity is split into mean and ¯uctuating parts, R[s, t] can be written as (Keane and Adrian 1992): R‰s; tŠ ˆ RC ‰s; tŠ ‡ RF ‰s; tŠ ‡ RD ‰s; tŠ ; …3† where RC is the correlation of the mean image intensity, RF the correlation of the mean with the ¯uctuating image intensity, and RD the correlation of the intensity ¯uctuations. The terms RC and RF vanish by subtracting the mean intensity from the image ®elds. The particle-image displacement is determined by measuring the location of the tallest correlation peak in RD. For digital PIV images, the displacement is usually split into integer and fractional pixel units (px), i.e., DX ˆ …m0 ‡ eX †dr and DY ˆ …n0 ‡ eY †dr ; …4† where [m0, n0] is the location of the correlation maximum, and dr the linear pixel dimension. The fractional part of the displacement is determined by interpolation of the correlation peak. The most common interpolation method is the so-called three-point Gaussian ®t (Willert and Gharib 1991), where eX is estimated by ln R 1 ln R‡1 ^eX ˆ ; …5† X X o^e o^e 2…ln R 1 ‡ ln R‡1 2 ln R0 †
r2DX  dr2
cov Ri ; Rj : with oRi oRj i j R

 RD ‰m0 1; n0 Š; R0  RD ‰m0 ; n0 Š; …6† R‡1  RD ‰m0 ‡ 1; n0 Š ; and vice versa for eY. The correlation RD can be split into mean and ¯uctuating parts; see Fig. 1. The mean correlation consist of a single correlation peak (for single-exposure/double-frame 1

Fig. 1. The cross-correlation (left) of the image intensity ¯uctuations is split into a mean part (middle) and a ¯uctuating part (right)

…7†

This expression includes nine terms, but only six are unique. The products of the partial-derivative terms for the Gaussian peak ®t are plotted in Fig. 2. This shows that the terms that do not contain the zero index are dominant. The covariance for the spatial cross correlation is given by (Priestley 1992; Westerweel 1993):

3 Image statistics The statistical properties of PIV images are characterized by the source density NS, de®ned as (Adrian 1984):
p NS ˆ CDz0 =M02 ds2 …10† 4 where C is the number density of the tracer particles, Dz0 the light sheet width, M0 the image magni®cation, and ds the mean particle-image diameter. The mean pixel gray value is given by (for a zero mean noise level):
: h/i ˆ cdr2 hIi; with: hIi ˆ s~00 CDz0 =M02 ˆ s~00
NS = p4 ds2 ; …11†

Fig. 2. The products of the partial-derivatives in Eq. (7) for the three-point Gaussian ®t as a function of the fractional displacement

where s^00 represents the mean exposure of a single particle image (Adrian 1988).1 The conditional ensemble mean of RD[s, t] for a given (near-uniform) velocity u is given by (Keane and Adrian 1992):

hRD ‰u; vŠjui ˆ_ cdr4 s^200 NI FI FO Fs

…12†

with su ˆ udr, tv ˆ vdr, where NI is the image density, FI and FO the in-plane and out-of-plane loss-of-correlation respectively (due to the ®nite size of the interrogation domain and of the light sheet depth), and Fs is the particle …8† image self-correlation. The expression is valid when the particle-image diameter is large with respect to the pixel ‡ hRD ‰m; nŠjuihRD ‰m ‡ s ‡ u; n ‡ t ‡ vŠjui dimension; for under-resolved images Fs has to be replaced with F~s , which is the convolution of Fs with the self‡K4 ‰m; n; s; t; u; vŠg correlation of the spatial pixel sensitivity (Westerweel where hRD ‰m; nŠjui is the mean correlation over the 1993, 1997). In Sect. 2 it was assumed that m0 ˆ n0 ˆ 0, ensemble of all possible realizations of the images for a which implies: FI ˆ 1 (Keane and Adrian 1992). Furthergiven velocity u; hRP ‰m; nŠi ˆ hRD ‰m; nŠju ˆ 0i the selfmore, it is assumed that the out-of-plane motion is small, correlation of the intensity ¯uctuations, and K4 a fourth- so that FO  1. The image density is de®ned as (Adrian order correlation term. If the intensity ¯uctuations have a 1984):
2 joint-normal pdf, then the last term in Eq. (8) can be 2 N  CDz =M DI I 0 0 dropped. This was assumed in an earlier analysis (West4 erweel 1993), but it is shown in the next section that the ˆ NS …DI =ds †2 ; …13† pdf for PIV images never assumes a normal form. This p explains why earlier results required a correction in order where D ( ˆ Nd ) is the dimension of the interrogation I r to obtain a numerical correspondence with simulation region. For u ˆ v ˆ 0 and DX ˆ DY ˆ 0 (i.e., u ˆ 0), results. No general expression for K4 exists when the pdf Eq. (12) yields the square of the intensity ¯uctuation has a non-normal form. Instead, the following model is level rI: substituted


covfRD ‰s; tŠ; RD ‰s ‡ u; t ‡ vŠg 1 XX  2 fhRP ‰m; nŠihRP ‰m ‡ u; n ‡ vŠi N m n

r2I ˆ_ s^200 CDz0 =M02 Fs …0†; p  1 ds2 …14† with Fs …0† ˆ 4  fhRD ‰s; tŠjuihRD ‰m ‡ s ‡ u; n ‡ t ‡ vŠjui (Keane and Adrian 1992), which can be expressed in non‡ hRD ‰s ‡ u; t ‡ vŠjuihRD ‰m ‡ s; n ‡ tŠjuig …9† dimensional form by normalization with the mean intenwhere r2I is the intensity variance, and j4 the fourth-order sity, i.e., rI cumulant of the image intensity ¯uctuations 1=2 ˆ NS ; n …15† (j4  l4 =r4I 3, where l4 is the fourth moment of the hI i intensity ¯uctuations). The terms between braces ensure that Eq. (8) commutes. The model was obtained by gen- which is commonly referred to as the image contrast. The eralizing the original derivation for Eq. (8) (Isserlis 1918), contrast is plotted in Fig. 3 as a function of NS for different values of ds =dr . The evaluation of Eq. (8) also requires a and its reliability was tested by means of Monte-Carlo simulations. Thus, it is found that the error for the measured displacement can be determined by substitution of 1The ^-symbol indicates that the light sheet intensity is incorthe second- and fourth-order statistical properties of the porated in s^00 , which is a slight difference with respect to the intensity ¯uctuations. original de®nition of s00 j4 K4 ‰m; n; s; t; u; vŠ  hRP ‰m; nŠi 2r2I

S5

Fig. 3. Contrast and reciprocal fourth cumulant of the PIV image intensity as a function of the source density (NS) for different particle-image diameters (ds) in pixel units (dr)

S6

value for j4 in Eq. (9), which is related to the kurtosis of the intensity pdf. At high source density …NS  1† the image is a speckle pattern for which the intensity pdf is exponential (Goodman 1985), but the intensity pdf for a low source density PIV image can be quite different. It is often assumed that the pdf of low source density PIV images is bimodal, which is fed by the notion that the images consist of bright particle images on a dark background (Guezennec and Kiritsis 1990). This would be the case when: (i) the tracer particles are identical in size; (ii) the tracer particles are illuminated by a uniform light sheet, and (iii) diffraction effects are negligible with respect to the particle-image size. This is unrealistic for most practical situations, where the light source is a laser so that the intensity pro®le is (approximately) Gaussian and the particle images are diffraction limited so that the particleimage intensity has (approximately) a two-dimensional Gaussian pro®le. Three additional models for the intensity pdf of a low source density PIV image are derived in the Appendix. The results are summarized in Table 1. The pdf models and examples of synthetic PIV images that match the model conditions in Table 1 are plotted in Fig. 4. The qualitative differences between the images in Fig. 4 appear to be very small. However, note that the pdf for the case of uniform particle images with Gaussian illumination still retains a bimodal shape, but that the case of Gaussian particle images with uniform illumination and the case of Gaussian particle images with Gaussian illumination only have a single dominant peak that is associated with the dark background. Also note that the value of j4 increases substantially when the model becomes more `realistic'. Further implications are discussed in Sect. 5.

Figure 5 shows the pdf of an actual PIV image (taken by Westerweel et al. 1997). The peak to the left represents the dark image background, which is broadened by the thermal noise of the CCD sensor. This noise peak has approximately a normal pdf. The right-hand side of the histogram represents the intensity pdf of the particle images. Note that there is a good agreement between the actual histogram and the pdf model for Gaussian images of tracer particles in a light sheet with a Gaussian intensity pro®le. This indicates that synthetic PIV images in which the tracer particles are illuminated by a light sheet with a Gaussian light sheet intensity pro®le and where the particle images have a Gaussian shape indeed represent actual PIV images in a realistic manner. The value of j4 was determined by means of generating synthetic PIV images that match the pdf in Fig. 5. The results for different particle-image diameters are plotted in Fig. 3. The data can be ®tted to a straight line, given by

j4  2:7=NS

…16†

In the case of incoherent illumination Eqs. (15) and (16) are valid for all values of NS. For coherent illumination Eqs. (15) and (16) are valid only for NS ˆ 1 (viz., speckle mode) the intensity ¯uctuations have an exponential pdf (Goodman 1985), so that f ˆ 1 and j4 ˆ 6.

4 Error scaling At this point it is illustrative to ®rst consider the estimation precision for the location of a single particle image. The displacement of a particle image is given by the difference of the particle-image locations over a small time delay, so that the estimation error for the displacement

Table 1. Models for the probability density function ¦(I ) of the image intensity for low source density images with different properties, the corresponding value for the normalization constant A in terms of the scaled mean image intensity A=hI  i, and the value of j4 for a source density NS = 0.075 (f = 3.7) Light sheet

Particle images

Ref.

¦(I*)/A

Uniform

Uniform

n/a

1 A d…I  † ‡ d…I  1† A p 1  

Gaussian

Uniform

Eq. (A3)

2I

A=hI  i 2 ln I

*

Uniform

Gaussian

Eq. (A4)

1/I

Gaussian

Gaussian

Eq. (A5)

…1=I  †

p ln…1=I  †

1 p 2p

j4 8.8

1 4

16.1

1 p 1 2 2p

25.2 38.3

S7

Fig. 4a±d. Schematic representations of the image intensity pdf and examples of synthetic PIV images that match the model conditions in Table 1: a uniform light sheet/uniform particle

images; b Gaussian light sheet/uniform particle images; c uniform light sheet/Gaussian particle images; d Gaussian light sheet/ Gaussian particle images

p would be 2 times the estimation error for the location. The minimum variance of the estimated particle-image location for a given set of pixel gray values, can be determined by means of the CrameÁr±Rao lower bound, or minimum variance bound (MVB) (Van den Bos 1982; Norton 1986). The MVB is irrespective of the estimation method, and can be used to determine scaling relations for the estimation error. It is assumed that the pixel grey-level values of a single particle image are given by Eq. (1). Then the MVB matrix for any unbiased2 estimator ^t…g† of a parameter vector t for a given set of measurements g, with: gm;n ˆ cdr2 I …Xm ; Yn †, is equal to the inverse of the Fisher information matrix F, given by:

F ˆ …og=ot †T V 1 …og=ot †;

…17†

where V is the covariance matrix of the noise m de®ned in Eq. (1). It is common to assume that the particle images Fig. 5. Histogram of pixel gray values for a digital PIV image have a Gaussian shape (Adrian 1984): " # thermal noise of the CCD sensor array and the ®nite 2a …X X0 †2 ‡ …Y Y0 †2 number of gray levels. Both can be modeled as normal ; I …X; Y; t † ˆ p 2 exp 8 2 white noise, i.e., d d s 4 s

V ˆg2 1 ; …19† 2 with t ˆ …a; X0 ; Y0 ; ds †T , where a is the total particle-image where 1 is the unit matrix, and g is the noise variance. Substitution of Eq. (19) in Eq. (17) yields a diagonal maexposure, ds is the e)2 particle-image diameter, and trix for F, with a component F X 0 ;X 0 for which the reciprocal (X0,Y0) the particle-image location. Wernet and Pline (1993) determined the MVB for the yields the MVB for the estimate of the X-component of the location of a single particle image for the case of low il- particle-image location (X0,Y0). For a Gaussian particle lumination intensity where the detection is limited by the image de®ned in Eq. (18) this component is given by: " #    m0 ‡1M  actual photon count per pixel, i.e., VDiag g. It is now 2cadr 2 X2 Xm X0 2 …Xm X0 †2 common to use high-energy pulsed Nd:YAG lasers, so that F X 0 ;X0 ˆ p 2 dr exp g 4 ds mˆm 1M ds2 =16 ds2 =16 the limitation is no longer the photon count rate but the 0 2 …18†

For a biased estimator ^ t…g† with a bias ^ t0 , e.g., in the case of socalled `pixel locking', it is common to de®ne the unbiased t…g† ^ t 0 , and perform the analysis for ^ t  …g†. estimator: ^ t  …g†  ^

2



n0 ‡12N

X

nˆn0

1 2N

"

exp

# …Yn Y0 †2 dr ds2 =16

…20†

more common to vary ds while maintaining a constant For simplicity, it is assumed that M ˆ N, so that the measurements are found in a square region with sides image density. Substitution of Eq. (13) in Eq. (24) yields L( ˆ Ndr). Under the assumption that ds =dr  1, the two rDX ds 1 DI =dr sums can be re-written as two integrals (Bobroff 1986):  c ; with c ˆ ; …25†

F X0 ;X0

" S8

dr

0 ‡L=2     XZ 2cadr 2 X X0 2  p 2 g 4 ds ds =4

 exp

# …X X0 †2 dX ds2 =16 ds =4 "

# …Y Y0 †2 dX : ds2 =16 ds =4

2c^ s00 =g N 1=2 I

i.e., the error is directly proportional to ds. This result is in correspondence to the scaling suggested by Adrian (1991). For DI =dr ˆ 32; NI ˆ 16; c^ s00 ˆ 256, and g ˆ 4, a proportionality constant of c0.06 is found. This is in good agreement with c0.05±7 found by Prasad et al. (1992) from the analysis of a PIV photograph.

X0 L=2

YZ 0 ‡L=2

dr

5 Results Y0 L=2 The numerical evaluation of the mathematical expressions (Note that dX and dY each absorbed a factor dr.) Substi- derived in Sect. 2 is validated against results from Montetution of u ˆ 4…X X0 †=ds and v ˆ 4…Y Y0 †=ds yields: Carlo simulations. For particle images with a diameter that is less than two pixel units, the measurement error is  2 ‡2L=d Z s dominated by the bias error (i.e., peak locking). Therefore,
2cadr 2 2 the analytical results and the simulation results have been F X0 ;X0  u exp u du g p4 ds2 determined primarily for ds =dr  2. 2L=ds The simulations were carried out for a 32 ´ 32-pixel ‡2L=d interrogation domain. The tracer particles in the simulaZ s
exp v2 dv: …22† tion domain were `illuminated' by a light sheet with a  Gaussian intensity pro®le. It was assumed that the con2L=ds tinuous particle images have a Gaussian shape. The intensity (viz., gray value) for each pixel was computed by The two integrals are bound for 2L=d  1, with limit s p p 1 values of 2 p and p respectively. Given that the MVB integration of the continuous image ®eld. The parameters that could be varied in the simulations were: (i) the image for the particle-image displacement would be twice the MVB for the particle-image location, the following result is density; (ii) the particle-image diameter; and (iii) the particle-image displacement. In addition, simulations were obtained: carried out in which the tracer particles were illuminated p   q  p=4 ds 1 rDX ˆ 2FX0 ;X0  ds : …23† with either a uniform or a Gaussian light sheet, and in which the particle images had either a uniform or a ca=g dr Gaussian shape. So, the measurement error for the particle-image disFor the evaluation of the analytical expressions in placement would scale proportional to the square of the Sect. 2, it was assumed that the continuous particle images particle-image diameter, and inversely proportional to the have a Gaussian shape, so that the discrete displacementpixel dimension. The error is also inversely proportional correlation peak hR ‰u; vŠjui would be given by: D to ca/g (which essentially represents a signal-to-noise 2 ratio). This result is different from the result by Wernet hRD ‰u; vŠjui ˆ r/ Fr …s DX=dr †Fr …t DY=dr † ; …26† and Pline (1993), who found that the MVB would be inversely proportional to the photon count rate and directly with r/ ˆ cdr2 rI , and proportional to the particle image diameter; the difference " # Z‡1 is due to the inherent coupling of the noise and signal …f w†2 …27† …1 jwj† exp 4 2 2 dw : strength in their analysis. The result in Eq. (23) is there- Fr …f† ˆ ds =dr fore complementary to the result by Wernet and Pline. It is 1 now assumed (without proof) that Eq. (23) would also be It appeared that the results for the optimal particle-image valid for the displacement-correlation peak, so that a diameter as given by the numerical results of the equations scaling relation can be determined for rDX as a function of in Sect. 2 and the Monte-Carlo simulations were in slight 1=2 p 2 the experimental parameters. With a  rI 4 ds ˆ s^00 NS , it disagreement. This could be resolved when the effective is found that particle-image diameter for the analytical evaluation was p rDX p=4 1 ds2 increased by a factor of two. The shape of the displace : …24† ment-correlation peak that was used in the numerical 2 1=2 dr c^ s00 =g N dr S evaluation is given in Table 2. The effective source density Hence, for constant source density NS the error scales was adjusted accordingly, so that the intensity statistics, proportional to the square of ds. It was shown in Sect. 3 i.e., f and j4 are identical for corresponding results of the that a constant NS implies that the relevant image statis- analytical solution and of the numerical simulation. tics, i.e., the image contrast and the fourth-order cumuThe result for rDX/dr as a function of e/dr is plotted in lant, are invariant. However, in simulation studies, it is Fig. 6. The analytical result predicts an almost linear 

exp

…21†

Table 2. Values for Fr(s) for different values of ds/dr that were used for the computation of the theoretical results s/dr

0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2.0 2.1 2.2 2.3 2.4 2.5 2.6 2.7 2.8 2.9

ds/dr 2.0

3.0

4.0

0.637 0.627 0.599 0.554 0.497 0.432 0.363 0.295 0.231 0.173 0.125 0.085 0.056 0.035 0.020 0.011 0.006 0.003 0.001 0.001 0.000

0.783 0.775 0.750 0.712 0.660 0.600 0.533 0.463 0.394 0.327 0.265 0.210 0.162 0.122 0.090 0.064 0.045 0.031 0.020 0.013 0.008 0.005 0.003 0.002 0.001 0.000

0.862 0.855 0.837 0.806 0.766 0.717 0.661 0.601 0.538 0.474 0.412 0.352 0.297 0.246 0.201 0.162 0.128 0.100 0.076 0.058 0.043 0.031 0.022 0.016 0.011 0.008 0.005 0.003 0.002 0.001

relation between the fractional displacement and the error amplitude for jej < 12. This is in excellent agreement with the simulation results. The same shape is also found for other values of ds and NS. The linear relation between the measurement error and the displacement can be used to reduce the effort for evaluating the interrogation perfor-

mance. Suppose that the error amplitude is exactly proportional to the fractional displacement, then the mean random error for jej < 12 is equal in practice to the random error for e ˆ 0.3. Thus, for the results below, the analytical expressions were only evaluated for a sub-pixel displacement of 0.3 pixels. The results for the error as a function of ds are shown in Fig. 7. The computations were carried out for constant source density, i.e., the contrast and j4 remain constant. The solid line is the analytical result for constant source density; it corresponds quite well with the simulation results (closed symbols). For ds/dr < 2 the error rapidly increases for decreasing ds, because the shape of the displacement-correlation peak is no longer Gaussian, and thus the interpolation with Eq. (5) is no longer appropriate. (Note that an additional bias error would contribute to the error for ds/dr < 2.) The error has a minimum for ds/dr2, which is in agreement with results obtained by others, e.g., Willert (1996). For ds/dr > 2 the error increases proportional to ds2 , which is in agreement with the result for the MVB (Eq. 24) in the case of constant NS (indicated by the dotted line). For comparison, the simulation results for the case of constant image density (open symbols) are also included in Fig. 7. These results scale in accordance with Eq. (25). Figure 8 shows the error as a function of the image density. Here ds/dr is maintained constant, so that NS increases in proportion to NI. The error decays very slowly when NI becomes larger than about 5, and then becomes almost constant. In an earlier study (Westerweel 1993), the fourth-order term in Eq. (8) was neglected, and it was found that the measurement error is independent of rI. This would mean that the measurement error would be independent of the concentration of particle images in the interrogation domain; see Eq. (14). Although such behavior would be acceptable for large NI, it is clearly against any practical understanding in which the measurement error should decrease when the information content (viz., the number of particle images) increases. This behavior is easily understood if one considers the scaling of the terms

Fig. 6. The rms displacement error as a function of the particle- Fig. 7. The rms displacement error as a function of the particleimage displacement (in pixel units) image diameter (in pixel units)

S9

S10

Fig. 8. The overall rms displacement error in pixel units as a function of the image density


 2 in Eq. (7): oe=oR i  O 1=rI , and cov Ri ; Rj 
O r4I …2 ‡ j4 † . Note that r2DX  O…2 ‡ j4 †, so that it is essential to include the fourth-order correlation term in the analytical expression in order to ®nd a result for the measurement error that is a function of the image density. The present analytical result (i.e., the solid curve in Fig. 8), which includes the fourth-order correlation term, displays behavior that is valid for both large and small values of the image density and is in good correspondence with the results from the Monte-Carlo simulations. It was mentioned in Sect. 3 that a proper representation of the PIV images is required in the generation of synthetic PIV images in order to obtain the correct statistical properties for the image intensity. This was further explored by generating synthetic PIV images in accordance with the properties of those in Fig. 4. Images of each type were generated with (sub-pixel) displacements from 0.1 to 0.5 pixel units and with image densities from 5 to 15. The results for the image statistics and the random error are shown in Fig. 9. The results illustrate that the synthetic PIV images that correspond to `simpli®ed' pdf models do not yield the correct image statistics, i.e., the fourth-order cumulant is underestimated. As a consequence, the simpli®ed images also underestimate the measurement error.

6 Discussion The analytical results and simulation results presented in this paper are valid for the correlation in Eq. (2) and the three-point interpolation in Eq. (5). The question is whether these results would also be valid for other, more advanced interrogation algorithms. Such algorithms are applied mainly in non-ideal situations, e.g., PIV images with strong inhomogeneous seeding, or ¯ows with strong velocity gradients. In the case of the algorithm proposed by Huang et al. (1993) the image ®eld for a ¯ow with a strong in-plane velocity gradient is deformed to compensate for the non-uniform displacement. Hence, the correlation between the primary pre-conditioned image and the secondary image matches the conditions for which the analytical expressions were derived. Thus, the results would also have relevance for such advanced algorithms. The expression in Eq. (24) was postulated as an analogy to Eq. (23) without proof. Such a postulate is not obvious in view of the differences that exist between the statistical properties of the displacement-correlation peak and those of the pixel gray values of a single particle images. However, the aim of the analysis was to obtain expressions for the scaling behavior of the error as a function of the experimental parameters, and the justi®cation for accepting the expression as a scaling relation lies in the correspondence in Fig. 7 between the simulation results and the expressions in Eqs. (24) and (25). Both the analytical result and the simulation result suggest that the measurement error can be reduced to about 0.02 px (for NI ˆ 15, ds/dr ˆ 2, DI/dr ˆ 32) under ideal circumstances. It should be noted that it may be very dif®cult to actually achieve such ideal circumstances. For example, all particle images in the simulation have one exact identical diameter. In a practical situation, it is more likely that there is some variation in ds. Simulations show that the error can easily increase by a factor of 2±4 for a 25% variation in ds. Also, the simulations and computations were carried out for images in which no other noise , is presource other than quantization noise, i.e., g ˆ p1 12 sent. In a practical situation, it is more likely that thermal sensor noise makes a signi®cant contribution; for example, the normal noise in the histogram in Fig. 5 corresponds to g ˆ 6.6. Therefore the results in Fig. 6 and Fig. 8

Fig. 9. The image statistics and rms displacement error for different types of synthetic PIV images: , uniform light sheet/uniform particle images; (, Gaussian light sheet/uniform particle images; 4, uniform light sheet/Gaussian particle images; d, Gaussian light sheet/ Gaussian particle images; Ð, theoretical prediction

were determined for ds/dr ˆ 4, which yields errors of 0.05±0.1 px. This is also the range of measurement errors found in actual PIV measurements. The present analysis showed that fourth-order statistics have to be included in order to obtain quantitative predictions for the measurement error. For coherent imaging, both the image contrast and j4 remain ®nite for all NS, but for incoherent imaging both f and j4 vanish for NS ! 1. However, the ratio j4/f2 maintains a constant value for NS ! 1 so that K4 in Eq. (8) never becomes negligible. It is found that the source density of a PIV image effectively determines the measurement precision of the interrogation analysis, by means of the relations that exists for f and j4 as a function of NS. Figure 8 demonstrates that the error depends only weakly on the image density (for NI > 5), but it is interesting to note that the probability of detecting the displacement-correlation peak correctly is determined by the image density, as demonstrated by the extensive simulation studies of Keane and Adrian (1992). Hence, a picture emerges in which the measurement precision is characterized by the source density, and the measurement reliability by the image density. What remains is to ®nd an explanation for the increase in the effective particle-image diameter in the analytical model. This difference cannot be attributed to the difference between the particle-image diameter and the width of the displacement-correlation peak; it is easily veri®ed that the e2 -width of the Gaussian model forpthe  displacementcorrelation peak in Eq. (27) is indeed 2ds . At this moment the most plausible explanation seems to be that the model in Eq. (9) for K4 , which is the dominant term in Eq. (8) may not be fully correct. Several possibilities for the modi®cation of Eq. (9) are currently being evaluated. Nonetheless, the present model with the modi®ed particleimage diameter yields results that show good quantitative agreement with simulation results.

It is found that the measurement error can be directly linked to the second-order and fourth-order moments of the image intensity probability density function, which are both simple functions of the source density. Hence, the source density is a parameter that actually characterizes the information content of a PIV image, whereas the image density is a parameter that is associated with an interrogation area, and therefore describes the performance of the interrogation algorithm. It is demonstrated that an oversimpli®cation of the image properties in generating synthetic PIV images for Monte-Carlo simulation studies systematically underestimates the measurement error for PIV. The results in Fig. 9 illustrate that the synthetic PIV images that do not include a realistic model for the shape of the particle images and of the light-sheet intensity pro®le yield an underestimation of the fourth-order image statistics, even when the ®rst-order and second-order image statistics, re¯ected by the image contrast, are identical. Therefore, it is important that a proper algorithm is used for the generation of synthetic PIV images in order to obtain images with the correct statistical properties. The present theory can be further extended to evaluate, for example, the effect of a ®nite number of gray levels (viz., the number of bits per pixel), and the effect of sensor noise. This only requires the correct speci®cation of the secondand fourth-order statistical properties of the image signal.

Appendix In this Appendix, the model expressions for the imageintensity pdf of a low source-density PIV image are derived. At low source density, the PIV images consist of individual particle images, so that the pdf can be derived from the image intensity of a single tracer particle. Consider the ensemble of all possible locations of a single tracer particle in the object domain. It is assumed that the intensity in a ®xed point in the image domain is given by I(f), where f is a 7 scalar index that represents the tracer location. For given Conclusion I(f), the probability that the image intensity is between I The improved theory provides satisfactory predictions of and I+dI is equal to the probability that the tracer-particle the measurement error as a function of the particle-image index has a value between f and f+df, i.e., displacement (Fig. 6), the particle-image diameter 1 (Fig. 7), and the image density (Fig. 8). It is correctly f …f† ; …A1† f …I †jdI j ˆ f …f†jdfj ) f …I † ˆ jdI=dfj predicted that the measurement error for very small displacements, i.e., at sub-pixel level, is directly proportional to the displacement (Fig. 6), which is the basis for high-precision measurement of the particle-image displacement (Westerweel et al. 1997). Also, the theory predicts that the minimum error is found for a particleimage diameter of about 2 pixel units (Fig. 7), which is in agreement with earlier conjectures based on observation and numerical simulations (Adrian 1991). Another conjecture that could be con®rmed is that the measurement error is proportional to the particle-image diameter. This was already demonstrated by Wernet and Pline (1993) for the special case in which the image noise level is determined by the photon count rate; in the present paper a result was obtained for the case of additive white noise. The result shows that the error is indeed proportional to the particle-image diameter when the image density is Fig. 10. Principle of evaluating the image intensity pdf according to Eq. (A1) kept constant.

S11

S12

where ¦(I) and ¦(f) are the probability density functions for I and f respectively. This procedure is illustrated in Fig. 10. Below, the intensity pdf is derived for three different cases, namely: (1) uniform particle images in a light sheet with a Gaussian intensity pro®le; (2) Gaussian particle images in a light sheet with a uniform intensity pro®le; (3) Gaussian particle images in a light sheet with a Gaussian intensity pro®le. (1) Consider tracer particles that are illuminated by a light sheet with a Gaussian intensity pro®le. The particle images are considered to be discs with a uniform intensity, so that the intensity of each particle image is determined only by its position within the light sheet, i.e.,



I …f† ˆ_ IL ‡ IH exp

 1 2 f ; with f ˆ z=14Dz0 ; 2

…A2†

where Dz0 is the light-sheet thickness. For a homogeneous distribution of the tracer particles, ¦(f) is constant, so that the procedure in Eq. (A1) yields the following intensity pdf:

1 1 : f …I  † ˆ A  p ; for 0 < I  < 1 ; …A3† I 2 2 ln I  with I   …I IL †=…IH IL †, and where A is a normalization constant. (Note that Dz0 has been absorbed in A.) The distribution in Eq. (A3) has two sharp peaks at I = 0 and 1, so that it is bimodal. (2) For an image that is composed of identical particle images with a two-dimensional Gaussian intensity distribution, I(f) is given by Eq. (A2), but with ¦(f) ˆ 2p f, and f ˆ r=14ds , where r is the radial distance to the particleimage centroid, and ds the particle-image diameter. The pdf for uniform illumination is then given by 1 f …I  † ˆ_ A  ; for 0 < I   1 : …A4† I (Again, A absorbs all constants.) Note that ¦(I ) has only one peak at I ˆ 0, and that the pdf is `truncated' at the maximum particle-image intensity. (3) The two preceding results can be used to determine the image intensity pdf for a Gaussian particle image of a tracer particle in a light sheet with a Gaussian intensity pro®le. For a given intensity level I only (Gaussian) particle images for which the maximum intensity is larger than I contribute to ¦(I ). The shape of the image intensity pdf is identical for all Gaussian particle images, except for the truncation at I ˆ 1 which is determined by the maximum intensity of the particle images; see Eq. (A4). The pdf for the maximum intensity of the particle images is obviously given by Eq. (A3). Hence, the following expression is found for the image intensity pdf of Gaussian particle images of tracer particles in a light sheet with a Gaussian intensity pro®le: Z1 1 1 1  p dI 0 f …I †ˆA _  0 I I 2 2 ln I 0 (A5) I r 1 1 ˆ A0  ln  for 0 < I  < 1 ; I I where A¢ absorbs any integration constants.

All pdf models derived above are singular for I ˆ 0, so that they cannot be probability density functions in the strict mathematical sense. However, multiplication with I eliminates this singularity, so that the ®rst moments of the pdf models are ®nite. This means that the constants A can be expressed in terms of the ®rst moment, which is equal to the mean image intensity hI i in Eq. (11) and relates the pdf to the source density.

References

Adrian RJ (1984) Scattering particle characteristics and their effect on pulsed laser measurements of ¯uid ¯ow: speckle velocimetry vs. particle image velocimetry. Appl Opt 23: 1690± 1691 Adrian RJ (1988) Statistical properties of particle image velocimetry measurements in turbulent¯ow. In: Laser anemometry in ¯uid mechanics ± III. Instituto Superior Tecnico, Lisbon, pp 115±129 Adrian RJ (1991) Particle imaging techniques for experimental ¯uid mechanics. Annu Rev Fluid Mech 23: 261±304 Bobroff N (1986) Position measurement with a resolution and noise-limited instrument. Rev Sci Instrum 57: 1152± 1157 Goodman JW (1985) Statistical optics. Wiley, New York Guezennec YG; Kiritsis N (1990) Statistical investigation of errors in particle image velocimetry. Exp Fluids 10: 138±146 Hesselink L (1988) Digital image processing in ¯ow visualization. Annu Rev Fluid Mech 20: 421±485 Huang HT; Fiedler HE; Wang JJ (1993) Limitation and improvement of PIV. Part II: particle image distortion, a novel technique. Exp Fluids 15: 263±273 Isserlis L (1918) On a formula for the product moment coef®cient of any order of a normal frequency distribution in any number of variables. Biometrika 12: 134±139 Keane RD; Adrian RJ (1992) Theory of cross-correlation analysis of PIV images. Appl Sci Res 49: 191±215 Norton JP (1986) An introduction to identi®cation. Academic Press, London Prasad AK; Adrian RJ; Landreth CC; Offutt PW (1992) Effect of resolution on the speed and accuracy of particle image velocimetry interrogation. Exp Fluids 12/13: 105±116 Priestley MB (1992) Spectral analysis and time series, 7th ed. Academic, San Diego Van den Bos A (1982) Parameter estimation. In: Sydenham PH (ed) Handbook of measurement science, vol I. Wiley, New York, pp 331±377 Wernet M; Pline A (1993) Particle displacement tracking technique and CrameÁr±Rao lower bound error in centroid estimates from CCD imagery. Exp Fluids 15: 295±307 Westerweel J (1993) Digital particle image velocimetry ± theory and application. Delft University Press, Delft Westerweel J (1997) Fundamentals of digital particle image velocimetry Meas Sci Technol 8: 1379±1392 Westerweel J (2000) Effect of sensor geometry on the performance of PIV interrogation. In: Adrian RJ et al (eds) Laser techniques applied to ¯uid mechanics. Springer, Berlin, pp. 37± 55 Westerweel J; Dabiri D; Gharib M (1997) The effect of a discrete window offset on the accuracy of cross-correlation analysis of digital PIV recordings. Exp Fluids 23: 20±28 Willert CE; Gharib M (1991) Digital particle image velocimetry. Exp Fluids 10: 181±193 Willert CE (1996) The fully digital evaluation of photographic PIV recordings. Appl Sci Res 56: 79