# Automatic dynamic mask extraction for PIV images containing an unsteady interface, bubbles, and a moving structure

David Dussol, Philippe Druault, Bachar Mallat, Sylvain Delacroix, Grégory Germain

► To cite this version:

David Dussol, Philippe Druault, Bachar Mallat, Sylvain Delacroix, Grégory Germain. Automatic dynamic mask extraction for PIV images containing an unsteady interface, bubbles, and a moving structure. Comptes Rendus Mécanique, 2016, 344 (7), pp.464-478. 10.1016/j.crme.2016.03.005 . hal-01304317

HAL Id: hal-01304317

https://hal.sorbonne-universite.fr/hal-01304317

Submitted on 19 Apr 2016

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Automatic dynamic mask extraction for PIV images containing an unsteady interface, bubbles, and a moving structure

David Dussol<sup>a</sup>, Philippe Druault<sup>a,*</sup>, Bachar Mallat<sup>b</sup>, Sylvain Delacroix<sup>b</sup>, Grégory Germain<sup>b</sup>

## ARTICLE INFO

**Article history:**

Received 16 September 2015

Accepted 21 March 2016

Available online xxxx

**Keywords:**

Particle image velocimetry

Unsteady free interface

Edge detection method

Morphological analysis

## ABSTRACT

When performing Particle Image Velocimetry (PIV) measurements in complex fluid flows with moving interfaces and a two-phase flow, it is necessary to develop a mask to remove non-physical measurements. This is the case when studying, for example, the complex bubble sweep-down phenomenon observed in oceanographic research vessels. Indeed, in such a configuration, the presence of an unsteady free surface, of a solid-liquid interface and of bubbles in the PIV frame, leads to generate numerous laser reflections and therefore spurious velocity vectors. In this note, an image masking process is developed to successively identify the boundaries of the ship and the free surface interface. As the presence of the solid hull surface induces laser reflections, the hull edge contours are simply detected in the first PIV frame and dynamically estimated for consecutive ones. As for the unsteady surface determination, a specific process is implemented like the following: i) the edge detection of the gradient magnitude in the PIV frame, ii) the extraction of the particles by filtering high-intensity large areas related to the bubbles and/or hull reflections, iii) the extraction of the rough region containing these particles and their reflections, iv) the removal of these reflections. The unsteady surface is finally obtained with a fifth-order polynomial interpolation. The resulted free surface is successfully validated from the Fourier analysis and by visualizing selected PIV images containing numerous spurious high intensity areas. This paper demonstrates how this data analysis process leads to PIV images database without reflections and an automatic detection of both the free surface and the rigid body. An application of this new mask is finally detailed, allowing a preliminary analysis of the hydrodynamic flow.

© 2016 Académie des sciences. Published by Elsevier Masson SAS. This is an open access article under the CC BY-NC-ND license

(http://creativecommons.org/licenses/by-nc-nd/4.0/).

## 1. Introduction

When performing PIV measurements in complex fluid flows with moving interfaces and a two-phase flow, the induced laser reflections may cause a lot of erroneous velocity vectors. For instance, the investigation of the bubble sweep-down phenomenon [1] observed on oceanographic research vessels [2–4] based on PIV measurements remains today a great challenging task. Indeed, in such a flow configuration, i) the presence of the moving hull surface in PIV image, ii) the presence

<footnote><sup>a</sup> Sorbonne Universités, UPMC Université Paris-6, CNRS, UMR 7190, Institut Jean-Le-Rond-d'Alembert, 75005 Paris, France</footnote>

<footnote><sup>b</sup> IFREMER, Marine Structures Laboratory, 150 quai Gambetta, 62321 Boulogne-sur-Mer, France</footnote>

<footnote>* Corresponding author.</footnote>

<footnote>E-mail address: philippe.druault@upmc.fr (P. Druault).</footnote>
of a free water-air interface in PIV image, and iii) the presence of bubbles in the measurement PIV plane and on the outside of the illuminated measurement plane cause several problems impeding a better exploitation of these measurements. These two surfaces act as a semi-reflective mirror inducing spurious laser reflections and ghost particles. Thus, before carrying out the analysis of the hydrodynamic flow, we need to solve these issues by automatically identifying both the area above the air-water interface and the hull area for each PIV frame. Instead of doing it by hand [3,5], one then develops a dynamic masking algorithm to isolate the hull structure and the free interface, which changes for each PIV frame. Previous research was mainly based on the intensity-based methods [6]. For example, some edge contour method operator like the one proposed by Canny [7] has been used as an image-processing method for capturing the surface [8]. However, in the presence of a flow with complex interfaces (deformable, two-phase, etc.), this edge detector method may not provide satisfactory results [9].

More specifically, for stationary solid fixed wall surfaces, the mask can be automatically applied to each PIV frame. However, near the wall and especially for high-velocity flows, the PIV post-processing treatment needs to be carefully performed to avoid erroneous PIV vector fields. For instance, Theunissen et al. [10] proposed an innovative interrogation method that allows one to increase the robustness and to enhance the spatial resolution.

For a moving solid object, Jeon and Sung [11] detailed a surface detection algorithm for a flow developing around a moving surface in which small texture patterns called textons are used to identify the surface. Such an algorithm may not be appropriate in the presence of a high deformable interface requiring some advanced PIV preprocessing method. To solve this issue, Park et al. [12] proposed to replace a classical PIV interrogation window with a window deformation adapted to the deformable surface.

