Theoretical description of PIV measurement errors

Expressions for particle image velocimetry (PIV) mean error and error variance are derived for iterative deformation method algorithms. The analytical expressions explicitly account for the role of in- and out-of-plane displacements, displacement gradients, particle image diameter, fill factor of the imaging sensor, image noise, light sheet intensity distribution, seeding particle concentration, the interpolation function used to deform PIV images, and the interrogation window size and weighting window. The newly derived analytical expressions show good agreement with errors estimated using synthetic image sets.


Introduction
The understanding of measurement errors for particle image velocimetry (PIV) implementations has been dominated by empirical evidence and rules of thumb with relatively few widely applicable theoretical results. A notable exception is the work of Westerweel (2000) who derived scaling arguments for some of the error sources relevant to conventional (non-iterative) PIV algorithms. Progress on the theoretical front may have been hampered by the inherent nonlinearity of the 3-point Gaussian sub-pixel displacement estimator and anomalous scaling of the error for top-hat weighted interrogation windows (Nogueira et al. 2001;Eckstein and Vlachos 2009) when used with basic PIV algorithms. The iterative deformation method (IDM) PIV algorithm avoids some of these challenges and can be presented, at least under idealised conditions, by linear equations which lend themselves to theoretical analysis. The IDM approach is described by Astarita (2007) who derived the complete expression for its transfer function, but the basic method evolved over several years with significant contributions from Nogueira et al. (1999), Scarano and Riethmuller (2000) and many others. The essence of the approach is to iteratively estimate displacement then deform the source images to remove the apparent displacement until a converged state is reached. Other benefits of the IDM approach include robustness to displacement gradients and the ability to manipulate the transfer function to deliver a flat frequency response over a large range of scales (Cameron 2011). In this study, an expression for PIV measurement error for IDM PIV algorithms is derived and its performance evaluated against synthetic PIV images.

Theoretical analysis
In this section, an equation for the cross-spectrum of a pair of interrogation regions relevant to IDM PIV algorithms is presented. Combining this equation with an expression for the cross-correlation function after the PIV algorithm reaches a converged condition leads to the explicit definition of PIV mean error and terms representing the variance of 'aliasing', 'change of brightness', and 'image noise' errors. The contribution of spatial gradients to the measurement error is then considered.

