Dynamic masking of PIV images using the Radon transform in free surface flows

Abstract

Particle image velocimetry (PIV) processing of free surface flow images often requires the use of digital masks to overcome the problems caused by the interface. In cases where a large number of particle images are collected it is essential that the time-varying boundary between the two phases can be tracked automatically to produce the binary masks. The Radon transform-based technique presented in this paper allows the automatic detection of the air–water interface in a stream of particle images acquired from a single camera. It is applied to time-resolved PIV measurements in the liquid phase of a stratified multiphase flow in a circular pipe. Accuracy estimations are provided using synthetic and real wave profiles. An extension to the more complex case of an overturning wave is also discussed.

1 Introduction

The presence of a wavy interface in particle image velocimetry images of free surface flows causes numerous issues. In most cases, the illuminating light sheet is oriented perpendicular to the air–water interface where undesirable reflections occur in the form of bright spots concentrated at the wave crests or troughs. Additionally, investigation of the liquid flow features requires the camera to be placed beneath the interface, which can act as a mirror, causing the particle images to appear symmetric on both sides. If their presence is not accounted for the ghost particles situated above, the interface will interfere in the PIV processing and result in a bias in the displacement estimation.

Lin and Perlin (1998) showed that these issues could be overcome by placing the camera at a certain upward inclination, defined as the Brewster angle, which eliminates all reflections from the water-to-air interface. However, such a solution demands some modification to the observation window and adds severe distortion to the images.

A more practical approach consists in minimizing the effects arising from interfaces in the processing of the images. In Theunissen et al. (2008), different methods are presented at the image preprocessing, correlation and data postprocessing stages of the PIV procedure to account for the presence of solid, fixed wall boundaries. Honkanen and Nobach (2005) proposed a background extraction method based on the assumption that noise sources can be discriminated from the particle images by their relative displacement, a condition which is not fulfilled in free surface flows. Digital masks are binary matrices containing ones at the pixel locations of interest and zeros elsewhere. They can be used directly to modify the images, in the present case by padding the area situated above the interface with zeros, or embedded in the correlation scheme (Gui et al. 2003). In situations where the boundaries delimiting the region of interest are fixed, the mask can be defined manually. However, in free surface flows, an automated identification procedure of the interface is needed.

Image-based tracking of surface wave profiles has been the focus of several recent studies, all of which rely on the detection of the vertical intensity gradient across the interface. In Zarruk (2005), a single camera positioned at a neutral angle was used for both PIV measurements and interface detection, leading to ambiguities in the location of the interface. Detection of the maximum intensity gradient was performed for each column of pixels, and outliers were removed using least-square regression and a statistical criteria. Mukto et al. (2007) recognized that a second camera was needed, placed above the interface at a downward angle so that the interface appeared as a clear boundary between a low-intensity region (the air phase) and a middle or high-intensity region (the liquid phase). An improved intensity threshold method accounting for inhomogeneous lighting was used to detect this edge and subsequent morphological operations yielded an accurate representation of the wave profiles. Their method was able to detect short wavelength capillary waves riding on top of the larger wind waves.

A similar idea was exploited by Hwung et al. (2009): they presented large-scale measurements of water wave profiles in a flume with a CCD camera placed above the interface in natural light. The effect of large overhead posture angles ranging from 15° to 60° was reported. In the three previous studies, the detection algorithms rely on the assumption that a single distinct edge corresponding to the interface location exists per column of pixels, i.e. that the wave is not overturning. Yao and Wu (2005) overcame this limitation with the use of an active contour model known as gradient vector flow (GVF) snakes. However, an edge map still had to be generated to initiate the process. Applications included a plunging breaker where the GVF snake successfully resolves the triple-valued water surface at the tip of the plunger.