To isolate a wavy free surface in a two-phase flow, previous works have proposed to analyze the mirror images of particles close to the air-water interface [13,14]. For instance, Sanchis and Jensen [13] used the Radon Transform for such an identification. Note that these methods cannot be applied to current flow measurements due to the presence of bubbles and their associated mirror particles. Another solution for extracting a wavy surface consists in using several cameras [9,15].

The use of a dynamic mask has also been applied to PIV measurements in the presence of bubbles. Different phase discrimination methodologies based on different concepts (edge detection [4,16], image intensity, filtering, size and intensity distribution, etc.) (see [17] and references therein) have been proposed. For instance, Cheng et al. [17] presented a filtering technique allowing one to clearly distinguish the two phases based on the associated sizes and intensity distributions.

In this paper, we propose a new dynamic masking algorithm applied to PIV measurements with unsteady free surface, with oscillatory hull motion and bubbles. It relies on the determination of a geometric mask based on an analysis of each individual raw PIV frame. Then, each mask is automatically applied to the corresponding PIV vector field. We do not modify locally the PIV image interrogation window to access the velocity vector near the boundaries.

The paper is organized as follows. The experimental set-up including PIV measurements is presented in section 2. Then, the algorithm allowing the detection of the hull oscillation surface in each PIV frame is given in section 3. Finally, in section 4, an original algorithm for the extraction of the free water-air interface from each PIV frame is provided.

## 2. Description of the experimental database

### 2.1. Experimental set-up

Experiments have been carried out in the IFREMER wave and current circulating tank. The flume tank is 18 m long, 4 m wide, and 2 m deep. A side observation window of 8 × 2 m² placed on one side of the tank allows users to observe the behaviour of the ship model during trials. A detailed description of the experimental apparatus can be found in [18]. A schematic representation of the experimental set-up is given in Fig. 1.

To reproduce the bubble sweep-down phenomenon observed by the French oceanographic vessel *Pourquoi Pas ?* during an acoustic survey performed in 2005 [2], a 1/30 scale ship model has been considered. Its length is 3.13 m between perpendiculars with a beam of 0.67 m and a draft of 0.183 m (see Fig. 2, left). In order to reproduce the sea conditions encountered during the sea survey, the trials have been carried out with a flow velocity of $U_\infty = 0.75 \text{ m s}^{-1}$ and regular waves for synchronization purposes. The wave amplitude and frequency are: $A = 33 \text{ mm}$ and $f_0 = 0.85 \text{ Hz}$. Details about the wave generator system (see Fig. 2, right) can be found in [19].

### 2.2. Description of the oscillatory hull motion

A 6 degrees of freedom motion generator is used to reproduce the hull motions [3,20]. In this work, the imposed motions are determined by means of a numerical model using computer program Fredyn [21]. In this program, a nonlinear strip theory approach is used to compute the hydrodynamic forces acting on the hull. In the time domain simulation, the ship motion components are determined from a set of six coupled differential equations of motion. The integration of the calculated accelerations provides the ship's, velocities, and a second integration leads to the ship's position at each time step. In this head sea configuration, the only significant motions are the heave and the pitch, whereas surge motions are insignificant (< 0.1 m at full scale). After converting the results to the present hull model scale by Froude scaling, we finally obtain 2° for the pitch and 20 mm for the heave. In the measurement plane, the imposed hull motions with a frequency of 0.85 Hz have been synchronized with the incoming waves.
<figure><img src="image_2.png" /><figcaption>Fig. 1. Top and side views of the schematic representation of the experimental set-up in a wave and current circulating tank.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Presentation of the 1/30 ship model with the ship motion generator system (left) and the wave generator system (right) generating regular waves with a current of 0.75 m s⁻¹.</figcaption></figure>

## 2.3. PIV measurements

Two-dimensional time-resolved PIV measurements have been carried out. Illumination has been provided by a standard, frequency-doubled, double-cavity Nd:YAG laser (NewWave, Gemini PIV), with a pulse energy of up to 120 mJ per pulse at 532 nm. An optical system mainly composed of cylindrical lenses has been used to generate a vertical light plane between two successive times with a very low interval time $t = 1200$ µs, which also corresponds to the time between the double frame images. Images have been treated with 8-bit image encoding. The seeding particles (referenced as Dantec S-HGS-10) used for the experiments are silver-coated hollow glass spheres 10 µm in diameter. Image acquisition has been carried out using a high-resolution Dantec Hi-Sense CCD camera of 1600 × 1200 pixels². The FlowMap software from Dantec Dynamics has been used for PIV image processing.
<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Details of the PIV measurement planes locations. In the graph on the right, the two dotted lines represent schematically the highest and lowest positions of the free surface during measurements.</figcaption></figure>