3
where I nFn is the Fourier transform of the sampled intensity distribution of the nth particle image, R p is the Fourier transform of the random noise intensity for the pth pixel, Ĝ is the frequency response of the image interpolation function used in IDM PIV methods at the location of the nth particle ( Ĝ n ) or pth pixel ( Ĝ p ), I n = I(z n ) is the brightness of the nth particle image, I(z) is the out-of-plane light sheet intensity distribution, w is the impulse response function of the PIV algorithm sampled at the location of the nth particle w n = w(x n , y n ) or pth pixel w p = w(x p , y p ) , (x n , y n ) are the in-plane coordinates of particle centres and (x p , y p ) the coordinates of the pth pixel, and z n is the out-of-plane particle coordinate. The 'a' and 'b' subscripts denote the first and second images of a PIV image pair, respectively, ' * ' represents the complex conjugate, while N is the total number of particles and P the total number of pixels in the sampling domain. The wavenumbers k x and k y are defined as k x = 1∕ x and k y = 1∕ y (e.g. Bendat and Piersol 2000), where x and y are wavelengths measured in pixels in the respective x and y directions. For convenience, the impulse response w for the algorithm has been split into two parts w = w 0.5 w 0.5 , with w 0.5 applied to weight the contributions of particle images in the first image and w 0.5 also applied to weight the contributions of particle images in the second image (Eq. 1). Together w(x, y)I 2 (z) define a three-dimensional sampling volume, where the I is squared reflecting the product of I a and I b in Eq. 1. Equation 1 is valid for image interpolation schemes that behave as linear filters, i.e. linear interpolation, sinc interpolation, Gaussian interpolation (Kim and Sung 2006) and many others.
, which is the case for particle images with a Gaussian intensity distribution (e.g. Raffel et al. 2013), F * xna and F xnb can be written as: (1) , |F x | is the modulus of the continuous optical transfer (or frequency response) function of the imaging system, Δ x is the particle image displacement (i.e. x nb = x na + Δ x ), and the ⊕ and ⊖ symbols denote aliased components which appear due to undersampling of the particle image. Contributions from even higher wavenumbers, e.g. |F x |(k x + 2) , are in most cases negligible and are not included in the analysis. Note that the optical transfer function is the Fourier transform of the point spread (or impulse response) function of the imaging system. For a diffraction limited system, the point spread function is often approximated by a Gaussian distribution (e.g. Raffel et al. 2013) for which |F x | can be written: where d p is the particle image diameter and the effect of spatial averaging of the particle intensity distribution by the sensor pixel is incorporated through the area fill factor f f 2 which is the proportion of the sensor pixel that is sensitive to light. The distribution of |F x | could also be measured directly for a given imaging system but this requires oversampling of the images to avoid the contaminating effects of aliasing. This can be accomplished, for example, by taking images of a fixed point light source with incremental sub-pixel shifts of the image sensor. Equivalent expressions can be written for F * yna , F ynb and |F y | by substituting y for x in Eqs. 2, 3 and 4 . The function |F x | is plotted for Gaussian-shaped particle images in Fig. 1 for d p = 1.5 . It is clear that for such small particle images |F x | is non zero for |k x | > N q (where N q = 0.5 pixels −1 is the Nyquist wavenumber). In the sampled spectrum, the high wavenumber components are folded back into the range |k x | < N q and incorporated in Eqs. 2 and 3 as |F ⊕ x | and |F ⊖ x | . These aliased terms directly lead to one source of PIV measurement error. The contribution of aliasing to the (pixels -1 ) - Fig. 1 Magnitude of the optical transfer function for Gaussian-shaped particle images with d p = 1.5 and f f = 1.0 error in estimating the centroids of single particle images has previously been realised by Alexander and Ng (1991) using similar considerations. Assuming separability of the image interpolation function, i.e. G(k x , k y ) = G x (k x )G y (k y ) to greatly simplify the analysis, the frequency response of the sampled interpolation function can be written: where |G x | and − k x (Δ x + x ) are the magnitude and phase components of the continuous frequency response of the image interpolation function, Δ x + x is the measured displacement, Δ x is the actual displacement, x is the measurement error, and ⊕ and ⊖ again indicate aliased components. Equation 5 is valid for uniform displacement fields (i.e. Δ xn = Δ x and Ĝ * na =Ĝ * a ) and for IDM implementations where both 'a' and 'b' images are deformed by half the measured displacement in opposite directions. An equivalent expression can be written for Ĝ * ya by substituting y for x in Eq. 5. Using the Hamming windowed sinc interpolation function as an example, the magnitude component |G x | is computed as the cosine transform of the impulse response of the interpolation function: where L p is a parameter that can be used to control the cutoff frequency, p controls the roll-off rate of the magnitude response (it is also the half width of convolution kernel when interpolation is carried out in the space domain), and H = 0.54 − 0.46 cos 1 p (x + p) is the Hamming window function. The function |G x | is plotted in Fig. 2 for a (5) 10-pixel-wide ( p = 5 ) Hamming weighted sinc interpolation function with L p = 1.0 . This interpolation function has significant amplitude in the aliased components |G ⊕ x | and |G ⊖ x | due to undersampling of the windowed sinc function, but this can easily be corrected by selecting a L p value less than one at the expense of additional image smoothing.
The background noise in PIV images arises due to thermal effects in the imaging sensor. By treating the random additive noise at each pixel as the product of the noise intensity level r p with a delta function (x − x p , y − y p ) , its Fourier transform R p is obtained as (e.g. Smith 2013):

Converged correlation function
After several iterations of the IDM PIV algorithm, displacement estimates converge and the cross-correlation function of a pair of interrogation regions from the deformed images is centred about the zero-displacement point. The converged condition of the correlation function can be expressed by equating the correlation coefficient at the +1 pixel lag to the correlation coefficient at the -1 pixel lag, i.e. in the x direction: or and in the y direction as: where C is the cross-spectrum defined in Eq 1 and: is the inverse Fourier transform.