The PIV image shown in Fig. 1 is part of co-current air–water flow experiments conducted at the Hydrodynamic Laboratory of the University of Oslo in a circular pipe. For specific values of the air and water flow rates, a wavy stratified flow regime was achieved where irregular, large-amplitude waves traveled down the length of the pipe (Andritsos and Hanratty 1987). Details about these experiments are provided in Sect. 2 and Sanchis et al. (2010). The intersection between the interface and the vertical laser light sheet appears as a discontinuous suite of bright spots caused by the light reflecting on the seeding particles floating at the surface. Above this line, ghost particles are visible due to the camera angle being inferior to the Brewster angle of the water-to-air interface. The upper blurry wave contour corresponds to the contact line between the surface and the curved pipe wall. Because of the pipe geometry and the 3-D character of the waves developing at the interface, the water elevation at the pipe walls differs significantly in phase and amplitude from the water elevation in the light sheet.

Fig. 1

PIV image of the liquid phase in air–water multiphase flow experiments

The sequence of bright particles displayed in Fig. 1 provides an unambiguous indication of the interface position in the PIV images, since these are situated at the intersection of the laser light sheet with the air–water interface. However, edge or contour detection methods would prove inadequate in this situation because the target particles are not easily distinguished from the other PIV particles, including the ghost particles situated above the interface.

A different method is presented in this paper based on the Radon transform, a mathematical tool well suited to the detection of linear features in noisy images and commonly applied in computer vision or seismic imaging (Murphy 1986). Section 2 provides information about the experimental set-up used to acquire the particle images. Section 3 describes the main properties of the Radon transform and details its implementation in the interface tracking algorithm. Section 4 discusses the accuracy and optimization of the method using simulated and real wave profiles. In Sect. 5, results are presented for the multiphase flow experiments presented previously, and an extension is proposed to the more complex case of an overturning wave. The main features of the Radon transform-based interface tracking (RTID) technique are outlined in Sect. 6.

2 Experimental set-up

The multiphase flow experiments were conducted in a horizontal 31 m PVC pipe with an internal diameter D = 10 cm. Pressure drop measurements with single-phase gas flow indicated that the relative roughness $\epsilon/D$ was in the range 1–2 × 10⁻⁴. The pipe consisted of adjacent sections with a length of 3.5 m connected by annular joints that ensured tightness. All joints were rigidly attached by collars to vertical beams that supported the whole structure. The vertical position of each joint was adjusted with care using an optical automatic level.

The fluids were air and water at atmospheric pressure with an average temperature of 21°C. Figure 2 shows the disposition of the pipe elements. Water was injected at the pipe bottom through a 5 cm I.D. tee branch. Honeycomb flow straighteners were placed right before and after the contact point between the liquid and gas phases. Depending on the liquid flow rate, a hydraulic jump formed at a distance varying between 1 and 10D downstream of the inlet. At the outlet, the pipe discharged into a separating tank at atmospheric pressure through a flexible plastic duct. Water and air were recirculated from the bottom and top exits of the tank, respectively. Twenty pipe diameters upstream from the outlet, the mean liquid level started to decrease in the pipe, as indicated in Fig. 2. Between this point and the hydraulic jump, the mean liquid level remained constant except at very low liquid flow rates.

Fig. 2

Schematic view of the circular pipe including details about the inlet and outlet conditions. The rectangular striped areas at the inlet represent honeycomb flow straighteners. A flexible duct connects the pipe outlet to the discharge tank

Water was circulated with a 1.4 kW centrifugal pump with a maximal volumetric flow rate of 90 m³/h. A frequency-regulated fan produced the airflow. The water and air mass flow rates were measured with a Endress Hauser Promass and an Emerson Micro Motion Coriolis flowmeters, respectively. Bulk velocities were calculated using a density of 997 kg/m³ for water and 1.2 kg/m³ for air, corrected with the temperature at the measurement point. Reynolds number in the gas phase based on the hydraulic diameter ranged from 4.8 × 10³ to 46 × 10³.