<table><caption>**Table 1**
Main parameters of both experimental PIV databases under consideration.</caption><thead><tr><th>Parameters</th><th>Database 1</th><th>Database 2</th></tr></thead><tbody><tr><td>Interrogation window size overlapping</td><td>32 × 32<br>25%</td><td>32 × 32<br>25%</td></tr><tr><td>Physical dimension (PIV plane)</td><td>378.1 × 283.1 mm<sup>2</sup></td><td>528.1 × 396.1 mm<sup>2</sup></td></tr><tr><td>Number of grid points</td><td>(n<sub>x</sub>, n<sub>z</sub>) = (66,49)</td><td>(n<sub>x</sub>, n<sub>z</sub>) = (66,49)</td></tr><tr><td>Spatial discretization</td><td>dx = dz = 5.75 mm</td><td>dx = dz = 8 mm</td></tr><tr><td>Pixel resolution</td><td>1600 × 1200</td><td>1600 × 1200</td></tr><tr><td>Pixel discretization</td><td>dp<sub>x</sub> = dp<sub>z</sub> = 0.24 mm</td><td>dp<sub>x</sub> = dp<sub>z</sub> = 0.33 mm</td></tr><tr><td>Pixel and spatial discretization</td><td>dx ≈ 24dp<sub>x</sub></td><td>dx ≈ 24dp<sub>x</sub></td></tr><tr><td>Sampling frequency f<sub>e</sub> (Hz)</td><td>8.5</td><td>8.5</td></tr><tr><td>Number of PIV snapshots</td><td>1530</td><td>1800</td></tr></tbody></table>

<table><caption>**Table 2**
Description of each PIV database carried out with both PIV measurement parameters given in Table 1.</caption><thead><tr><th>Label</th><th>Location</th><th>Experimental conditions</th></tr></thead><tbody><tr><td>(a)</td><td>Symmetry plane (y = 0)</td><td>Swell with fixed hull</td></tr><tr><td>(b)</td><td>Plane (y = 100)</td><td>Swell with fixed hull</td></tr><tr><td>(c)</td><td>Symmetry plane (y = 0)</td><td>Swell with forced hull oscillations</td></tr><tr><td>(d)</td><td>Plane (y = 100)</td><td>Swell with forced hull oscillations</td></tr></tbody></table>

The PIV measurement system has been implemented by using two distinct parameters allowing an investigation of both the near flow field (database 1) and of the flow field in a more extended area (database 2). Database 2 is deduced from the first one by repositioning the camera and repeating the experiments, because of the field width required. Indeed, database 2 allows the visualization of both the bubble sweep-down phenomenon and the flow field in a large area including the free surface and the hull. Note that during the measurements related to database 1, the free surface may be above the measurement area (see Fig. 3).

The associated parameters are given in Table 1. Furthermore, for each database, we perform four PIV measurements campaigns based on different flow conditions and plane locations (see Fig. 3). A description of these conditions is provided in Table 2. In the following, each database will be referenced by using the label detailed in these tables. For instance, the PIV database denoted ‘2d’ corresponds to i) the PIV parameters mentioned in Table 1 for database 2 and ii) the flow conditions and the measurement plane location indicated in Table 2 with label (d).

To process a hydrodynamic analysis and for a future investigation of the bubble generation, the dynamic positions of the free surface and of the hull surface in each PIV image must be determined. In the following, we describe the algorithm allowing an instantaneous characterization of an adaptive dynamic mask applied to each PIV vector field. Thus, an instantaneous PIV vector field will be set to zero in the region within the mask hull structure and also over the detected free surface.
<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Schematic representation of the successive steps of the algorithm allowing the determination of the free surface.</figcaption></figure>

### 3. Detection of the hull edge contours in the PIV image

The removal of the hull surface in the PIV frame is accomplished after determining the edge contours of the hull in the frame. The hull's presence in the PIV image generates quite high-intensity areas related to laser reflections. Some illustrations are given in Fig. 5, where these reflections and the edge of the hull are clearly indicated. The edge contours of the hull are manually determined in the first PIV image of each measurement campaign.

When dealing with flow configurations without hull motion (database ‘a’ and ‘b’), this mask is duplicated for all other PIV images.

When dealing with flow configurations with hull motion (databases ‘c’ and ‘d’), knowing the imposed hull motion and the PIV measurement frequency, the mask related to the hull edge contours is directly deduced from the one related to the first PIV image. This dynamic mask is then applied to the rest of the PIV images in order to remove the region of the PIV image containing the hull.

An illustration of the latter is given in Fig. 13, where the dashed-red line represents the hull edges extracted from some selected PIV images.

### 4. Detection of the free surface

The PIV raw images obtained in grayscale intensity are first considered. The variations of the intensity observed in each raw PIV image are associated with PIV particles, bubbles (high-intensity areas) and with different types of background noise (reflections, particles or bubbles in or out of measurement plane, or other sources of noise). Some examples of PIV raw images are displayed in Fig. 5, one in the presence of bubbles and bubble-surface reflections, and the other one with only surface reflection. Indeed, the curved water-air interface acts as a semi-reflective mirror modifying locally the intensity in PIV images. Moreover, the presence of bubbles generated by breaking bow waves in PIV raw images located out of the measurement plane also induces some gradient magnitude. This increases the difficulty in determining accurately the free surface. Another difficulty encountered concerns the time-space deformable nature of the free surface.

As a first investigation, the classical Canny [7] or Prewitt [22] edge detection methods are considered to extract directly the free surface from the PIV image. The edge detection algorithms from the Image Processing toolbox in MATLAB are used.

