Influence of individual variations of particle image intensities on high-resolution PIV

Abstract

Individual variations of intensity of tracer particles, e.g., due to out-of-plane displacements between exposures, strongly limit the achievable accuracy of correlation-based PIV processing. The RMS error originated by this effect correlates with the spatial resolution that can be achieved with the processing algorithm making especially high-resolution algorithms like iterative image deformation affected by this error. Both aspects are shown, the gain of resolution by iterative image deformation and the loss of accuracy due to individual variations of particle intensities.

1 Introduction

Particle image velocimetry (PIV) has become the prime choice for processing image-based flow measurements in fluid dynamics experiments. The basic algorithm of digital PIV processing (Utami et al. 1991; Willert and Gharib 1991; Keane and Adrian 1992; Westerweel 1993) utilizes the cross-correlation of image sub-spaces (interrogation areas) for local displacement estimation from two consecutively acquired images of a tracer-particle-laden flow.

A variety of image processing techniques has been developed in the past to significantly improve both the accuracy of the particle displacement measurement beyond the nominal resolution of the optical sensor and the spatial resolution beyond the nominal averaging size of image sub-spaces to be correlated. These include image-deformation techniques (Jambunathan et al. 1995; Tokumaru and Dimotakis 1995; Nogueira et al. 1999, 2005a, 2005b; Lecuona et al. 2002; Scarano 2004; Astarita 2007, 2008; Schrijer and Scarano 2008), where the entire images are deformed accordingly to the assumed velocity field before the sub-division into interrogation areas to be correlated. Here the deformation’s degree of freedom is related to the grid of velocity estimates, independent of the interrogation area size. With a high overlap of neighboring interrogation areas in combination with an iterative correction of residual displacements of the estimated velocity, the spatial resolution is governed by the grid spacing without loosing the robustness of the large interrogation areas. Therefore, this method is gained to improve the achievable spatial resolution of the PIV processing. Instabilities of this technique, occurring for high overlaps of interrogation areas due to negative responses in certain frequency ranges (Nogueira et al. 1999; Scarano 2004) can be avoided either by applying appropriate spatial filters to the estimated velocity field or the application of appropriate windowing functions to the interrogation areas, which then have frequency responses with only positive values. Investigations of stability and spatial resolution of iterative image deformation applying either spatial filters or window functions can be found in Lecuona et al. (2002), Scarano (2004), Nogueira et al. (2005a, b), Astarita (2007), Schrijer and Scarano (2008).

In Nobach and Bodenschatz (2007), (2009) the variation of particle intensities has been identified as a major source of error in planar PIV, when e.g., particles with out-of-plane motion change their intensity or particles are present in only one of the consecutive PIV images or the illumination profile changes between the exposures, e.g., due to misalignments. To demonstrate the dominating influence of the intensity variation on the accuracy of correlation-based PIV algorithms, an iterative window shift method with bi-cubic spline interpolation, widely accepted as one of the best methods so far (Raffel et al. 2007; Stanislas et al. 2008), is used exemplarily in a computer simulation. All particle images get an Airy disk profile with a random maximum intensity, equally distributed between zero and 1,000 photo electrons (see comments about the noise below), corresponding to e.g., different sizes or reflectivity. The maximum intensity does not change between the exposures for only in-plane motion. With an out-of-plane motion, the particles change their position relative to the light sheet yielding different illumination of each individual particle in the two exposures. In this simulation, a top-hat profile of the light sheet is simulated, where the illumination changes only, if a particle enters or leaves the light sheet. Therefore, with an out-of-plane displacement of 1/4 of the light sheet thickness 25% of the particle images loose their associated counterpart in the other PIV image, while 75% of the particles do not change their individual intensity between the exposures. In Fig. 1, the total RMS error over the particle image diameter (defined by the first zero value of the Airy disk, obtained from a series of simulated PIV images, see Sect. 2.1) is shown for three test cases: (1) For only in-plane motion (without noise), the RMS deviation goes down to about 0.03 pixel for particle image diameters of 3 pixels and further decreases with increasing particle image diameter. For smaller particle images, the effect of under-sampling limits the achievable estimation accuracy to more than 0.1 pixel for a particle image diameter of 1 pixel. (2) For only in-plane motion, but with strong photon noise (1,000 photo electrons for the brightest particles giving about 32 electrons noise), read-out noise (RMS of 20 electrons) and quantization noise (10 electrons per count, yielding a mean gray value of 102 for the aforementioned 1,000 photo electrons), the RMS error is significantly larger than in the previous case, but still below 0.05 pixel in the range of particle image diameters of 3 pixels and beyond. (3) An out-of-plane component (here 1/4 of the light sheet thickness) has a much stronger influence than the noise. The out-of-plane motion, even without noise, limits the achievable accuracy to about 0.1 pixel at the optimum of the particle image diameter of about 3 pixels.