PIV velocity measurements in a vertical plane were performed in a channel section located 260D downstream from the pipe inlet and 50D upstream from the outlet. A Nd:YLF Darwin-Duo pulsed laser of 20 mJ provided the illumination, and images were recorded with a Photron FASTCAM-APX CMOS camera at a rate of 1,000 frames per second. The camera was placed 20 cm under the pipe centerline at an upward angle of 15° in order to map the liquid flow phase. Distortion caused by the curved pipe wall was minimized by placing the pipe measurement section inside a hollow square box filled with Isopar$\texttrademark$.

In our set-up, only one camera was used for both the subsurface velocimetry and interface detection purposes, contrary to Mukto et al. (2007). This choice was guided in part by the presence of a large number of water droplets at the higher gas flow rates investigated. Preliminary attempts to detect the interface from above the water level showed that these were picked up by the edge detection algorithm and disrupted the measurements. Also, at sufficiently high gas flow rates, the liquid cross-section in the pipe became crescent shaped with the waves washing along the sides of the pipe (Spedding and Spence 1993). A permanent film of water formed up the pipe walls and blurred the view to the pipe centerline. In these experiments, optical access to the air–water interface proved therefore only possible from under the water surface, i.e. through the liquid phase.

3 Radon transform theory and implementation

The Radon transform of a greyscale image is the projection of the image’s intensity along specified directions using line integrals (Bracewell 1995). Denoting f(x, y) as the image intensity function, its Radon transform $\mathcal{R}$ is defined by

$$ {\mathcal{R}}(\theta,\rho)=\int\int f(x,y)\delta(\rho-x \cos \theta-y \sin \theta) \hbox{d}x \hbox{d}y, $$

(1)

where δ is the Dirac delta function. The second term in the integrand forces the integration of f(x, y) along the line

$$ \rho-x \cos \theta-y \sin \theta=0 $$

(2)

as shown in Fig. 3, θ being the orientation and ρ the offset relative to a parallel line passing through the image centre. A simple example is provided in Fig. 4: the Radon transform to the right shows a maximum at ρ = −3 pixels and θ = 45° corresponding to the alignment between the two bright dots in the input image. The radon.m Matlab function available in the Image Processing Toolbox (Matlab$\textregistered$ R2010a, The MathWorks Inc., Natick, MA, USA) was used to compute the Radon transform values. A normalization procedure was added to take into account the aspect ratio of rectangular input images.

Fig. 3

Geometry of the Radon transform

Fig. 4

Application of the Radon transform to a simple case. Left input image; right (θ, ρ) plane

The Radon transform is well suited to the detection of linear features in noisy images, since intensity fluctuations due to noise tend to be canceled out by the process of integration. However, one limitation is that the length of the linear feature of interest should approach or exceed the original image dimension, otherwise the associated Radon transform value would decrease sharply due to the fact that line integrals are performed across the whole image.

The general approach followed in the interface detection problem outlined in Sect. 1 is to divide the original images into N _s vertical segments and apply the Radon transform on each to detect the local alignment of the bright particles lying at the interface, hereafter referred to as the “target particles”. The segment length is constrained by two opposing criteria. It should be sufficiently short so that the target particles appear aligned along a straight line, leading to a sharp peak in the (θ, ρ) plane. Additionally, each segment should contain enough of the target particles so that their alignment stands out against the random alignments of the background particles. The second condition is particularly important if there are varying illuminating conditions across the image and because the target particles are placed at irregular intervals along the interface, augmenting the risk of having one segment containing too few target particles, or even none. The significance of N _s in terms of accuracy is discussed in detail in Sect. 4.1.

Figure 5a shows a 100 pixels wide segment of a PIV image similar to the one displayed in Fig. 1. Because the segment width corresponds to only a small fraction of the wavelength, the target particles appear to be aligned along a straight line that provides a reasonable estimate of the interface location. In the Radon transform shown in Fig. 6a, the brightest peak at ρ = 50 pixels and θ = +4° corresponds to the upper blurry feature marking the contact line between the interface and the near pipe wall. The alignment of the target particles is represented by a secondary peak at ρ = 18 pixels and θ = +6°. A number of local maxima also appear due to the random alignment of the seeding particles. In this case, the range of θ has been limited to ±20° in order to reduce the computing cost, based on the maximum expected wave slope.

