Elimination of unsteady background reflections in PIV images by anisotropic diffusion

A novel approach is introduced that allows the elimination of undesired laser light reflections from particle image velocimetry (PIV) images. The approach relies upon anisotropic diffusion of the light intensity, which is used to generate a background image to be subtracted from the original image. The intensity is diffused only along the edges and not across the edges, thus allowing one to preserve, in the background image, the shape of boundaries as laser light reflections on solid surfaces. Due to its ability to produce a background image from a single snapshot, as opposed to most methods that make use of intensity information in time, the technique is particularly suitable for elimination of reflections in PIV images of unsteady models, such as transiting objects, propellers, flapping and pitching wings. The technique is assessed on an experimental test case which considers the flow in front of a propeller, where the laser light reflections on the model’s surface preclude accurate determination of the flow velocity. Comparison of the anisotropic diffusion approach with conventional techniques for suppression of light reflections shows the advantages of the former method, especially when reflections need to be removed from individual images.


Introduction
Particle image velocimetry (PIV) images are often affected by unwanted light reflections occurring when the laser light impinges on a solid surface. The intensity of those reflections can be one order of magnitude larger than that of the particle images, causing a high auto-correlation peak in the correlation function. Such a peak can be much higher than the particle image's displacement cross-correlation peak, thus precluding accurate determination of the flow velocity in the proximity of a solid surface.
Several approaches have been devised to avoid laser light reflections when conducting the measurements. Whenever possible, it is good practice to cover the model with matt black paint, so that most of the laser light impinging on a surface is absorbed instead of being reflected (Gui et al 2001). Fluorescent paint (e.g. rhodamine) can be applied to the model surface to change the wavelength of the reflected light from green to red (Depardon et al 2005). A bandpass filter mounted on the camera lens allows one to reject the red light from the surface, thus retaining only the green light scattered from the particles. However, in many cases the model cannot be painted, either because of the presence of wall tapping for pressure measurements, or so as not to alter the surface roughness and consequently the boundary layer properties. For flat surfaces, Kähler et al (2006) report that tangential model illumination allows a dramatic suppression of undesired wall reflections. Nevertheless, in the presence of more complex model geometries, where the model surface presents a curvature, tangential illumination cannot be achieved on the entire surface, but only Measurement Science and Technology