Canny's method corresponds to a signal-processing optimization problem allowing the edge determination when the gradient magnitude rises to a maximum in the gradient direction based on two user-defined threshold values (to detect strong and weak edges). This method includes the weak edges in the output only if they are connected to strong edges.

Prewitt's method detects horizontal and vertical edges in an image by computing an approximation of the gradient of the image intensity magnitude. The application of this method requires to also select a user-defined threshold value (called ‘thresh’ in Matlab). It is related to a specific gradient magnitude, generating an output with only magnitudes higher than this threshold.

For the present application, different threshold values have been successively tested for both methods. Both provide quite similar results and fail to directly isolate the free surface, especially in the presence of numerous bubbles and reflections near the free surface, confirming previous investigations about the difficulty in retrieving well the edges of a PIV raw image on different scales [23]. As the classical surface detection methods are not robust enough for such an investigation [9], a new algorithm is developed, taking into account the complexity of the PIV images under consideration. This new procedure is detailed in the following, and the successive steps are schematically represented in Fig. 4, from the edge and particles detection to the localization of the free surface.
<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Step 1. Example of two raw PIV images for database ‘1b’ (a) and database ‘1a’ (c). Corresponding images after the application of the edge detection method (b) and (d).</figcaption></figure>

## 4.1. Step 1: Edge detection of gradient magnitude

### Step 1-1: Edge detection contour

As it is not possible to determine directly the location of the surface, we propose to extract first the particle candidates that include the bubbles, the conglomeration of PIV particles and unwanted bright noise or the surface hull reflections.

This first step consists in using the Canny [7] or Prewitt [22] edge detection methods to detect the high gradient magnitude in the image. Both methods provide quite similar results and, in the following, we only present results based on Prewitt's method. Different values of the threshold are successively tested for individual images. Finally, by looking at threshold automatically chosen by the algorithm applied to images with or without bubbles, it seems that this automatic threshold is able to discriminate edge contours in agreement with what is expected. Indeed, we expect to avoid any detection of too many individual pixels related to the particle or to very small reflections above and close to the free surface. The threshold's values vary from ≈ 0.03 (for PIV frames without bubbles) to ≈ 0.06 (for PIV image with bubbles). As Prewitt's method is based on the gradient of the image intensity magnitude, when one pixel corresponds to a high gradient magnitude, it extracts the nine pixels located around the latter. However, in the case of a low gradient magnitude, only one pixel can be recovered. At this edge detection stage, the gray-scale image is now binarized.

### Step 1-2: Morphological closing algorithm

Now that the image is binarized, the morphological closing algorithm (function *imclose* in MATLAB solver), which consists of morphological dilatation followed by erosion, is applied close to the edges of the detected white area. Morphological dilatation thickens the white area boundary by adding pixels in the neighborhood of the latter. Here, a disk-structuring element of 2 pixels radius is used to preserve the circular nature of the particle. Note that it is also possible to merge neighboring white area lines that might not be connected when using Prewitt's edge detection method. An illustration of the Prewitt's application followed by the application of the morphological procedure is displayed in Fig. 5 with a zoom view in Fig. 6. Note that these two figures are representative examples. The proposed free surface detection algorithm has to consider now not only the simplified flow configuration (like in Fig. 5a) but also the more complex one related to the bubble sweep down phenomenon (like in Fig. 5b). When dealing with PIV images without bubbles (Fig. 5d), the presence of hull reflections above the free surface is clearly indicated. In this case, the water-air interface seems to be easily detectable
<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Step 1-2. Database '1b'. Zoom views of the area delimited with green and black lines in Fig. 5 (top left). Raw PIV image ((a) green area and (c) black area) and corresponding image after the successive application of the Prewitt's method and of the morphological one ((b) green area and (d) black area).</figcaption></figure>

outside the hull. Conversely, for images containing bubbles and their reflections (see Fig. 5b), the determination of the interface is not an easy task, even by hand. However, this first step will facilitate the discrimination between the air and the water interface in the following.

## 4.2. Step 2: Detection of the particles

### Step 2-1: Extraction of particles, bubbles and reflections.

In a previous investigation of bubbles in a two-phase flow, Deen et al. [24] proposed a mask applied to PIV frames to extract the bubbles. This mask is based on a separate analysis of two frames: the first one, not containing any bubbles, is subtracted from the frame with bubbles to be analyzed and the resulted image is scaled. Such background subtraction allows us to remove both the noise and the bubble edges [24]. In a similar way, but without the same objective, we propose to extract the conglomerate of particles and the particles themselves. This permits to filter bubbles and their associated reflections. The latter stem not only from bubbles, but also from the hull. Then, an area in the PIV frame that does not contain any bubbles is extracted (see Fig. 7). It is located near the bottom of the measurement plane, where a homogeneous seeding distribution is obtained. This area is not supposed to contain any bubbles and is far away from the free surface, avoiding the presence of reflections.