Errors for uniform displacement fields
Equations 9 and 10 can be solved for x and y , respectively, after substituting Eqs. 1, 2, 3, 5, and 7 combined Magnitude of the frequency response function for the Hamming weighted sinc interpolation function. p = 5 , L p = 1.0 with equivalent expressions for F * yna , F ynb , and Ĝ * ya and making use of small angle approximations cos(2 k x x ) = 1 and sin(2 k x x ) = 2 k x x . The resulting expression for is a function of random particle coordinates ( x n , y n , z n ) and the random noise level at each pixel r p . Mean error ( E[ ] ) and error variance E[ 2 ]-(E[ ]) 2 , where E is the expected value operator, can be evaluated by integrating the error probability distribution over uniformly distributed particle coordinates.
An expression for the mean error is obtained as: where T 1 , T 2 , T 3 , T 4 , and T D1 are defined in Table 1 and can be evaluated numerically after discretising at sufficiently small dk . Equation 12 indicates that the mean error is periodic in Δ x with contributions from wavelengths of both 1 and 2 pixels. The 2 pixel wavelength contributions to the mean error are generated exclusively from aliasing of the image interpolation function frequency response. If the interpolation function is designed such that there is negligible aliasing, T 1 , T 2 , and T 4 can be neglected, only T 3 will contribute, and the mean error will be sinusoidal with a period of 1 pixel. The variance of the measurement error for uniform disp l a c e m e n t f i e l d s c a n b e a p p rox i m a t e d a s can be referred to as the 'aliasing' error, the 'change of brightness' error and the 'image noise' error, respectively.
They are obtained as: where c p (particles per pixel 3 ) is the volume concentration of particles, r is the standard deviation of additive white noise in the PIV image, and the terms T 5 , T D2 , T xs , T ys , T xn , and T yn are defined in Table 1. Note that 'pixel' is used here as a unit of length rather than as a discrete image element, as is common in the PIV literature. The equivalent noise bandwidth ( E b ) incorporates the role of interrogation region size in scaling the error and was first used in a PIV context by Astarita (2007) but has not yet been widely adopted. It is defined as: where |T PIV |(k x , k y ) is the modulus of the PIV frequency response function (the Fourier transform of the PIV impulse response function). It can be calculated for IDM algorithms using the method described in Astarita (2007).
The term 2 A x can be referred to as 'aliasing' error as it arises due to undersampling of the particle images, i.e. it results from the |F ⊕ x | and |F ⊖ x | terms in Eqs. 2 and 3 . The aliasing error is a periodic function of in-plane displacement and its variance scales linearly with the inverse of particle concentration. Such scaling with particle concentration has previously been predicted by Westerweel (2000).
The term 2 B x can be referred to as 'change of brightness' (COB) error as it arises if particle images change brightness between 'a' and 'b' frames as described by Nobach and Bodenschatz (2009). The COB error involves cross-particle correlations, e.g. the product of particle n with particle n + 1 in Eq. 1, and is a function of the out-of-plane displacement ( Δ z ) but is independent of particle concentration. The terms T xs and T ys ( Table 1) are functions of the separation distance between particles ( x s , y s ) and approach zero for x s , y s >> d p (the integral limits in Eq. 14 can be reduced accordingly for numerical evaluation).
The image noise error 2 N x stems from the product of particle images and the background noise in Eq. 1. It scales proportional with the variance of noise intensity 2 r and inversely proportional to the particle concentration c p . The terms T xn and T yn , similar to T xs and T ys approach zero for x n , y n >> d p .
The terms I z1 , I z2 , I z3 are functions of the light sheet intensity distribution and for arbitrary intensity distributions are defined in Table 1. For the specific case of a Gaussian distribution described by I(z) = I max e −8z 2 ∕S 2 , where I max is peak brightness and S is the light sheet thickness at the 1∕e 2 of peak intensity level, I z1 , I z2 and I z3 simplify to: For a top-hat intensity distribution, it follows that: where S for the top-hat distribution is the full thickness of the light sheet. It is interesting to note that I z2 , which scales the error with out-of-plane displacement, is smaller for a Gaussian light sheet distribution than for a top-hat distribution if Δ z ∕S is less than about 0.3. For larger out-of-plane displacement, the top-hat distribution performs better than the Gaussian. The theoretical measurement errors given for the x displacement component by Eqs. 12, 13, 14, and 15 represent the dominant contributions to the IDM PIV measurement error for uniform displacement fields. It is straightforward to obtain equivalent expressions for the y displacement component. The theoretical distributions are approximate in nature, primarily due to the small angle approximations used in their derivation but also due to a number of minor terms that are neglected for simplicity. The theoretical error equations are evaluated against simulated PIV image sets in Sect. Results and discussion.