Elimination of unsteady background reflections in PIV images by anisotropic diffusion
at one specific location of the surface. The influence of the imaging angle was investigated by Lin and Perlin (1998). The authors report that, for measurements in water flows with a free surface, tilting the camera to the air-water Brewster angle has the effect of removing most of the reflections from the water-free surface. More recently, Kähler (2009) investigated the effect of the model material and surface treatment on the intensity of light reflections; the author found that aluminum models with highly polished surfaces have the minimum diffusive reflectivity among the tested materials (steel, carbon fiber reinforced plastic, glass, PMMA).
Despite the efforts above, in many cases laser light reflections are still present in the PIV recordings and need to be treated in the pre-processing phase (image restoration). The objective of image restoration is to remove the unwanted background from the images while keeping the particle image's signal. In some cases, the background removal is simply achieved by recording a background image without tracer particles, and then removing it from the PIV recordings. Even if the background image has not been acquired, in stationary problems (no moving interfaces or pulse-to-pulse light intensity variation) it is possible to generate such background images by image statistics, computing the minimum or the average of the light intensity at each pixel location (Adrian and Westerweel 2011). However, it is well known that the subtraction of time-average intensity may lead to removal of particle images from the recordings, especially in the low-velocity regions where the particles displacement is the minimum. For slowly moving light reflections or light intensity variations at a frequency much lower than the acquisition frequency, subtraction of a time-varying background, generated by sliding-average or sliding-minimum light intensity over a short kernel, can be employed. Theunissen et al (2008) assessed different pre-processing techniques for background reflection removal using Monte Carlo simulations. The authors reported that a combination of local minimum intensity subtraction and equalization of mean intensity is effective in removing reflections that are constant in time. Alternatively, Sciacchitano and Scarano (2014) proposed the use of a temporal high-pass filter to remove unsteady background reflections while retaining the particle image intensity. The approach is based on the decomposition of the signal in the frequency domain and the removal of the low-frequency content, representative of the unwanted light reflections. The underlying assumption of the approach is that the contribution of the reflection (low-frequency) is well separated in the frequency domain from the contribution of the particle images (high-frequency), meaning that the reflection must reside for a few time instants in a pixel location. Recently, Mendez et al (2017) proposed a proper orthogonal decomposition (POD)-based background removal, which can in principle also eliminate reflections on moving surfaces, provided that a sufficiently large ensemble of images is available for convergence of the POD modes. However, even this approach requires that the reflection resides for several recordings in the same pixel locations, which is often not the case in the presence of towed or moving models, flapping or rotating wings.
In the latter cases, background removal cannot rely upon image statistics, but must be conducted on individual raw images. Most approaches rely on the consideration that the particle images have a shorter length scale (typically 1-5 pixels) than the background reflection, which may cover tens of pixels. Hence, the contribution of the particle images can be isolated from that of the background by applying a spatial high-pass filter, where the filter kernel should have a linear size at least as large as the particle image diameter. Several filters have been proposed in the past, including the top-hat slidingaverage filter, the Gaussian filter, the median filter (Adrian and Westerweel 2011) and the min/max filter (Westerweel 1993). Nevertheless, the use of an isotropic filter, which has the same effects in all directions, typically yields low performance in the proximity of sharp reflections, causing reduction of the signal level. Honkanen and Nobach (2005) proposed a simple background extraction approach for double-frame PIV images, where the second frame of the image pair is subtracted from the first frame. The idea is that everything that stays stationary in the image pair, namely the background, is removed from the first image. However, this approach may lead to particle cancellation in cases of high source density or in regions where the flow velocity is the lowest. Another approach proposed by Deen et al (2010) to eliminate moving reflections relies upon the combination of several image processing techniques, such as intensity normalization, background subtraction and masking. The approach was successfully employed by the authors to remove the undesired correlation peak due to non-stationary bubbles in a two-phase flow. However, the requirement to use several techniques makes this approach computationally expensive for most practical applications. Mejia-Alvarez and Christensen (2013) modified the algorithm proposed by Honkanen and Nobach (2005) by computing the normalized local intensity with respect to the difference between sliding median and minimum intensities. Although this algorithm is able to suppress the residual background reflections which are not eliminated by Honkanen and Nobach's approach, its performances have been demonstrated only for diffused reflections due to an irregular rough wall, and not in the presence of sharp reflections occurring: e.g. when the laser sheet impinges on a solid surface.
The discussion above shows that an effective methodology for the removal of the unwanted laser light reflections from individual PIV images, thus without making use of image statistics, is currently lacking. Such an approach would find its application in PIV measurements where the laser light reflection is unsteady, e.g. due to the presence of transiting objects, propellers, flapping or pitching wings.
In the image processing community, approaches for edge detection based on anisotropic diffusion have been widely used over the last three decades since the seminal paper of Perona and Malik (1990). The idea is to compute a slidingaverage of the intensity of an image on an anisotropic kernel, which accounts for the intensity gradient. The approach has been successfully employed to enhance edges with respect to background noise. Further improvements to the technique have been proposed by Chao and Tsai (2006) for restoration of astronomical images. In their case, the image of the Henize 70 nebula was obscured by sparking stars. The approach was employed to segregate the stars, which had a low length scale of the order of a few pixels, from the nebula characterized by a large length scale and low intensity gradient.
In PIV images, the light reflections are sharp and usually have higher intensity levels than the particle images, as opposed to nebula images. Hence, in this paper, we further develop the approach of Perona and Malik (1990) and Chao and Tsai (2006) for isolating the contribution of the particle images from that of unwanted light reflections.