Fig. 1

Total RMS error of the displacement estimates obtained from a series of simulated PIV images (see Sect. 2.1) as a function of the particle image diameter (50% overlap of interrogation areas, particle number density 0.013 pixel⁻²) for in-plane motion only, in-plane motion with simulated noise and motion with an out-of-plane component (1/4 of the light sheet thickness, without noise)

In Nobach and Bodenschatz (2007), (2009) the influence of individual intensity variations has been studied for classical FFT methods including windowing, direct correlation methods and iterative methods with sub-pixel shift and deformation of the interrogation windows. In this study, the strong influence especially on high-resolution PIV methods such as iterative image deformation is investigated. Both aspects will be shown, the gain of resolution by iterative image deformation and the loss of accuracy due to individual variations of particle intensities.

2 RMS error

2.1 Simulation

A series of 100 PIV image pairs with 128 × 128 pixels size and a homogeneous displacement field over the entire image with random displacement in the range of ±1 pixel in both in-plane directions and with a varied out-of-plane component has been generated. The simulated light sheet has a top-hat profile. The particle images have random maximum intensities, equally distributed between zero and 1,000 photo electrons, and Airy disk profiles with 3 pixels diameter (defined by the first zero value of the Airy disk). The images have been analyzed with an iterative window shift and first-order deformation technique (Scarano 2002) with 32 × 32 and 16 × 16 pixels window size and an iterative image deformation with a triangular weighting applied to each PIV window of 32 × 32 pixels size. Except for the interrogation area size, the window function is identical to that in Nogueira et al. (1999), who apply the square of the triangular window to the product of the two PIV windows. To isolate the effect of decreasing the effective window size by weighting, the triangular weighting function has also been applied to the iterative window shift and deformation with a 32 × 32 pixels window. All methods use 10 iteration.

From the estimated displacement field (grid points), the total RMS error has been calculated. Note, that the validation procedure has a strong influence on the absolute values of the empiric probability of outliers as well as on the empiric RMS values. Here, a universal outlier detection (Westerweel and Scarano 2005) with a threshold of 0.5 pixel plus 2 times the found median RMS derived from the eight neighboring vectors of each velocity grid point has been used. In Fig. 2, the total RMS errors and the outlier probabilities of the four investigated estimation procedures are shown as a function of the out-of-plane displacement with a fixed 8 × 8 pixels velocity estimation grid for all methods, corresponding to 75% overlap for the 32 × 32 pixels windows and 50% overlap for the 16 × 16 pixels window, and as a function of the overlap for a fixed out-of-plane displacement of 1/4 of the light sheet thickness.

Fig. 2

a, c Total RMS error for different estimation procedures obtained from a series of simulated PIV images with a homogeneous displacement field and b, d outlier probability; a, b as a function of the out-of-plane displacement (75% overlap for 32 × 32 pixels windows and 50% overlap for 16 × 16 pixels window) and c, d as a function of the overlap (out-of-plane motion: 1/4 of the light sheet thickness, particle number density 0.013 pixel⁻²)

For all investigated methods, with increasing out-of-plane motion the RMS error increases approximately exponentially. Then, at different positions for the different methods, also the probability of outliers increases. The estimators become unreliable and outliers dominate the RMS error. For increasing overlap, the error and the probability of outliers stays almost constant over a wide range of overlaps. Then, for the window deformation methods, first the error increases rapidly and then also the probability of outliers. The maximum overlap of these methods with acceptable RMS error or outlier probability depends on the number of particles in the area that is defined by the mesh of grid points of estimated displacement vectors. A minimum of about 2 particles within one grid area is necessary for acceptable RMS errors and a minimum of about 0.5 particles per grid area for an acceptable outlier probability. For the image deformation method, the RMS error increases slightly at overlaps below 50% and stays almost constant beyond. The probability of outliers is almost constant in the entire range of overlaps. In contrast to the window deformation methods, the image deformation does not show the rapid increase in the RMS error or the outlier probability at large overlaps.