Errors related to spatial gradients
Thus far we have considered measurement errors associated with uniform displacement fields; however, real flows are characterised by spatial gradients which contribute to measurement uncertainty. While the IDM algorithm compensates for large-scale in-plane gradients, small-scale and out-ofplane gradients are unresolved and remain un-corrected by the iterative process. The unresolved displacement field is described by the convolution: where w 2 (x, y, z) = w(x, y) (z) , Δ x is the streamwise displacement field, Δ x (x, y, z) is the unresolved displacement field, ⊛ represents the convolution integral, and is a delta function indicating that all scales in the out-of-plane direction are unresolved. If the range of unresolved displacements within the sampling volume is small compared to the particle image diameter, the measured displacement will include a contribution from the weighted average of the unresolved displacement sampled at the random locations of the tracer particles. The variance of this weighted average is described by: where E b I z1 ∕c p can be interpreted as the inverse of the effective number of particles within the sampling volume w 3 (x, y, z) = w(x, y)I 2 (z) , and 2 G x can be referred to as the variance of the 'gradient' induced error. In the context of the measured displacement variance, the gradient error represents a partial recovery of the unresolved displacement field and is therefore not really an error in the traditional sense. However, as the energy of the partially recovered scales is allocated to the centre coordinate or the interrogation region rather than their true position, spectra of measured displacements will deviate from their true value. Equation 24 will be 1 3 evaluated against simulated PIV images in Sect. Results and discussion.

Results and discussion
Simulated data sets were obtained by processing computer generated PIV images which were created with randomly seeded particle coordinates. The synthetic particle images were sampled from a continuous Gaussian intensity distribution convolved with a 1-pixel-wide top-hat filter (i.e. a fill factor f f = 1 ). This approach to estimating PIV errors has been used extensively in other studies (e.g. Timmins et al. 2012;Raffel et al. 2013) and is the basis for much of the current understanding of PIV algorithm performance. The images were processed with a conventional implementation of the IDM PIV algorithm with 64 × 64 pixel Blackman weighted interrogation windows and using a 10 × 10 pixel Hamming weighted sinc function ( p = 5 , L p = 1 ) to interpolate the images (e.g. Astarita and Cardone 2005;Astarita 2007). Note that sinc interpolation, also known as Whittaker interpolation, is amongst the most widely used interpolation schemes for IDM PIV algorithms (e.g. Scarano and Riethmuller 2000;Foucaut et al. 2004;Astarita and Cardone 2005). The combined effect of the interrogation region size, weighting window, and velocity field smoothing between interrogation iterations results in a noise bandwidth for the algorithm of E b = 1∕25.72 2 (pixels -2 ). Here '25.72' pixels can be thought of as the equivalent top-hat weighted interrogation region size, significantly smaller than the actual size of 64 pixels due to the Blackman weighting window and the de-convolution effect of the IDM algorithm. Figure 3 shows good agreement between the mean error obtained using Eq. 12 and simulated data for a particle image concentration of c p S = 0.2 particles per pixel 2 . Note that c p S is the area concentration of particle images and can be estimated directly from PIV images by counting particles. The obtained error distributions are comparable in shape and magnitude to those of Kim and Sung (2006);Foucaut et al. (2004); Astarita and Cardone (2005) for similar PIV algorithms. The result confirms that the performance of interpolation schemes used for PIV images is fully described by |G| , at least for those schemes that behave as linear filters. We observed during testing a weak dependence of mean error on c p , particularly for small particle image diameter, with mean errors for c p S = 0.02 approximately 15% lower than for c p S = 0.2 with d p = 1.0 pixel particle images. This behaviour is not accounted for by the theoretical mean error distribution and is likely related to the small angle approximations used in their derivation.
Error variance as a function of particle image diameter, in-plane displacement, and out-of-plane displacement are presented in Figs. 4, 5, and 6, respectively. Figure 4 indicates a balance between the aliasing error which decreases with increasing particle image diameter, and the COB error which increases with diameter. The net result of the superimposed errors is an apparent minimum in the total error variance. The position of the minimum is sensitive to the out-of-plane displacement and will additionally depend on the particle concentration which scales the aliasing error. Note that for the case Δ z ∕S = 0 the COB error is zero and the total error variance is equal to the aliasing error. The aliasing error is periodic in Δ x and follows reasonably well the expected (1-cosine) distribution (Fig. 5) although the error is generally slightly under predicted by the theoretical distributions. The COB error (Fig. 6) follows fairly well the theoretical distribution which approaches the empirically determined exponential scaling of Nobach and Bodenschatz (2009) only for large out-of-plane displacements. The region Δ z ∕S < 0.2 where the error is significantly reduced can be exploited using small time separations between exposures. However, the optimum relative error will be achieved at some nonzero value of Δ z ∕S depending on the magnitude of the aliasing (and other) errors. Figure 7 shows the effects of additive image noise with standard deviation r are reasonably well described by Eq. 15. The simulations, however, indicate approximately 20% higher error variance than predicted over a range of d p . Smaller particle image diameters are found to be particularly sensitive to image noise, but for d p > 2 and typical image noise levels for PIV cameras ( r < 1 ), the image noise is unlikely to contribute significantly to total measurement error.
Equation 24 is evaluated with a 1-dimensional sinusoidal displacement field described by Δ x = A sin(2 y∕ ) with results presented in Fig. 8. Both theoretical and simulated errors were obtained by averaging in the y direction, i.e. over all phases of the sine wave. Again, a reasonable agreement between the simulated and theoretical distributions is obtained, with the error variance rapidly increasing through the partially resolved scales 100 < < 25 and then becoming constant for the unresolved scales < 25 . The error variance scales with A 2 over a range of A < 0.4 ; however, larger values of A can lead to outliers with significantly higher error levels and should be avoided by managing laser pulse separation time.
Generally, the theoretical error distributions for the aliasing, change of brightenss, image noise, and gradient errors replicate quite well the trends of the simulated data. The small differences that we do observe appear to primarily be due to the small angle approximations used to greatly simplify the derivation of the theoretical distributions.