Proposed methodology
To explain the technique, consider a raw image I R (x, y ) where both unwanted laser light reflections and tracer particle images are present. A background image I BG (x, y ), ideally containing only the unwanted laser light reflections and no tracer particle image, is often computed via the sliding-average of the intensity of I R (Adrian and Westerweel 2011). The latter operation is typically conducted by convolution of I R with a kernel G, which is usually top-hat or Gaussian. As pointed out by Koenderink (1984) and Hummel (1987), I BG can also be interpreted as the solution I(x, y , t) of the diffusion equation: with the initial condition I(x, y , t = 0) = I R (x, y ). Note that in equation (1) the intensity is diffused isotropically, with no preferential direction (i.e. diffusion occurs at the same speed in all directions). Also, equation (1) treats in the same way particle images, which cover only a few pixels, and reflections, which typically affect several pixels. When equation (1) is applied to typical PIV images with reflections (figure 1(a)), a clear smoothing of the unwanted laser light reflection is noted (figure 1(b)). As a consequence, the solution of equation (1) is not a good estimate of the true background image. Following Perona and Malik (1990), equation (1) can be rewritten into the anisotropic diffusion equation, so that diffusion occurs along the edges and not across the edges: (2) Perona and Malik (1990) proposed to choose the diffusion coefficient c as a function of the magnitude of the intensity gradient: with g being a suitable monotonic function. The authors used the following expression for g: where K is a positive constant termed as a threshold parameter. When g is chosen as a monotonically decreasing function, little diffusion occurs in the direction of the high intensity gradient, e.g. at the interface between the laser reflection and the fluid region. Conversely, the diffusion mainly occurs in the direction of low intensity gradients, i.e. along the light reflection. The approach proposed by Perona and Malik (1990) is very effective in avoiding smoothing of the edges, but it only considers the magnitude of the intensity gradient, and not the local intensity. As a consequence, it does not cause diffusion of small intense particle images, which remain in the estimated background image (figure 1(c)).
Modifications to equation (4) have been proposed by Tsai (2006, 2010) to account not only for the intensity gradient, but also for the local intensity variance. In the present work, the diffusion coefficient is also computed as a function of the local normalized intensity I n to enable the distinction between reflections, which cover several pixels in an image, and pointwise bright spots such as the particle images. I n is evaluated as the local intensity normalized with respect to the local mean of the intensities (computed with respect to 12 neighbors in a diamond shaped kernel; the neighbors are defined by D 8 distance = 1 and D 4 distance = 2, as described in Gonzalez and Woods (2002)): The particle images are typically characterized by large values of the normalized local intensity I n compared to the magnitude of the intensity gradient |∇I|, whereas the reflections feature small values of I n with respect to the corresponding |∇I|. It is to be noted that the normalized local intensity I n and the magnitude of the intensity gradient |∇I| can be compared directly, since |∇I| is defined in the discretized form as the difference between the intensities of the neighboring pixels. Thus, the diffusion coefficient is large for the particle images and small for the reflections, as shown in figure 2. This choice enables high diffusion for the particle images, and as a result, the particle images are no longer present in the background image (figure 1(d)).

Numerical implementation
Following Perona and Malik (1990), equation (2) is discretized as follows: where (i, j ) are the pixel locations along the y and x directions, respectively, 0 ⩽ λ ⩽ 0.25 for numerical stability, the subscripts N, S, E and W represent North, South, East and West, and ∇ indicates the nearest-neighbor differences: In this work, λ = 0.2 is used in all analyses. The diffusion coefficients are updated at each time instant as a function of the local intensity gradient and normalized intensity level:

Selection of the threshold parameter and number of iterations
To solve the anisotropic diffusion equation (2), the value of two relevant parameters must be selected: the threshold parameter K and the number of iterations t f . A parametric study is conducted to determine which combination of K and t f is the most effective for background removal in PIV images. First, the effect of the threshold parameter is studied by considering different values of K for the same number of iterations. Background and pre-processed images of a typical PIV raw image are illustrated in figure 3 for t f = 300 iterations and K equal to 5, 10 and 50, respectively. It is observed that for small values of the threshold parameter (K = 5), the particle images are not diffused sufficiently and therefore are still present in the background image (first column in figure 3). Conversely, a large value of K (K = 50 in the example), causes diffusion of the sharp reflection along with the particle images. Hence, the reflection is not eliminated sufficiently in the pre-processed image obtained by subtracting the background image from the original raw image (last column in figure 3). The results can be explained based on the definition of the diffusion coefficient (equation (6)), where the large value of K makes the diffusion coefficient approach unity. In the latter case, the diffusion process becomes isotropic, as expressed in equation (1). It is observed that an intermediate value of K (K = 10) yields better results than those for K = 5 and K = 50, by diffusing the particles sufficiently and by retaining the sharp reflection in the background image (middle column in figure 3). The results in figure 3 thus suggest using an intermediate value (K = 10) for the threshold parameter in the proposed anisotropic diffusion approach.
To investigate the effect of the number of iterations, different values of t f are considered, keeping the threshold parameter constant (K = 10). The background and pre-processed images are shown in figure 4 for the three cases of t f equal to 10, 300 and 1000. It is observed that for a small number of iterations (t f = 10), the particle images are not diffused completely in the background image (first column in figure 4), yielding a pre-processed image where the signal is strongly attenuated. Conversely, for a large number of iterations (t f = 1000), the reflection is diffused in the background image (last column in figure 4), and therefore remains partly present in the pre-processed image. When t f = 300 iterations is employed, the reflection remains sharp in the background image, whereas the particle images are diffused, yielding better background removal without loss of signal from the tracer particles (middle column in figure 4). Based on the above, the combination (K, t f ) = (10, 300) is more suitable than the other combinations in removing the background of the sharp reflections.
The effect of the number of iterations on the intensity levels of particle images and reflections is further illustrated in figure 5 for different values of K. Two windows of 10 × 10 pixels are considered representative of the reflections region and of the particle images region, respectively, as illustrated in figure 5 (left). After each iteration of the anisotropic diffusion  algorithm, the intensity in each of the two windows is computed and plotted in figure 5 (right). The intensity of the reflection is calculated as the mean intensity in window 1, whereas the particle's intensity is calculated as the maximum intensity in window 2, since the maximum intensity level represents the particle image peak intensity. Figure 5 (right) shows that the rate of diffusion is high for K = 50, causing the reflection to diffuse along with the particle images, which is not desirable. Instead, for K = 5 the diffusion is very slow and it takes more iterations to attenuate the particle image's intensity compared to the other two values of K. Thus, K = 10 is found to be a good choice for the threshold parameter. The plots for K = 10 show that the particles are removed sufficiently after about 300 iterations; further increasing the number of iterations does not produce any improvement in the background image. In contrast, the reflection intensity reduces slightly with an increasing number of iterations. Hence, a larger number of iterations has the effect of causing diffusion of the laser light reflections, returning an output image that is not representative of the actual background. For this reason, the number of 300 iterations in combination with a threshold parameter K = 10 is considered a good choice to generate the background image. It is to be noted that the values are not the optimum values, although they are shown to be effective. The reader is advised to plot curves, as in those of figure 5-right, for a pair of PIV images to find out the suitable values of K and t f . Then, the same values could be used for pre-processing all the images of a set. A rule of thumb is to select the value of K such that it provides a high slope for the curve of the particles' intensities and a low slope for that of the reflections' intensities; the number of iterations is selected as the minimum value of t f for which the particles' intensities are below a certain threshold (e.g. five counts).
It should be noted that the algorithm is not very sensitive to the choice of the processing parameters K and t f , in the sense that a variation in these parameters by 10%-20% would in practice yield the same background image. The computational time is proportional to t f and is comparable to that of other standard filters.