Principally, the RMS error and the probability of outliers are larger for smaller interrogation areas. So do the iterative window shift and deformation methods with 32 × 32 and 16 × 16 pixels rectangular windows. The two methods with 32 × 32 pixels triangular windows lead to results between the two other methods,where the RMS error of the iterative image deformation tends to values closer to the window shift and deformation with 16 × 16 pixels rectangular window.

2.2 Experiment

Experimental data have been taken with randomly arranged markers in a glass block moving through a light sheet generated by a video projector (Fig. 3). The glass block moves along one axis of the translation stage, while the plane of illumination is tilted with respect to the axis of motion. During the translation of the glass block with a constant velocity through the observation area of the camera, a series of 80 images of 480 × 480 pixels size has been taken. By choosing the number of frames between the two images to be correlated, different out-of-plane components can be imitated. For details of the experiment see Nobach and Bodenschatz (2009). The images are available at http://pivproc.nambis.de/data/glassblock.html. For a comparison to the previous simulation, the images taken with a 4 mm wide top-hat illumination profile with a slope of 0.75 have been re-processed in this study.

Fig. 3

a Sketch of the experimental setup, b illumination profile and c observed image (contrast enhanced)

Unfortunately, the precision of the translation stage and the motion of the glass block are not satisfactory. An a priori analysis discovered a frame to frame variation of the displacement. Additionally, a small perspective error has been found generating a velocity gradient in y direction. To compensate the displacement variations and the velocity gradient within the observation field, for each image pair, an a priori analysis with two large interrogation areas (352 × 192 pixels) with 50% overlap in y direction has been taken as a reference to derive the mean displacement and the velocity gradients in y direction. The second PIV analysis is done with standard interrogation areas (32 × 32 or 16 × 16 pixels respectively) in a 352 × 288 pixels large window, centered within the original observation area of 480 × 480 pixels. This area coincides with the area that is taken for the reference estimation. Based on the difference between the PIV analysis with standard interrogation areas and the reference estimation with large interrogation areas, the RMS error is calculated. To suppress effects from the edges, the RMS analysis uses only valid vectors from a further reduced window (160 × 96 pixels). The validation procedure is identical to that of the simulation above. For better statistics, all validated displacement vectors from all image pairs with the same number of frames between them, selected from the original series of 80 images, have been averaged.

Figure 4 shows the experimentally obtained RMS values and the outlier probabilities for the four investigated PIV procedures as a function of the out-of-plane displacement with a fixed 8 × 8 pixels velocity estimation grid for all methods, corresponding to 75% overlap for the 32 × 32 pixels windows and 50% overlap for the 16 × 16 pixels window, and as a function of the overlap for a fixed out-of-plane displacement of 1/4 of the light sheet thickness. The results agree with the previous simulation (Fig. 2), verifying both the effect of the intensity variations and the simulation procedure. Remaining deviations are possibly originated in cross-illumination of markers leading to a higher particle number density and a different distribution of particle image intensities.

Fig. 4

a, c Experimentally obtained total RMS error of the displacement estimate for different estimation procedures and b, d outlier probability; a, b as a function of the out-of-plane displacement (75% overlap for 32 × 32 pixels windows and 50% overlap for 16 × 16 pixels window) and c, d as a function of the overlap (out-of-plane motion: 1/4 of the light sheet thickness)

3 Resolution

To prove the gain of resolution by image deformation, a series of 100 pairs of PIV images with 512 × 512 pixels each has been generated with a random in-plane displacement on a pixel-resolution (Gaussian distribution for each component and for each pixel with an RMS value of 0.5 pixel) and no out-of-plane motion. The particle images have random maximum intensities, equally distributed between zero and 1,000 photo electrons, and Airy disk intensity profiles with 3 pixels diameter (defined by the first zero value of the Airy disk). To obtain a high spatial resolution and estimation stability, the particle number density is chosen higher (0.05 pixel⁻²) than in the previous simulation. The images have been analyzed with the four procedures as above.

From the individual displacement estimates, which are interpolated with bi-cubic splines and re-sampled at all pixel positions, and the simulated displacement, which originally is given for all pixel positions, a two-dimensional coherent frequency response

$$ C_{ij}={\frac{\left\langle U^\ast_{{\rm est},ij} U_{{\rm sim},ij} + V^\ast_{{\rm est},ij} V_{{\rm sim},ij}\right\rangle} {\left\langle U^\ast_{{\rm sim},ij} U_{{\rm sim},ij} + V^\ast_{{\rm sim},ij} V_{{\rm sim},ij}\right\rangle}} $$