The application of a specific threshold value allows the determination of the edges or contours around the high gradient magnitude. Then, we fill the area inside the detected edges or contours. For instance, if only one pixel exhibits a high gradient magnitude, the edge detection method isolates its eight neighbors and then only nine pixels are really extracted and colored in white. For such a filling process, the *regionprops* command of Matlab is used to obtain properties such as the white areas from the binarized image. It allows the determination of the shape of the connected white pixels. Once this operation is done, we count the number of pixels contained in each white area. Fig. 7 (right) displayed the number of white area as a function of their surfaces (number of white pixels). A value of 9 is mainly recovered, meaning that the main detected area corresponds to exactly one pixel area. This value of 9 should not be related to the possible peak-locking bias error associated with PIV measurements. Moreover, as stated above, this value is directly linked to the specified threshold in Prewitt's method. For other threshold values, the maximum value could be reached for area of one pixel, even in the
<figure><img src="image_8.png" /><figcaption>**Fig. 7.** Step 2-1. Top: Database '1b'. Bottom: Database '1a'. Left: Illustration of the extracted green delimited area. Right: Surface area (expressed in pixels) of the white area detected in the green delimited area.</figcaption></figure>

presence of bubbles. It remains difficult to explicit the origin of the smaller areas detected. They could be related to either
noise, particles, small bubbles, or reflections.

Step 2-2: Removal of the high intensity large area

The next step consists in applying a filtering procedure to remove the white areas that are superior to 9 pixels. An illustration of the observation of such white areas is displayed in Fig. 6. These large white areas are linked either to a conglomerate of PIV particles (in the center of the graph (d)) or to the edge detection and filling processes (on the left of the graph (d)). Note that the white areas of 9 pixels may not be related to only PIV particles in other regions of the PIV frame. Indeed, in the bubble cloud, some bubbles of similar small areas can be observed. We compute the number of 9-pixel white areas contained in a window of 24 × 24 pixels, provided a convenient measure of such an area distribution. This number computed in such a window located inside a bubble cloud area may be 1.5 times higher than that computed in a similar window, but located at the bottom region of the PIV frame without bubbles [25].

This filtering operation permits to remove the large-scale bubbles and their associated reflections (see Fig. 8). This gives a first approximation of the free surface location, except where some particle reflections (9-pixel white area) still occur. These reflections associated with the small white areas (less than 9 pixels) must be suppressed.

4.3. Step 3: Delimitation of the particle-containing area

Step 3-1: Geometric conversion into the PIV mesh grid

As the final mask will be applied to a PIV vector field available on a regular mesh of 66 × 49 points, we propose to convert
the previous investigations into physical coordinates corresponding to the regular PIV grid mesh (see Table 1). As the PIV
mesh discretization corresponds to a 24 × 24-pixel window (databases 1 and 2), we propose to analyze all consecutive
24 × 24-pixel windows. If one contains some white areas, the PIV mesh grid point is fixed to the white color; otherwise, it
is black (see Fig. 9).
<figure><img src="image_9.png" /><figcaption>**Fig. 8.** Step 2-2. PIV raw image after filtering the large high-intensity area for database ‘1b’ (left) and database ‘1a’ (right).</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Fig. 9.** Step 3-1. Database 1b. Conversion to the PIV mesh grid. Analysis from consecutive windows of 24 × 24 pixels (left) and of 32 × 32 pixels² (right). In each graph, the x- and y-axes correspond to the PIV mesh grid ({i_x, i_z} ∈ {1:n_x, 1:n_z} = [66 × 49]).</figcaption></figure>

To overcome the impact of the window's size on data treatment, and knowing that the PIV interrogation window is set to 32 × 32 pixel² (see Table 1), we performed similar treatments for all consecutive windows of 32 × 32 pixels². The results are displayed in Fig. 9. Even if these figures are slightly different when compared, the next steps of the algorithm provide eventually similar free surface detection. In the following, only the results deduced from the investigation of 24 × 24-pixel windows are presented.

### Step 3-2: Filling the region containing particles

Based on a neighbor-point analysis performed from the top-left corner to the bottom-right one, an interpolation between each white area is made to isolate the whole area related to these conglomerations (see Fig. 10, right). Neighbor-point analysis is conducted for each PIV grid point as follows: if at least three of its eight neighbors are white colored, then this grid point is colored in white. Two areas (white and black) are then obtained, delimiting both regions in the PIV frame (see Fig. 10, left). In Fig. 10 (right-hand side), the associated roughly free interface is represented by a black line.

### 4.4. Step 4: Free surface determination

#### Step 4-1: Rough approximation of the free surface

A first approximation of the interface is done by performing a second-order polynomial interpolation. The resulting interpolated curve is indicated with a red line in Fig. 10 (right). This curve is connected to the boundary domain using two extreme points represented by two green squares in Fig. 10. The locations of these two points indicate the tendency of the free surface. As a function of the PIV image, the connecting line between these points (from left to right) increases or decreases. For instance, Fig. 10 (right) illustrates an example of the decreasing line (green dotted line).

#### Step 4-2: Removal of particle reflections

When the connecting line decreases (like in Fig. 10, right), we sweep the PIV mesh grid along each horizontal line from the top right to the bottom left. For each PIV grid point, the analysis of its eight neighbors is carried out. If at least 4 of them
<figure><img src="image_11.png" /><figcaption>**Fig. 10.** Step 3-2 (left) and Step 4-1 (right). Database 1b. Extraction of the mesh grid point containing particles for database '1b' (Left). Rough smooth approximation (red-line) differentiating the two black and white areas. Decreasing line (green dotted line) used for the next step. In each graph x- and y-axes correspond to the PIV mesh grid ($[i_x, i_z] \in [1:n_x, 1:n_z] = [66 \times 49]$).</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Fig. 11.** Steps 4-2 and 4-3. Illustration of the detection of the resulted free surface for database '1b' (top left) and database '1a' (bottom left). The blue line is generated by step 4-2 and the red line by step 4-3. Illustration of the mask generation; the black area is set to zero. The x- and y-axes correspond to the PIV mesh grid ($N_x \times N_y = 66 \times 49$).</figcaption></figure>