Conclusions
Analytical descriptions of the PIV measurement error derived for iterative deformation method (IDM) PIV algorithms show good agreement with results obtained from synthetic PIV images. The theoretical error distributions can be evaluated almost instantaneously compared to the hours or days to run synthetic simulations with a sufficient number of repetitions to achieve satisfactory statistical convergence. The derived equations indicate that the mean error is periodic with in-plane displacement and arises due to aliased components, both in the image spectrum and in the frequency response of the image interpolation function. The total error variance is split into four components referred to as the aliasing error, the change of brightness error, the image noise error, and the gradient induced error. The derived error distributions fully account for the role of particle image diameter, fill factor of the image sensor, light sheet intensity distribution, seeding particle concentration, additive image noise, imposed displacement gradients and interrogation region size and weighting window.
Although the distribution of measurement errors as functions of these parameters have been known empirically for some time, the understanding gained from the newly derived theoretical error equations may lead to better experimental results through optimisation of experimental parameters. For example, the analytical error assessment can be used for a given experimental condition to estimate a range for the time separation between image exposures that minimises the relative error, or to find an optimum particle image diameter that can be configured by adjusting camera lens aperture. There is also potential to tune parameters of the PIV analysis algorithm to minimise errors by, for example, adjusting the cut-off frequency of the image interpolation function to balance aliasing errors.
The derived equations may also find a role in the assessment of measurement uncertainty in real experiments, for example, using an approach similar to the "uncertainty surface" method of Timmins et al. (2012) which uses synthetic images to estimate errors as a function of a wide range of experimental parameters. This type of assessment, however, requires as a first step knowledge of several experimental parameters such as light sheet thickness S, particle concentration c p , the particle image size d p , and the image noise r ; parameters of the analysis algorithm such as the frequency response of the image interpolation function G, and the noise bandwidth E b ; and properties of the displacement field. Although it is not yet standard to report these parameters for a given experiment, our analysis suggests that all are necessary and important for complete characterisation of an experimental setup.
Further work is required to verify the relationship between measurement errors in real PIV setups and simulated or theoretical estimates. However, this remains a challenge as the PIV measurement uncertainty for any realistic implementation cannot generally be determined. It may also be possible to extend the analysis to PIV algorithms other than IDM.