Fig. 5

a a 100 pixels wide segment of Fig. 1; b the position and slope of the interface found by the RTID method

Fig. 6

a Radon transform of the segment displayed in Figure 5; b enhanced image after histogram truncation and renormalization

Selection of the local maximum that corresponds to the interface location is conducted in two steps: first, the histogram of the Radon transform image is truncated and re-normalized (Fig. 6b) in order to eliminate the noise produced by the seeding particles. It is crucial at this point that enough of the target particles are present in the PIV image segment; otherwise, this procedure would remove their contribution. Secondly, of the remaining peaks, the closest to the image bottom boundary is selected as the one corresponding to the target particles, an assumption that can be made based on the camera position and viewing angle.

At each segment, the algorithm yields a pixel location and an inclination angle that define the line in Eq. 2, as shown in Fig. 5b. Reconstruction of the interface is done through piecewise cubic Hermite interpolation bewteen adjacent segments, using the pixel location and slope of the interface in each segment as input. In Fig. 13a, the detected interface is superposed to the PIV image shown in Fig. 1. Ten images segments were used, and the corresponding nodes are displayed as white circles.

The computational steps composing the RTID algorithm are summarized in Fig. 7. Sorting the nodes is a trivial task when the interface position corresponds to a single value in each segment of the image. Likewise, if the segment size is chosen carefully, it is unlikely that the interpolation procedure should lead to singularities such as self-intersections of the interface. Steps 4 and 5 were thus simply bypassed in the processing of the multiphase flow images shown in Sect. 5.1. However, they are essential to the analysis of more complex wave profile geometry such as the overturning wave presented in Sect. 5.2.

Fig. 7

Flowchart describing the computational steps used in the RTID method

4 Error analysis and limitations

4.1 Influence of N _s on the random error level

The accuracy of the present method was assessed using simulated wave profiles as in Mukto et al. (2007). Those were defined by:

$$ \eta_i(x)=a \sin\left(\frac{2 \pi}{\lambda} x+ \phi_i\right) \quad \hbox{with} \; i=1.400, $$

(3)

where the amplitude a and wavelength λ were set to values representative of the waves captured in the multiphase flow experiments described in Sect. 5.1, respectively 40 pixels (1 cm) and 1,500 pixels (30 cm), while the phase ϕ varied randomly. Firstly, each profile was represented by a continuous line of white pixels in a black background image 120 pixels high and 1,024 pixels long, as shown in Fig. 8a. The positioning of the pixels along the wave profiles defined by Eq. 3 is subject to round-off inaccuracies that were found to generate a r.m.s error level of 0.28 pixels. The linear feature detection method was run on the 400 test images using a number of image segments N _s varying between 4 and 14. For each value of N _s , the r.m.s difference $\epsilon_{rms}$ between the simulated wave profiles η _i and the detected wave profiles y _i was calculated according to Eq. 4:

$$ \epsilon_{\hbox {rms}}=\frac{1}{400}\sum_i\sqrt{\sum_x(y_i(x)-\eta_i(x))^2}-0.28 $$

(4)

where the round-off error has been substracted in order to estimate the random error imputable to the RTID algorithm only.

Fig. 8

A 500 pixels wide portion of the simulated wave profile images used as input for the error analysis. a Continuous line of white pixels, b randomly scattered white dots and c 5 × 5 Gaussian particles and PIV image background

Additionally, we noted the fraction of samples where a maximum absolute error of more than 2 pixels was detected. This threshold was set arbitrarily to determine whether the method failed to detect the simulated wave profile in one or more segments of a sample. In the real images, ±2 pixels also corresponds to the uncertainty in the position of the interface due to the overexposure of the target particles.