(1)

is calculated, where U _sim,ij and V _sim,ij are the two-dimensional Fourier transforms of the simulated u and the v displacement fields, U _est,ij and V _est,ij are the estimated counterparts, the asterisk denotes the conjugate complex and 〈〉 denotes the ensemble average. The products and the coherent frequency function are calculated element-wise for the two-dimensional functions. From the two-dimensional coherent frequency response function, a common (one-dimensional) one is derived by iteratively optimizing a one-dimensional function c _i so that the component-wise products c _i c _j fit best the two-dimensional function C _ij with minimum L₂ norm.

Figure 5a shows the frequency response function for the four investigated estimation procedures all using a velocity estimation grid of 8 × 8 pixels. With a rectangular weighting window, the frequency response clearly drops below zero at 1/(32 pixel) or 1/(16 pixel) corresponding to the interrogation area size of 32 × 32 or 16 × 16 pixels, respectively. The triangular weighting window applied to a 32 × 32 pixels interrogation window leads to a frequency response function with only positive values, while the resolution increases beyond the nominal resolution of the interrogation window size, reaching almost an effective window size of half the nominal window size. The image deformation technique can further improve the spatial resolution, which then is limited by the velocity grid of 8 × 8 pixels. Figure 5b shows the obtained bandwidth (−3 dB limit) as a function of the overlap of interrogation areas. Clearly, the image deformation technique gains most by increasing the density of the velocity estimation grid. Note, that the overlap of interrogation areas for a given grid of velocity estimates changes with the size of the interrogation areas yielding a shifted overlap for the method with the 16 × 16 pixels interrogation area compared to the other methods.

Fig. 5

a Coherent frequency response for the different estimation procedures for a velocity estimation grid of 8 × 8 pixels corresponding to 75% overlap for 32 × 32 pixels windows and 50% overlap for 16 × 16 pixels window; b bandwidth as a function of the window overlap obtained from a series of simulated PIV images with a random in-plane displacement field (no out-of-plane motion; particle number density 0.05 pixel⁻²)

4 RMS error versus resolution

Combining the results from Figs. 2c and 5b, the correspondence between obtainable resolution and the RMS error can be derived. Therefore, the first simulation has been repeated with a higher (0.05 pixel⁻², matched to Fig. 5b) particle number density without and with out-of-plane motion. In Fig. 6a, b, the total RMS error and the outlier probability of the four investigated estimation procedures are shown as a function of the overlap similar to Fig. 2c and d with only in-plane motion and a particle number density of 0.05 pixel⁻². Both, the error and the outlier probability are significantly smaller than in the previous simulation. The influence of the overlap on the results of the window shift and deformation methods is similar to the previous simulation: Below a particle number density of about 2 particles per grid area the error and below about 0.5 particles per grid area also the probability of outliers increase rapidly. The error of the image deformation method increases only slightly and outliers could not be detected at all in this simulation.

Fig. 6

a Total RMS error and b outlier probability for the different estimation procedures as a function of the window overlap obtained from a series of simulated PIV images with a homogeneous displacement field; c total RMS error against the bandwidth (only in-plane motion; particle number density 0.05 pixel⁻²)

Figure 6c shows the obtained total RMS errors against the bandwidth with overlaps of interrogation areas varied between 0 and 87.5%. The various methods cover different ranges of obtainable bandwidths and RMS errors yielding a lower bound of about 0.02 pixel, slightly increasing with the obtainable bandwidth.

For the window shift and deformation techniques with rectangular window functions, the achievable bandwidth basically depends on the size of the interrogation areas. The bandwidth increases slightly with the overlap up to about 50%. For higher overlaps, the bandwidth stays constant. As seen in the previous cases, for further increased overlap the RMS error rapidly increases at about 2 particles per grid area. For these methods, the optimum choice of PIV processing is given for 50% overlap with a minimum of 8 particles per interrogation area.

For window functions applied, the bandwidth further increases, e.g., for the triangular window function up to about 75% overlap, reaching almost the bandwidth of the window shift and deformation with a rectangular window function of half the size corresponding to an “effective” window size, which is half as large as the nominal size. Again, the RMS error increases rapidly at about 2 particles per grid area. For this method, the optimum choice of PIV processing is given for 75% overlap with a minimum of 32 particles per interrogation area.