are equal to zero, this grid point is set to a zero value. Conversely, if at most 3 of the neighbors are equal to zero, we set this grid point to 1. Then, by sweeping the whole PIV mesh grid, we obtain the results presented in Fig. 11 (blue line).

In case of an increasing tendency as a function of the x-domain, a similar procedure has been applied by sweeping the PIV mesh grid along each horizontal line, but now from the top left to the bottom right.
<figure><img src="image_13.png" /><figcaption>**Fig. 12.** Step 4-3. Raw PIV images superimposed on the free surface (blue dotted line) for database ‘1b’ (left) and database ‘1a’ (right).</figcaption></figure>

### Step 4-3. Free surface determination

The last step consists simply of a fifth-order polynomial interpolation of the result of the previous step. This interpolation (red line) is superimposed to the latter (dotted line) (see Fig. 11). The masking mesh applied to the PIV velocity vector field is then generated (Fig. 11) where a black mesh grid point is set to a zero value and a white grid point is put to a value of 1. An illustration of this application is presented in Fig. 12, in which the final interface is plotted on the PIV raw images.

The whole procedure is applied to the whole PIV database available (see Tables 1 and 2). The algorithm is applied separately to each PIV image independently of the mask’s information related to other PIV images. Recall that for PIV vector field determination, a pair of raw images is required. This algorithm has been tested successively for the two images of the latter and a similar water-air interface has been retrieved. For each considered database, the free interface detection algorithm is only applied to the first image of the pair, and the resulted interface is imposed to the computed PIV vector field. An illustration of the application of the free surface detection algorithm is presented in Fig. 13 during a phase where a bubble cloud is generated (database ‘2c’). Note that the free surface is determined in the PIV mesh grid with a pixel discretization of 24, in each direction (see Table 1), leading to a local error of ±24 pixels in each direction. Thus, when plotting this surface on the PIV raw images discretized in pixels, possible but very small bias can be locally observed.

### 4.5. Free surface properties

When looking at PIV database ‘1’ (independently of the flow conditions), we observe that the wave surface is sometimes above the region of interest. In this sense, it is not possible to properly analyze the water-air interface signal, as it is not continuous in time. The results presented in the following concern exclusively those deduced from database ‘2’, for which the water-air interface can always be detected.

First, Fig. 14 represents the time evolution of the fluctuating surface signal extracted at a fixed upstream location ($x = 76$ mm). It is observed that the amplitude of the free surface is of ±50 mm, which is in a great agreement with previous free surface measurements [2] performed in similar inlet regular wave conditions with an imposed wave amplitude of 33 mm.

Furthermore, to ensure that the periodicity of the incoming wave is recovered, a Fourier analysis of the free surface signal is performed. Fig. 15 displays the spectral representation of the free surface signals obtained in each database for which the imposed frequency ($f_0 = 0.5$ Hz) at the inlet is always recovered. While for databases ‘2a’ and ‘2b’ (without hull motion), some low-frequency peaks are present, we observe that for the other databases (‘2c’ and ‘2d’, with forced hull oscillations), a peak appears at $f = 2f_0 = 1.7$ Hz. This can be directly related to the interaction of hull motions, also generated at the frequency $f_0$, and of the incoming wave.

Moreover, to examine the effect of hull oscillations on the bow wave, the Root Mean Square of the free surface signals is computed. It is then observed that the level of the RMS remains globally the same when far away from the hull (for databases ‘2a’ and ‘2c’ in the symmetry plane). However, the levels of the RMS when close to the hull increase notably by about 9% for database ‘2a’ (without hull motions) and more than 15% for database ‘2c’. This result shows the interaction between the hull motions and the free surface. This leads to increase the free surface variabilities and then the frequency of occurrence of bubble clouds.

## 5. Preliminary analysis of PIV vector fields

In the following, a preliminary analysis of the hydrodynamics is performed, showing the application of the new dynamic masking algorithm detailed above. This mask is applied directly to the PIV vector field. Indeed, the post-processing of the PIV image does not change as a function of the free surface location and of the hull interface. That is the PIV image interrogation window is not locally modified near the interface boundaries. Thus, each instantaneous PIV velocity vector located over the water-air interface is set to zero.
<figure><img src="image_14.png" /><figcaption>**Fig. 13.** Step 4-3. Database '1c'. Illustration of the free surface detection applied to consecutive raw PIV images when the bubble sweep-down phenomenon occurs.</figcaption></figure>