As shown on Fig. 9, $\epsilon_{\hbox {rms}}$ decreases monotically as the number of segments N _s , and hence the spatial resolution, increases. For N _s  = 4, the 1,024 pixels wide image is divided into four vertical segments, yielding one data point every 256 pixels which represents 6 points for every wave period. At N _s  = 14, the resolution has increased to 20 points per period and the r.m.s difference converges towards zero after subtracting the round-off error. Likewise, the fraction of samples where a maximum absolute error of more than 2 pixels was recorded dropped to zero for N _s  ≥ 5.

Fig. 9

R.m.s difference (excluding the round-off error) between the detected and simulated wave profiles as a function of the number of segments N _s . Open squares continuous white line; filled circles scattered white pixels; and open circles scattered 5 × 5 gaussian particles over PIV image background

A second test was performed where only 80 white pixels were spread randomly along the wave profiles on the same black background images, yielding an average horizontal spacing of 13 pixels between adjacent white dots (Fig. 8b). This time, increasing the number of segments not only enhanced the spatial resolution but also reduced the average number of white dots contained into each segment, resulting in a less precise estimation of the interface position and augmenting the probability that a segment contained no points at all.

This is reflected in Figs. 9 and 10 where an optimum number of segments can be found balancing spatial resolution and signal loss. For N _s  = 7, only one sample out of 400 presents a maximum absolute error greater than 2 pixels and the r.m.s difference (excluding round-off inaccuracies) is 0.16 pixels.

Fig. 10

Fraction of samples presenting a maximum absolute error of more than two pixels as a function of the number of segments in simulated wave profiles tests. Refer to Fig. 9 for symbols

In a third test, the black background was replaced by a portion of PIV particle images taken in the liquid phase (Fig. 8c). Synthetic particles replaced the white dots used in the second test. Their intensity followed a 5 × 5 gaussian pattern with a standard deviation of 2 and peak of 255, simulating overexposed particles similar to those marking the interface in the real images. Due to the additional noise of the background, the performance of the method degraded faster as the number of segments increased and a minimum in the two error criteria defined above was observed for N _s  = 6, with 3% of the samples presenting maximum absolute errors greater than 2 pixels and $\epsilon_{rms}=0.39$ pixels. However, the interface is usually better defined in the real images, as seen in Fig. 1, than in Fig. 8c. The spreading of the synthetic particles in the latter case is representative of areas of the PIV images containing less particles at the surface and/or with insufficient illumination. The error levels presented in Fig. 9 can therefore be regarded as conservative.

4.2 Systematic (bias) error in the wave crest region

A bias error occurs in the crest region of the waves where the curvature of the water surface is maximum and the interface may no longer be approximated as a suite of adjacent straight segments. Unless the number of segments N _s is increased, which in turn would lead to the errors discussed in Sect. 4.1, the RTID method is bound to systematically underestimate the interface position at the waves’ crests. To assess the magnitude of this offset, fourteen images containing some of the highest waves found in our data set were selected. The “exact” location of the air–water interface was determined by manually picking a number of points using a graphical input function and fitting a smoothing spline. The wave profiles were then detected with a number of segments N _s varying between 4 and 14. The mean bias error between the “exact” and detected wave profiles in the whole image and the vertical offset under the waves’ crests were computed for each image and then averaged over the fourteen samples. Results are presented in Fig. 11. For N _s  ≥ 7, the mean bias error was found to stabilize at −0.25 pixels, but the vertical offset under the wave crests continued to decrease to reach −1.1 pixels for N _s  = 14. The offset vanished for all waves at N _s  = 20, however, using such a high value for a large number of images would compromise the robustness of the method. A single case is illustrated in Fig. 12 for a wave of amplitude 30 pixels (6.2 mm) and front slope 15°.

Fig. 11