With the image deformation technique and strong overlap of interrogation areas, the bandwidth can be increased further. The RMS error increases much less then with the other methods and follows an extrapolation of the results of the other methods. The probability of outliers stays small even for strong overlaps of the interrogation areas.

A comparison of the various methods shows that the RMS error corresponds to the achievable bandwidth, where the error slightly increases with the bandwidth independent of the method. The bandwidth itself is determined by the PIV method and the size of the interrogation areas. Clearly, the image deformation is the best choice. It realizes highest bandwidths at low RMS errors and outlier probabilities. The bandwidth can be arbitrarily chosen by the overlap of the interrogation areas, even beyond the limit of 2 particles per grid area of the window shift and deformation methods. Of course, a reasonable limit is given by the computational load that comes with high overlaps.

The picture changes completely in the case of an out-of-plane component. Figures 7 and 8 show the results similar to Figs. 5 and 6, but for an out-of-plane component of 1/4 of the light sheet thickness. The frequency response and also the bandwidth in Fig. 7 do not change with the out-of-plane component compared to Fig. 5. However, the RMS error and the outlier probability are much larger than in the case of only in-plane motion (Fig. 8a, b compared to Fig. 6a, b). Independent of this, again a rapid increase of the RMS error can be seen below 2 particles per grid area and a rapid increase of the outlier probability below 0.5 particles per grid area for the window shift and deformation techniques. For the image deformation, the RMS error and the probability of outliers increase slightly, mainly around 50% overlap and keep almost constant otherwise.

Fig. 7

a Coherent frequency response for the different estimation procedures for a velocity estimation grid of 8 × 8 pixels corresponding to 75% overlap for 32 × 32 pixels windows and 50% overlap for 16 × 16 pixels window; b bandwidth as a function of the window overlap obtained from a series of simulated PIV images with a random in-plane displacement field (out-of-plane displacement of 1/4 of the light sheet thickness; particle number density 0.05 pixel⁻²)

Fig. 8

a Total RMS error and b outlier probability for the different estimation procedures as a function of the window overlap obtained from a series of simulated PIV images with a homogeneous displacement field; c total RMS error against the bandwidth (out-of-plane displacement of 1/4 of the light sheet thickness; particle number density 0.05 pixel⁻²)

Figure 8c shows that also here the RMS error increases with the achieved bandwidth. However, with the out-of-plane component the errors are much larger than with only in-plane motion and, additionally, the results for the various methods do not overlay any more. Without changes compared to Fig. 6c, the achievable bandwidth and the RMS error is larger for a smaller interrogation area (16 × 16 pixel compared with 32 × 32 pixel) using the window shift and deformation method with a rectangular window function with an optimum of about 50% overlap. Again, the window shift and deformation with a 32 × 32 pixel interrogation area and a triangular window function achieves almost the same bandwidth as with the 16 × 16 pixel large rectangular window function. The RMS error here is larger than in the previous simulation with only in-plane motion, but significantly smaller than the error with the 16 × 16 pixel large rectangular window function. The optimum still is at about 75% overlap. Again the image deformation reaches the highest bandwidth at strong overlaps, but here for the prize of a larger RMS error in the entire range of overlaps and bandwidths compared to the other methods.

Therefore, in the case of an out-of-plane component, instead of the image deformation the window shift and deformation with a window function (here the triangular window with 75% overlap) is the best choice. It reaches almost the resolution of the window shift and deformation with a rectangular window with half the size, while the RMS error is smaller than that of the other methods. Generally, the application of a window function can be recommended, since it comes with a smaller error than decreasing the size of a rectangular window.

5 Conclusions

The effect of particle image intensities varying individually between the two consecutive images on the obtainable accuracy of different types of PIV image processing with different spatial resolution has been investigated. The RMS error strongly correlates with the achievable spatial resolution for all investigated methods of PIV image processing. Especially image deformation PIV methods, which are able to achieve highest spatial resolution, are strongly affected by the error due to changing particle intensities. Since no mean against this error has been found so far, the achievable resolution and accuracy must be balanced by either minimizing the out-of-plane displacement between exposures or using large enough averaging areas. The later, however, contradicts the purpose of high-resolution PIV methods as e.g., the iterative image deformation. The best compromise seem to be iterative window shift and deformation techniques with weighting functions. This type of PIV image processing has a low RMS error, while the achievable spatial resolution can be set in a wide range by choosing an appropriate overlap and size of the interrogation windows. If one finds a method to reduce the PIV error due to the individual intensity variations of particle images, these results must be reviewed.