To observe the effect of an hull motion and/or incoming wave on the hydrodynamics, we focus on the representation of
the instantaneous vertical velocity components, w, which is the most sensitive to the hull and wave motions. An illustration
is given in Fig. 16 for database '2c'. This figure displays two instantaneous w components computed from the two raw
PIV images labelled (i) and (f) in Fig. 13. They are represented before and after the application of the mask. Moreover,
the corresponding instantaneous streamlines are also given in this figure. These two figures are retained because they
are representative of the mask application in the presence or the absence of bubbles in the raw PIV images (i) and (f)
respectively. We then observe the downward flow motion (graphs (f)) related to the downward motion of the hull and also
to the wave beginning to get off. Conversely, the vertical velocity component associated with the other instant (graphs (i))
<figure><img src="image_15.png" /><figcaption>**Fig. 14.** Time evolution of the fluctuating free surface signal extracted at a fixed upstream point (x = 76 mm). Top to Bottom: Databases '2a', '2b', '2c', and '2d'.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Fig. 15.** Fourier spectra (expressed in m/Hz) of the free surface signals determined in each database. In each graph, the y-axis corresponds to the frequency axis (expressed in Hz).</figcaption></figure>

exhibit an upward motion of the flow just behind the hull. Note that these instantaneous pictures also emphasize (especially the (f) image containing a lot of bubbles) that it would be necessary to take into account the presence of bubbles in the PIV algorithm. Indeed, the presence of bubbles may modify locally the instantaneous flow field. Such investigations are ongoing.

Furthermore, Fig. 17 presents the time evolution (over more than ten periods of the hull oscillations) of the instantaneous vertical velocity components (w) extracted at a similar point located just above the bulb of the hull, in databases '2a' and '2c'. Note that when dealing with database '2a', this point remains at a fixed position. As for database '2c', in order to consider a similar location, a moving point is extracted at each instant. An illustration of the location is provided in Fig. 16 with a black circle. The vertical velocity component signal is superimposed onto the fluctuating wavy surface deduced from the masking operation and extracted at the similar streamwise position. It is first observed that the amplitude of the wavy surface for database '2c' is globally higher than the one of database '2a'. This is due to the hull motion, which leads to increase the wavy amplitude. Besides, the representations of the instantaneous w component underline the rule of the hull motion onto the hydrodynamics. Indeed, the present analyses seem to demonstrate that the hull motion increases the magnitude of the velocity w. More analyses have to be performed to confirm these preliminary analyses. For instance, the phase average procedure based on a Proper Orthogonal Decomposition and previously implemented [5,26] would allow a better investigation of the influence of the ship's motions on the intensity of the vertical velocity component, which is directly related to bow wave generation.

## 6. Conclusion

The flow at the bow of an oceanographic research vessel is experimentally investigated from PIV measurements with 1/30 ship model. The presence of an oscillating unsteady water-air interface, the hull, some bubbles and many reflections
<figure><img src="image_17.png" /><figcaption>**Fig. 16.** Database 2c. Instantaneous velocity components w without the application of the dynamic mask (top) and after the application of the mask (center). Bottom: Associated streamlines. The labels (f) and (i) indicate the two instants related to PIV images (f) and (i) provided in Fig. 13.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>**Fig. 17.** Time evolution of the instantaneous vertical w velocity component for database '2a' (top) and database '2c' (bottom), extracted at a similar point located above the bulb of the hull. The locations of this point are indicated with a black circle in Fig. 16. The free wavy surface (deduced from the masking operation) is also superimposed to the velocity signal.</figcaption></figure>
leads to generate numerous laser reflections and therefore spurious PIV velocity vectors. In this note, a post-processing method applied to PIV raw images is implemented to detect instantaneously the hull edge contours and the unsteady interface oscillations. A masking procedure is then developed to generate a geometric mask for every pair of PIV images (that is every vector field), and thus eliminate all the vectors that are outside the flow area.

As the presence of the solid hull surface induces laser reflections, the hull edge contours are simply detected in the first PIV frame of each measurement campaign. Then knowing the oscillatory sinusoidal motion of the hull during the measurements, the edge contours are dynamically estimated for all other PIV frames. As for the free surface's determination, a specific algorithm is proposed. It is based on the following process: 1) an edge detection of the high gradient magnitude in the PIV frame in grayscale intensity (this stage is based on the Prewitt's method followed by the morphologic closing method), 2) the detection of particles performed by filtering the high-intensity large area related to bubbles and/or hull reflections, 3) an estimation of a rough region of the PIV image allowing us to recover the particles and their reflections, 4) a rough approximation of the free surface when removing the particle reflections in the PIV image. The free surface is finally obtained through a fifth-order polynomial interpolation. The resulted free surface is validated successfully using Fourier analysis by recovering the initial wave frequency content and visualizing selected PIV images containing spurious high-intensity areas.

Finally, each mask is automatically applied to the corresponding PIV vector field. A preliminary analysis of the resulted PIV vector fields is then conducted, showing the effect of ship motions on the intensity of the vertical velocity component, which is directly related to bow wave generation. Moreover, to improve the present PIV measurements, especially those in the presence of bubble clouds, it would be necessary to estimate the effect of these bubbles on the PIV algorithm and then on the resulted PIV vector fields. Such analysis is ongoing.

## Acknowledgements

D. Dussol would like to thank IFREMER for the financial support of his internship. The authors would like to thank the DGA (French Government Agency for Defence) and the Nord-Pas-de-Calais Regional Council for the financial support of the Ph.D. theses of both S. Delacroix and of B. Mallat. The authors are also grateful to T. Bacchetti, J.-V. Faq and B. Gaurier for their involvement in the PIV measurements.