Experimental assessment
The performance of the proposed anisotropic diffusion approach is assessed via PIV images acquired for the investigation of the propeller blade vortex interaction (Yang et al 2016, Yang 2017. This particular experimental test case is   The flow is seeded with micron-sized water-glycol particles produced by a SAFEX Twin Fog Double Power smoke generator. A detailed description of the experimental setup can be found in Yang (2017). Figure 7 shows a raw image pair and the corresponding displacement field (from a single camera) obtained in this experiment. As can be seen in figure 7-left, the raw images are affected by strong laser light reflections, especially at the leading edge of the propeller blade. Since the propeller is spinning at 2500 rpm, the propeller tip moves about 2 mm between the two image frames. Note that the measurement plane is located about 12 mm upstream of the blade leading edge; as a consequence, the fluid displacement in the measurement plane differs from the blade displacement. The cross-correlation analysis on the raw images returns a displacement field (figure 7-right) that is highly affected by the laser light reflections. The flow displacement in front of the propeller blade is highly over-estimated due to the presence of the blade reflection, which moves between the two frames.
It should be noted that, since the PIV acquisition was not synchronized with the rotation of the propeller blade, the position of the latter varies among different recordings (see figure 8). As a consequence, standard background removal approaches based on the statistical analysis of the sequence of images (e.g. subtraction of the time-average or time-minimum intensity) fail in removing the background reflection. Furthermore, even more advanced approaches based on image statistics, such as the POD filter (Mendez et al 2017), are not effective in this specific case, due to the limited number of recordings (250) per set of images, which causes the POD modes to not reach statistical convergence.
When image pre-processing is performed, the relative intensity of the particle images with respect to the laser reflection can be highly enhanced. Figure 9 shows a comparison among the image pre-processing by sliding-average subtraction, median filter subtraction, median-based-normalization subtraction and the proposed anisotropic diffusion approach. In the first method, the background image (figure 9-first row-first column) is built as a sliding-average (i.e. isotropic diffusion) of the image intensity in a kernel of 3 × 3 pixels in 30 iterations. A Gaussian weighting is applied to the intensity within the kernel. In the second method, a  median filter of a kernel of 5 × 5 pixels is applied to the raw image to generate the background image (figure 9-second row-first column), which is then subtracted from the raw image to get the pre-processed image (figure 9-second rowsecond column). The next method is based on Mejia-Alvarez and Christensen's (2013) approach, where normalization is performed with respect to the local median and minimum intensities in a kernel of 9 × 9 pixels. The two frames are then subtracted from each other to eliminate the background and further normalization is applied with respect to local maximum and minimum intensities. Finally, the intensity values are stretched according to the global maximum and minimum intensities in the original raw images. The background and pre-processed images obtained with this medianbased-normalization algorithm are shown in figure 9 in the third row, first and second columns, respectively. As can be seen in figure 9 (rows-1, 2, 3), these three methods (i.e. pre-processing by sliding-average subtraction, median filter subtraction and median-based-normalization subtraction) yield background images where the particle image intensity is highly reduced. However, the light reflection on the propeller blade is also diffused with respect to the raw images. As a consequence, when the pre-processed images are evaluated from the difference between the raw images and background, they still feature laser light reflections, which yield erroneous vectors in the displacement field (figure 9-rows-1, 2, 3-last column). Instead, when the background image is built with the proposed anisotropic diffusion approach, the particle image intensity is strongly diffused, whereas no significant diffusion occurs in the light reflections on the propeller blade (figure 9-last row-first column). As a result, the background image by anisotropic diffusion is much more representative of the true background. In the pre-processed image, the intensity of the laser light reflections on the propeller blade becomes lower than that of the particle images (figure 9-last row-second column). Hence, the computed displacement field does not feature any erroneous vector associated with unwanted laser light reflections on the solid surface (figure 9-last row-last column).
For a quantitative assessment of the performance of the anisotropic diffusion filter, the cross-correlation analysis is conducted in two interrogation windows of 65 × 65 pixels, shown in figure 7 (left). The interrogation window 1 is located in front of the blade and features strong laser light reflections, whereas interrogation window 2 is in a region free of any unwanted reflections. The results of the cross-correlation analysis are illustrated in figures 10 and 11, and the corresponding pixel displacements and cross-correlation signal-tonoise ratios (SNR) are reported in table 1.
In interrogation window 1, the strong reflection on the propeller blade yields a high peak in the cross-correlation functions obtained from raw images and pre-processing by sliding-average subtraction or median filter subtraction (figures 10(a)-(c)). The position of the peak corresponds to the displacement of the propeller blade between frame 1 and frame 2. As mentioned before, such displacement is not the same as the fluid displacement, because the plane of the propeller does not coincide with the measurement plane. Such a peak is much larger than the true particle displacement peak (which can be seen around Δx = 6 pix, Δy = 4 pix) yielding a correlation SNR smaller than one. As a consequence, an erroneous displacement vector is estimated. Conversely, in the case of pre-processing by median-based-normalization or anisotropic diffusion, the correlation peak due to the blade movement is considerably attenuated or not even visible (figures 10(d) and (e)). Hence, the particle image's displacement peak is correctly identified, yielding a valid vector estimation. However, a comparison of the correlation SNRs from table 1 shows that the particle's signal is attenuated more with the median-based-normalization algorithm (SNR = 1.5) than with the anisotropic diffusion approach (SNR = 4.0).
When the cross-correlation analysis is conducted in a region free of any reflections (i.e. interrogation window 2), all pre-processing methods correctly identify the displacement peak (figure 11), leading to a displacement estimate accurate within 0.1 pixels. However, it is noted that the image pre-processing with the median-based-normalization algorithm strongly attenuates the particle's signal, returning a relatively low SNR which results in slightly inaccurate displacement measurements. On the contrary, the SNR obtained with the anisotropic diffusion filter is approximately the same as that achieved with the raw images, which indicates that the   A comparison among image pre-processing by sliding-average subtraction (kernel of 3 × 3 pixels and 30 iterations), median filter subtraction (kernel of 5 × 5 pixels), median-based-normalization subtraction (kernel of 9 × 9 pixels) and the proposed anisotropic diffusion approach (diamond shaped kernel, K = 10 and 300 iterations). Background images (first column; red: first recording and green: second recording), pre-processed images (second column) and corresponding displacement fields (third column).  approach has no detrimental effect in regions where no reflections are present.

Conclusions
A novel approach is proposed to suppress undesired light reflections from PIV images. The approach relies upon generating a background image by anisotropic diffusion of the intensity distribution of the raw image. The principle is that, by means of anisotropic diffusion, the image intensity is diffused only along the edges and not across the edges, maintaining sharp reflections in the background image. The latter is then subtracted from the original image, yielding a pre-processed image where no reflection is present and only the contribution of the particle images is retained. Contrary to most approaches for background removal that require the analysis of an image sequence (e.g. subtraction of time-average or time-minimum image intensity, POD filter, high-pass filter in the frequency domain), the proposed approach is applicable to individual images, and is therefore suitable for all cases where the reflection is unsteady, or when a short image sequence has been acquired, yielding lack of convergence in the statistical analysis.
A parametric study has been conducted to evaluate the effect of two key parameters of the approach, namely the threshold parameter K and the number of iterations t f . The threshold parameter K governs the rate of diffusion: high values of K yield isotropic diffusion, typically over-smoothing the reflections; conversely, low values of K slow down the diffusion process. The number of iterations t f determines the number of neighboring pixels involved in the diffusion process. It is found that values of K = 10 and t f = 300 are effective for the PIV images used in this work. Readers are advised to perform the parametric study for a pair of images to determine the suitable values of K and t f . An automatic estimation of the optimal value of these parameters goes beyond the aim of this investigation, and is left for future work.
The proposed approach is applied to real PIV images acquired for the study of the blade vortex interaction, characterized by sharp and unsteady reflections of the propeller blades. Due to the unsteady character of the reflections, background removal approaches based on statistical analysis of the entire sequence of images are not effective. The results of the anisotropic diffusion background removal are compared with the conventional pre-processing methods of isotropic diffusion (sliding-average) filter, median filter subtraction and median-based-normalization filter. The comparison shows that the proposed approach is effective in removing the unsteady reflections, allowing the estimation of the particle's displacement, even in close proximity to the reflection region. In regions of the image not affected by reflections, the use of the anisotropic diffusion filter retains approximately the same image quality as in the raw images.
In the present work the performance of the method has been demonstrated for the case of sharp reflections, occurring e.g. when the laser light impinges on a solid surface. In the presence of diffused reflections, the anisotropic diffusion coefficient assumes approximately the same value in all directions, and the anisotropic filter behaves in practice as an isotropic filter (sliding-average).