Bias error between the “exact” and detected wave profiles as a function of the number of segments N _s . Open symbols mean value over the whole image, filled symbols vertical offset at the wave crest

Fig. 12

The RTID algorithm applied to a single crest. Dashed line N _s  = 6, continuous line N _s  = 20. The PIV image is displayed in the background. The y-axis has been stretched with respect to the x-axis with a ratio of 3:1 for clarity

Finally, it should be noted that due to the inherent limitation on the number of nodes at which information about the slope and height of the interface is collected, the RTID method is not suited to the detection of small-scale disturbances such as capillary waves. These were beyond the scope of the study for which the present algorithm was designed.

5 Applications

5.1 Multiphase flow PIV images

The set-up described in Sect. 2 was applied to the study of interfacial waves developing in a two-phase, stratified flow regime. Disturbances formed at the pipe inlet and evolved into large amplitude waves. At a sufficient distance along the pipe, a balance was reached between the stabilizing effect of gravity and the destabilizing effect of air pressure variations and a fully developed flow field was obtained.

PIV images of the liquid phase were acquired in a test section situated 260 pipe diameters from the inlet (see Fig. 2) at constant values of the water, and air mass flow rates Q _l and Q _g . 600 image pairs were collected for a given value of Q _l and Q _g . The RTID method was used to automatically detect the position of the interface in those images. Examples are shown in Fig. 13 for two sets of flow rates. In the latter case, the higher air mass flow rate lead to a decrease in the mean liquid level and steeper wave slopes. At low and moderate values of Q _g , the interface was detected successfully for all images in each data set. The number of segments N _s was adjusted to the different cases, with shorter and steeper waves requiring more segments to be properly represented. For the higher air mass flow rates, measurements were affected by the entrainment of air into the liquid phase caused by roll waves or slugs. The presence of air bubbles in the laser light sheet produced strong reflections, which prevented the alignment of the target particles from being detected.

Fig. 13

The RTID method applied to multiphase flow PIV images in a 10 cm diameter circular pipe. a $Q_l=1.1\,\hbox{kg}/\hbox{s},\; Q_l=0.0107\,\hbox{kg}/\hbox{s}$; b $Q_l=1.1\,\hbox{kg}/\hbox{s},\;Q_l=0.0133\,\hbox{kg}/\hbox{s}$

The position of the interface was used to generate digital masks containing zeros at all pixel locations situated above the interface and ones underneath. In the subwindows enclosing part of the interface, the masked out pixels were excluded by the DigiFlow PIV processing software ($\copyright$Dalziel Research Partners) to avoid them influencing results in the displacement estimation process. As an example, Fig. 14 displays an instantaneous velocity vector field obtained in a 10-cm liquid layer underneath the interface, with the original particle image shown on top.

Fig. 14

Instantaneous velocity vector field in physical coordinates with the interface detected by the RTID method shown in grey. The y-axis is dilated for clarity. Top original particle image with a data aspect ratio of 1

The RTID algorithm performed reasonably fast: with N _s  = 10 and θ varying between −20° and +20° by increments of 0.5° the interface detection process applied on a 1,024 × 300 pixels portion of an image took 1.32 s on a single Intel$\textregistered$ Core$\texttrademark$ E3600 (1.86 GHz) processor.

5.2 Overturning wave

In this section, the RTID method is applied to the more complex case of an overturning wave. The image displayed in Fig. 15a was acquired from an experimental study on breaking waves on a beach. Similar incident waves are reported in Jensen et al. (2003). As in Fig. 1, the interface is characterized by the alignment of bright particles at the intersection with the vertical laser light sheet; however, a number of additional difficulties arise due to the wave geometry:

• a higher number of image segments N _s is needed in order to correctly represent the wave curvature. The range of possible angles of inclination θ also has to be extended to [−90 90] to cope with nearly vertical portions of the interface.

• the interface is hardly visible in the wave trough where lighting conditions are poor.