## References

[1] G. Deane, D. Stokes, Scale dependence of bubble creation mechanisms in breaking waves, Nature 418 (2002) 839–844.

[2] S. Delacroix, Caractérisation de la génération et de la propagation de bulles autour de la carène des navires scientifiques, Ph.D. thesis, Université de Bretagne occidentale, 2014.

[3] S. Delacroix, G. Germain, J. Billard, B. Gaurier, J. Faq, T. Bacchetti, Experimental study of bubbles sweep-down occurrence on oceanographic research vessels, in: OMAE – ASME, 33rd International Conference on Ocean, Offshore and Arctic Engineering, San Francisco, CA, USA, 2014.

[4] S. Delacroix, G. Germain, L. Berger, J. Billard, Bubbles sweep-down occurrence characterization on research vessels, Ocean Eng. (2016) 34–42.

[5] P. Druault, E. Bouhoubeiny, G. Germain, POD investigation of the unsteady turbulent boundary layer developing over porous moving flexible fishing net structure, Exp. Fluids 53 (1) (2012) 277–292.

[6] W. Melville, F. Veron, C. White, The velocity field under breaking waves: coherent structures and turbulence, J. Fluid Mech. 454 (2002) 203–233.

[7] J. Canny, A computational approach to edge detection, IEEE Trans. Pattern Anal. Mach. Intell. PAMI-8 (1986) 679–698.

[8] Z. Huang, S. Hsiao, H. Hwang, Observation of coherent turbulent structure under breaking waves, Int. J. Offshore Polar Eng. 19 (1) (2009) 15–22.

[9] M. André, P. Bardet, Velocity field, surface profile and curvature resolution of steep and short free-surface waves, Exp. Fluids 55 (1709) (2014) 1–19.

[10] R. Theunissen, F. Scarano, M. Riethmuller, On improvement of PIV images interrogation near stationary interfaces, Exp. Fluids 45 (2008) 557–572.

[11] Y. Jeon, H. Sung, Piv measurement of flow around an arbitrarily moving body, Exp. Fluids 50 (2011) 787–798.

[12] J. Park, S. Im, H. Sung, J. Park, Piv measurements of flow around an arbitrarily moving free surface, Exp. Fluids 56 (2015) 1–16.

[13] A. Sanchis, A. Jensen, Dynamic masking of PIV images using the radon transform in free surface flows, Exp. Fluids 51 (2011) 871–880.

[14] A. Ayati, J. Kolass, A. Jensen, G. Johnson, A PIV investigation of stratified gas–liquid flow in a horizontal pipe, Int. J. Multiph. Flow (2014) 129–143.

[15] G. Gomit, L. Chatellier, D. Calluaud, L. David, Free surface measurement by stereo-refraction, Exp. Fluids 54 (6) (2013) 1540–1550.

[16] S. Delacroix, G. Germain, J. Billard, B. Gaurier, J. Faq, Étude expérimentale du phénomène d'entraînement de bulles à l'étrave des navires scientifiques, in: 14es Journées de l'hydrodynamique, Val-de-Reuil, France, 18–20 November, 2014.

[17] Y. Cheng, S. Pothos, F. Diez, Phase discrimination method for simultaneous two-phase separation in time-resolved stereo PIV measurements, Exp. Fluids 49 (2010) 1375–1391.

[18] E. Bouhoubeiny, G. Germain, P. Druault, Time-resolved PIV investigations of the flow field around rigid cod-end net structure, Fish. Res. 108 (2–3) (2011) 344–355.

[19] B. Gaurier, P. Davies, A. Deuff, G. Germain, Flume tank characterization of marine current turbine blade behaviour under current and wave loading, Renew. Energy 59 (2013) 1–12.

[20] S. Delacroix, G. Germain, B. Gaurier, J. Billard, Experimental study of bubble sweep-down in wave and current circulating tank. Part I: Experimental set-up and observed phenomena, Ocean Eng. (in revision).

[21] J.D. Kat, J. Paulling, Prediction of extreme motions and capsizing of ships and offshore vehicles, in: OMAE Conference, 20th Int. Conf. Offshore Mech. and Arctic Eng. No. OMAE2001/OFT-1280, 2001.

[22] J. Prewitt, Object Enhancement and Extraction, Academic Press, New York, 1970.

[23] F. Guo, Y. Yang, B. Chen, L. Guo, A novel multi-scale edge detection technique based on wavelet analysis with application in multiphase flows, Powder Technol. 202 (2010) 171–177.

[24] N. Deen, P. Willems, M.V.S. Annaland, J. Kuipers, R. Lammertink, A. Kemperman, M. Wessling, W.V.D. Meer, On image pre-processing for PIV of single- and two-phase flows over reflecting objects, Exp. Fluids 49 (2010) 525–530.

[25] D. Dussol, Étude de l'écoulement autour de l'étrave des navires océanographiques (in French), Master's thesis, Université Pierre-et-Marie-Curie, Paris-6, 2015.

[26] E. Bouhoubeiny, P. Druault, G. Germain, Phase-averaged mean properties of turbulent flow developing around a fluttering sheet of net, Ocean Eng. 82 (2014) 160–168.