• for each vertical segment of the image either one or three values of ρ and θ can be expected to correspond to the interface position. The task of reconstructing the interface is therefore complicated.

Fig. 15

The RTID method applied to an overturning wave. a Original image; b initial solution with the nodes shown in white circles; and c final solution after filtering

As indicated on the flowchart displayed in Fig. 7, the nodes are first connected to each other based on a nearest neighbour principle, starting in this case with the lower left corner of the image. The result is shown in Fig. 15b with N _s  = 35 segments. In the wave trough where the interface cannot be distinguished, the method failed to find any valid point. A number of spurious nodes corresponding to the random alignment of seeding particles were detected in this area and included in the solution. Fortunately, the angle of inclination associated with these nodes is uncorrelated with the angle of inclination of adjacent nodes so that the interpolation process leads to aberrations such as double point singularities or excessive curvature.

A filtering procedure (step 5 in Fig. 7) is then initiated whereby the interface is being represented as a plane parametric curve x(t), y(t) whose curvature is defined by

$$ \kappa(t)=\frac{|x^{\prime}y^{\prime\prime}-y^{\prime}x^{\prime\prime}|} {(x^{\prime\;2}+y^{\prime\;2})^{3/2}}. $$

(5)

Spurious points are identified as being closest to self-intersections, cusp (hairpin) singularities and/or areas where κ exceeds some preset threshold value. Only one point is removed, and the remaining nodes are reorganized into a new sequence from which the interface is reinterpolated. The procedure is repeated until the regularity criterias described previously are satisfied. Figure 15c shows the final interface position after four passes with the control points displayed as white circles. In the wave trough, the filtering procedure succeeded in eliminating spurious data points; however, bad lighting conditions prevented the interface from being detected at all, and the nearest valid points are too far to provide a satisfactory interpolated solution. The rest of the wave profile is correctly described including the wave tip where curvature is highest.

Tests were conducted on other overturning wave images from the data set of Jensen et al. (2003) with similar results. Contrary to the multiphase flow case presented in Sect. 5.1, it was not possible to find a unique value of N _s for which all images could be successfully treated. Due to the complex, rapidly evolving shape of the wave and nonuniform light conditions, each image had to be analysed separately. The processing effort also increased dramatically compared to the multiphase flow images, with each sample requiring more than 13 s on a single processor.

6 Conclusion

In this paper, the automatic detection of an air–water interface in free surface flow PIV images is addressed. The present method relies on the use of the Radon transform to detect the alignment of seeding particles floating at the interface and illuminated by the light sheet. For this purpose, the original PIV image has to be decomposed into a number of segments where the local interface position and inclination is found. Hermite cubic interpolation between adjacent segments permits the reconstruction of the interface across the whole image. A number of tests performed on synthetic wave profiles indicate that the technique is able to track the interface with an accuracy of ±0.67 pixels in the worst case (low density of target particles and noisy background), with 3% only of the sample images presenting a deviation of more than 2 pixels from the solution. Error analysis also showed that the number of segments N _s chosen to decompose the PIV image could be optimized depending on the concentration of target particles and the wave field characteristics.

The RTID method was initially developed to automatically generate digital masks for sub-surface, high-speed PIV velocimetry in multiphase flow experiments. A single camera was placed below the pipe centerline at an upward angle to acquire the PIV images. Data sets of 600 images were recorded at different values of the air and water mass flow rate yielding a stratified wavy flow regime. For each set, the method managed to detect correctly the 600 interfaces in one pass without any user input beyond the initial choice of N _s .

Modifications were suggested to apply the RTID technique to the more complex case of an overturning wave. The original algorithm was adapted to treat the triple-value problem caused by the presence of the interface at multiple locations in a single vertical segment. A filtering procedure was designed to eliminate data points created by the random aligments of seeding particles in the liquid phase. Except in regions of the image where lighting is insufficient, the overturning wave profile was detected successfully.