Visualizing large-scale flow using synthetic aperture PIV

We discuss the application of synthetic aperture particle image velocimetry for measuring the flow around human swimmers using small bubbles as tracer. We quantify the two-dimensional projection of the velocity field in planes perpendicular to the viewing direction of an array of six cameras. With help of simulations, modelled after the experiment, we address questions about depth selectivity and occlusion in dense bubble fields. Using vortex rings in the swimming pool, we provide a proof of principle of the method. It is further illustrated by the vorticity field produced by a human swimmer.


Introduction
In its simplest guise, synthetic aperture particle image velocimetry (SA-PIV) is a computationally cheap and efficient way to focus on planes throughout a three-dimensional (3D) volume and obtain 2D velocity information in these planes. In more elaborate applications (Mendelson and Techet 2015), a complete three-dimensional velocity field was provided by SA-PIV. In this paper, we will explore this technique in its simplest form to quantify large-scale flow in a swimming pool. The purpose of our research is to visualize and quantify flow around a swimmer. We describe a system that has been integrated in a swimming pool, uses small air bubbles as tracer and synthetic aperture particle image velocimetry to obtain the velocity field in planes perpendicular to the viewing direction. The distance of these planes to the camera system (their depth) can be selected by refocusing the images (movies) after they have been taken. This may help in identifying flow structures produced by a swimmer. For example, the strength of vortex rings is most obvious in planes that dissect these rings perpendicularly. Depth selectivity was achieved with six cameras that simultaneously record images of a 1 m-thick curtain of small bubbles that emanate from five line sources in the (false) bottom of the swimming pool. The resulting depth resolution is approximately 8 cm. The bubbles are illuminated by ambient lighting. The standard way to visualize flow involves the use of a light sheet from a pulsed laser and solid particles as flow tracer. Clearly, this is not an option in the training swimming pool, where our system was implemented. 1 Although it is not our prime goal, our method could be used to obtain the 3D velocity field. Other, superior, volumetric methods exist to obtain this goal. In defocusing digital particle image velocimetry, a mask with multiple apertures (usually three) before the lens is used to obtain multiple images from the tracer particles. With pattern recognition, the points belonging to a single tracer particle can be found. The spread between those points is related to the third spatial component. Because the particle field must be reconstructed from multiply-imaged tracer particles, the main limitation of this technique is the low particle image density [typically ≈ 0.034 particles per pixel (ppp)] (Pereira and Gharib 2002).
In tomographic PIV, the full three-dimensional volumetric velocity field is measured via the reconstruction of the three-dimensional particle distribution and calculation of the volumetric cross correlations (Elsinga et al. 2006). Particle locations are reconstructed by applying optical tomography on images from cameras with different points of view. The seeding density that could be reached in this technique is larger ( ≈ 0.024−0.08 ppp ). The viewable depth (out-of-plane) is typically three to five times smaller than the in-plane dimensions. Very high particle image densities ( 0.125 ppp ) were accommodated in a 3D particle tracking strategy for tomographic velocity field measurements (Schanz et al. 2016). At this moment, the computational effort for these 3D techniques prevents us from providing an almost instantaneous view of the velocity field in planes that can be selected a posteriori.
Synthetic aperture refocusing is a 3D vision method that does not rely on the identification of individual particles and reconstruction of their spatial coordinates. The idea is to achieve depth selectivity from the parallax of multiple cameras by judiciously shifting camera images and summing them. Bubbles that reside at the chosen depth overlap in the summed (refocused) image with a resulting image intensity multiplied by the number of cameras, and other bubbles only contribute to a low-intensity background. By imposing an intensity threshold on the summed image, the depth plane can be selected in the post processing stage (Belden et al. 2010(Belden et al. , 2012Langley et al. 2014;Mendelson and Techet 2015). Since depth selection is fast and simple, this is our method of choice and can be used in providing almost instantaneous feedback to swimmers. The particle image density reported is approximately 0.015−0.125 ppp (Belden et al. 2010). Because this method is less sensitive to particle occlusions, it also enables a larger depth of view compared to tomographic PIV.
In Sect. 1.1, we briefly discuss the SA-PIV technique. We describe the experimental setup in the swimming pool in Sect. 2. In Sect. 3.1, we discuss the application of SA-PIV to a synthetic bubble distribution that is modelled after the actual situation in the swimming pool. Issues considered are the depth selectivity of the method, and the influence of ghost bubbles. Finally, in Sect. 4, we present the result of experiments on the bubble sheet and on the motion of vortex rings in the swimming pool. Figure 1 illustrates the principle of SA-PIV; it is tailored to our situation. Bubble images are recorded by six cameras simultaneously. Due to parallax, the bubble images arrive at different locations at the camera image planes. By applying shifts which depend on depth, the camera images of a bubble for a particular depth Z can be made to coincide. These image shifts follow from a calibration. The depth of field of each camera should cover the entire measurement volume, while their fields of view (FOV) should overlap as much as possible to maximize the FOV of the resulting refocused image.

Synthetic aperture refocusing
By summing, the images of the different cameras are combined into one refocused image with a narrow depth of field on a chosen focal plane. The intensity contrast between in-focus and out-focus bubbles is determined by the number of cameras; the image contrast becomes better with more cameras. Thus, in our application, the camera images are combined as (Belden et al. 2010) (1) the refocused image on the kth plane, N cam the number of cameras, and I Z k ,i the transformed image of the ith camera for the kth plane. From the pinhole model in Fig. 2, we can derive a relation between the depth selectivity and the camera parallax. The position of the image of a bubble with world coordinates (X A , Z A ) on the image plane of camera i is with s i the distance between pinhole and image plane of camera i, s 0 the Z coordinate of the pinhole and X C i the X coordinate of the pinhole of camera i. To correct for parallax of two cameras, the image of camera j must be shifted with respect to that of camera i by Equation 2 embodies the principle of refocusing. In practice, however, the shifts are inferred from a calibration procedure which allows for higher order corrections, such as lens errors, refraction at the interface between camera housing and water, and viewing angle. Now, assume that we move a bubble from depth Z A to Z A + Z , but with its X coordinate unchanged. Refocusing this bubble requires a change in image shift with D ij = X C i − X C j the camera parallax. When d ij decreases and becomes equal to the diameter of a bubble image, the depth difference Z can no longer be resolved. If b is the actual bubble diameter, it is approximately b s i ∕(s 0 + Z A ) on a camera image. The equation , with d ij given by Eq. 3 determines the depth resolution Z . To first order in b ∕D ij , its solution is This simple argument can be refined by considering the intensity distribution of a bubble, and when applied to PIV, the correlation function of bubbles at different depths. These refinements will be addressed using a faithful numerical simulation of our setup in Sect. 3. (3)

Experiment
In this section, we discuss our choices for the camera arrangement, its calibration and the properties of the used flow tracers. A detailed account can be found in van Houwelingen (2018).

Cameras
We used six monochrome cameras (Sony: XCG-CG240) with a full frame resolution of 1920 × 1200 pixels at a frame rate of 41 fps. To obtain higher frame rates of 50 fps, the height of the active region on the CMOS sensor is cropped to 1900 × 956 pixels. Larger frame rates would have demanded more advanced illumination of the scene and larger bandwidth of the camera to computer connection. The cameras are arranged as a hexagon with sides of 0.3 m.
All cameras have a 16 mm lens (KOWA: LM16HC F 1.4 f 16 mm). This camera-lens combination results in a depth of field of approximately 2 m , causing the complete depth of the measurement volume to be in focus, as required for SA-PIV. The midplane of the measurement volume is at Z = 5.95 m ; there, the field of view is approximately 3.1 × 1.5 m 2 , which is sufficient to capture the flow around a complete swimmer. In this setup, one pixel corresponds to approximately 1.6 mm in world units.
In a typical PIV experiment, interrogation windows of 32 × 32 pixels are used with a 50 % overlap region. This would result in a spatial resolution in the viewing plane of approximately 25 mm . Considering the PIV design rule for in-plane motion (in-plane shifts smaller than 32∕4 = 8 pixels, velocities of 0.65 m/s can be resolved with ease (Adrian and Westerweel 2011).
Close-up images in a small tank showed that the typical bubble diameter is 4 mm, so that Eq. 4 predicts a depth resolution of Z ≈ 8 cm . The chosen camera spacing d cam is a compromise between depth resolution, which increases with increasing d cam , and field of view, which decreases with increasing d cam . Summarizing, the spatial velocimetry resolution X, Y, Z is 2.5, 2.5, and 8 cm , respectively.

Calibration
Calibration of the camera images is necessary for accurate shifting (and warping) before summing them at the chosen (refocused) depth Z. The calibration frame supports a banner of 3.0 × 1.5 m 2 . The black banner was filled with white dots of 10 mm in diameter every 80 mm. Three dots in the centre of the grid were made identifiable to recognize origin and orientation of the calibration field. The frame was traversed in the water, through the measurement volume above the bubble system, at specified distances Z k from the camera array ranging from 5450 to 6450 mm with increments of 125 mm.
This sets the planes of interest for the refocusing method. At each distance, a record is made with the six cameras. Because every camera views the same measurement volume from a slightly different location, the images are shifted (and warped) compared to each other. At each position, a function is determined to map known world coordinates to camera pixel values. The calibrations also take care of lens aberrations originating from the camera lens and refraction at the window of the underwater casing. The positioning of the calibration frame-by handsuffers from inaccuracies. These inaccuracies are not relevant when refocusing at the calibrated distances Z k . As explained below, alignment errors may be corrected through a pinhole camera model.
The world coordinates (X, Y) of the dots on the grid are known. The pixel coordinates (x, y) of the dots in the images are automatically detected using a high-pass filter on the image, followed by a least-squares fit to find the centroid of the intensity. The dots in the image are linked to those on the calibration grid. From this association, a mapping function For N cam = 6 cameras and M = 8 calibration planes, this yields N cam × M = 48 mapping functions. Each of these functions is in the form of a second-order polynomial relating pixel coordinates at camera i to planar world coordinates at depth k: with coefficients a i,k lm , b i,k l,m , l = 0, 1, 2;m = 0, 1, 2 determined in a least-squares procedure.
The accuracy of the mapping was checked by comparison of the known world coordinates of the calibration grid markers, and the ones inferred from the transformed camera images. The average error was ≈ 0.45 pixels. Belden et al. (2010) argue that the error in the mapping functions should be less than 0.45 pixels for a full 3D reconstruction of particle locations. Our implementation of SA-PIV is more modest: we are interested in 2D velocity fields in planes, whose distance to the camera array in the viewing direction can be selected.
An interesting corollary of the quadratic mapping functions is that camera locations, camera viewing angles and alignment of the calibration grid can be reconstructed using a pinhole model. The coefficients of the zeroth-and first-order polynomial terms determine camera placement X c , Y c and magnification, from which the Z coordinate of the camera location can be found. When the grid is not placed precisely at the intended Z-positions Z k , this affects the calibrated magnification of all cameras simultaneously.
In turn, this leads to a correction on the used Z k values. In addition, the camera viewing angles and grid orientation can be determined from the quadratic polynomial coefficients. A change of these angles with position Z k provides direct information on misalignment of the calibration grid.
The calibration functions associated with the positions Z k , k = 1, … M allow refocusing on these planes. Inbetween depths can be refocused by interpolating between the corresponding polynomial coefficients. Concluding, with our mapping functions the planes of interest can be reconstructed with an accuracy up to 1 pixel, except for the corner regions of the recordings, where the influence of lens deformations is largest. The disadvantage of our calibration procedure is that there is no obvious and simple relation between the polynomial coefficients and viewing geometry. More advanced methods use projective geometry, which obfuscates the necessity of precise knowledge of the calibration grid location (Tsai 1987;Zhang 2000). Specifically, a calibration method for cameras imaging through refractive interfaces was described by Belden (2013); in addition, in this method, precise knowledge of the location of the interface and calibration world points is not needed. Higher accuracies may be possible by the use of volumetric self-calibration techniques (Wieneke 2008;Svoboda et al. 2005;Faugeras et al. 1992).

Image preprocessing
Since SA-PIV relies on intensity information in images, image preprocessing is necessary before application of Eq. 1 (Belden et al. 2010). Briefly, the following steps are taken: (1) remove background by subtracting an average background recording; (2) obtain a flat field image from an average bubble recording and divide all images by this flat field; (3) transform the image to world coordinates and crop to the boundaries of the refocused frame; (4) equalize histograms of all six frames; (5) enhance the contrast of bubbles by stretching intensities; (6) obtain the refocused image by summing the images of the different cameras; and (7) threshold the refocused image.

Bubbles
To observe flow motion, bubbles are used as tracer particles. The bubble system is integrated in the bottom of the swimming pool (see Fig. 3). It covers a distance from 5 to 10 m from the starting platform, so that the FOV of the cameras is in the center of the generated bubble curtains, which is at Z = 5.95 m from the side wall. The bubble system consists of five parallel PVC tubes with length of L t = 4.5 m , placed 0.25 m apart. Each tube has a series of small holes with a separation distance of 0.02 m along its length. The holes are laser-drilled and conically shaped with smallest diameter of 0.2 mm on the inside of the PVC tube.
The air flow is supplied by a compressor (Leonardo 201). Four filters (Metal Work Pneumatic: 20 μm, 5 μm, 0.01 μm , and an active coal filter) are used for bubbles with clean breathable air. The air flow through the bubble system is Fig. 3 Schematic plan view of the bubble system placed in the swimming pool (van Houwelingen 2018). The camera system is implemented in a special double wall of the swimming pool. To increase the image contrast, a black canvas is placed on the wall opposite the camera system. On the right side, a schematic frontal view of the camera array is shown. The dashed-dotted line indicates the FOV of the cameras. The camera spacing D i,j ranges from 300 to 600 mm, but for an estimate of the depth resolution, we take d cam = 300 mm Five homogeneous bubble curtains are produced with an air flow rate of Q air ≈ 0.01 m 3 ∕min per tube. The air bubbles have a diameter of approximately 4 mm. Taking a rising velocity of approximately 0.3 m/s (see Sect. 4.2), we compute an average vertical spacing of 1.4 cm. From these numbers, we estimate a total of 1.7 × 10 5 bubbles in the ( 4.5 × 1.9 × 1.0 m 3 ) volume above the bubble system and 9 × 10 4 bubbles in a camera field of view. The volume fraction of the air in the bubble field is 6.2 × 10 −4 ; the decrease of the water mass density is not perceptible for a swimmer.
Although the experimental bubble field appears dense, the number of bubbles per pixel is 5 × 10 −2 . However, what matters for our application is the number of bubbles in a refocused slice of the bubble volume, which is approximately one order of magnitude smaller. With these numbers, the number of bubbles in a refocused PIV interrogation window of 32 × 32 pixels is a mere 5.

Simulation
It is useful and instructive to study the feasibility of synthetic aperture PIV in a simulation. The simulation is designed after the experimental setup and provides an idealized world in which all bubbles have the same intensity profile. This makes refocusing, which is based on intensity discrimination, tractable. We study the ability of refocusing to reproduce bubble density variations in a direction perpendicular to the image plane. We will do the same in the experiment in Sect. 4. Next, we will ask ourselves whether SA-PIV can find isolated velocity structures hidden in a dense field of bubbles.
Refocusing is thwarted by the emergence of ghost bubbles and occlusions. Ghost bubbles are bubbles that live at other depths than the refocused depth, but accidentally fall on the proper location in several camera images. Ghost bubbles become more numerous when the bubble density increases. In Sect. 3.1, we will present a statistical model for the chances to find a ghost bubble, and compare its prediction to the simulation.
The simulation uses the exact same geometry as the experiment, with cameras at exact the same positions. However, the imaging of a camera is done using a simple pinhole model. We randomly sprinkle N bubbles with radius r in the three-dimensional measurement volume. The idealized intensity profile of each bubble i is where the bubble center is assumed in the origin. The maximum intensity of these bubbles is 1, and the integrated intensity equals the bubble area A b . The images of bubbles can overlap: bubbles closer to the camera may obscure bubbles that are further away. This was modelled by sorting the bubble coordinates with respect to their distance Z to the camera. Camera images are formed by mapping the bubble world coordinates (X, Y) onto the camera pixel map, and painting bubble images in sorted order, so that bubbles closer to the camera obscure those further away. (The bubble images are not diffraction limited.) Unlike in the experiment, no image processing was done before refocusing. A camera image is shown in Fig. 4.
Much as in the experiment, we use the synthetic aperture technique to resolve the variation of bubble density as a function of depth Z. We count the number of bubbles (5) I i (x, y) = 1 − (x 2 + y 2 )∕r 2 for x 2 + y 2 < r 2 , and 0 otherwise , N ′ in the refocused image by integrating its intensity and dividing by the intensity A b of a single bubble. Due to spatial discretization, A b is, on average, slightly smaller than the nominal A b = 1 2 r 2 . In the experiment, we will measure the density profile of single-and multiple bubble lines (Fig. 5).
To first order in the bubble density, the number of bubbles seen by a single camera is the total number of bubbles, L x , L y the width and height of a camera FOV, where the term quadratic in N reflects occlusion. As shown in Fig. 8, this simple formula agrees well with the number of bubbles observed in the simulation. For a number of bubbles N ≤ 10 5 , less than 20% of the bubbles are occluded, but this number rapidly increases with increasing N.

Ghost bubbles
A key problem of synthetic aperture refocusing is the emergence of ghost bubbles. The idea of ghost bubbles is illustrated in Fig. 6a. A ghost bubble is formed by bubbles outside the refocused plane that accidentally appear at the proper location in several camera images simultaneously. If this is so for all N cam camera images, it is not possible to eliminate the ghost bubble. If this is so for some (m) camera images, the ghost bubble may be removed by imposing an intensity threshold. Therefore, being able to impose an intensity threshold on the refocused image is key for eliminating ghost bubbles with m < N cam . The illumination in our simulation is ideal: all bubbles have unity maximum intensity. Therefore, an intensity threshold t SA admits ghost particles that appear in m camera images. Thus, we have introduced the concept of m−ghosts: a ghost bubble involving m camera images only. Clearly, to kill all mortal m < N cam ghosts, the threshold has to be chosen t SA = (N cam − 1)∕N cam with respect to the maximum (unity) intensity in our simulation. In this context "mortal" refers to ghost bubbles that can be eliminated by imposing the appropriate threshold. Whether this is feasible or not in practical situations depends on the quality of the illumination.

Statistics of ghost bubbles
We will now discuss a simple statistical model of ghost bubbles, and compare the outcome to the result of the simulation. The model is used to compute the chances of finding a ghost bubble. Each camera image of a bubble projects a tube into the bubble volume with cross-sectional area A b (nominally A b = 1 2 r 2 ). In-focus bubbles are located in the intersection of all N cam tubes and the focal plane.
The chances of finding ghost bubbles can be computed from the following argument. Since the parallax of the cameras is small, the volume of a camera tube is A b L z , with L z the depth of the measurement volume. Then, N cam A b L z is the bubble tube volume for all cameras together. The chance for a bubble to end up in this joint volume is p = N cam A b ∕(L x L y ) . The chance P that out of N bubbles, k bubbles end up in this joint tube volume is binomial: For m camera tubes to be populated simultaneously, we have to multiply P with the probability that each of the m tubes contains at least one bubble. These probabilities are again binomial, with (1 − 1∕m) k the probability that none of the k bubbles falls in a tube, and thus 1 − (1 − 1∕m) k the probability that at least one of the k bubbles falls in a given tube, and (1 − (1 − 1∕m) k ) m the probability that this occurs for m tubes simultaneously. Summarizing, the probability for a ghost bubble seen by m cameras in a collection of N bubbles is As N is large, and p is small, the sum in this equation extends over a few k only (Fig. 7). Because in our simulation the positions of the bubbles are known, we can quantify exactly the depth resolution and the chances of ghost bubbles. To this aim, we compute the correlation between the refocused image at depth Z 0 and the exact intensity field I i of bubbles i in bins [ − ∕2, + ∕2]: We refocus at depth Z 0 = 0 and normalize the correlation function by the integrated squared intensity ∫ ∫ I 2 (x, y) dxdy of a single bubble image times the number of bubbles in the bin with width at depth Z. Alternatively, we can normalize by setting C(Z = 0) = 1 . If no ghost bubbles are present, C(Z) is a peak with maximum 1 sitting on a background 0. The width of the peak is the refocusing depth resolution. Ghost bubbles cause a background with height P (m) ghost . The result for several values of N and m is shown in Fig. 6.
Equation 8 predicts that the chances that a bubble image is a ghost rapidly increase with increasing number of bubbles N as P (m) ghost ∝ N m . This behavior can be observed in Fig. 6, where the background dependence on N approximately follows straight lines with slope m in the log-log plot. The agreement between Eq. 8 and simulation is not perfect, because the bubble images are not uniformly illuminated disks. Due to spatial discretization, their maximum intensity is smaller than 1 most of the time, so that m deviates from m = N cam t SA + 1 . In addition, at large numbers of bubbles N occlusion occurs, which is not contained in the statistical model described by Eq. 8.

Resolving density variation
In the experiment, we use SA-PIV to resolve the variation of the bubble density while scanning in Z through the rising plume of a single bubble line. Similarly, we study the bubble density in the case that all five bubble lines are active. In the simulation, we model this by concentrating the bubbles in a slab with width = 0.05 m at Z = 0 , or by concentrating them in five slabs with the same width, at the locations of the bubble lines. We then cut trough the Z-dependent density field by refocusing at varying depth Z and count the number of bubbles N ′ in the refocused image by integrating its intensity and dividing by the intensity A b of a single bubble image. The resolution of density modulations is shown in Fig. 8a, b for a refocusing intensity threshold t SA = (m − 1)∕6 with m = 2 , which is the smallest value that can detect a density modulation. Apart from the effects of occlusion and the widening of the FOV, for a threshold with m = 1 , all bubbles are always seen in the refocused image, independent of Z. The increase of the width of the refocused slab from 0.05 to 0.12 m is consistent with the estimate Eq. 4 of the depth resolution Z ≈ 0.08 m . The refocused number distribution of an array of slabs is shown in Fig. 8b. Because of the finite depth of focus, the number of bubbles in a single slab is increased through leakage out of its neighboring slabs, which results both in widening and an increase of the height of the peaks.

Particle image velocimetry
In our experiment, we use bubbles to measure the components of the flow velocity in planes perpendicular to the line of view; these planes are depth resolved. In particle image velocimetry, two subsequent refocused images are correlated, and the location of the peak of the correlation function denotes the shift of the particles.
When the velocity field depends on Z, the bubbles that make a ghost bubble each have a displacement that is different from the other bubble in the next snapshot. When they again make a ghost bubble, this ghost is unrelated to the first one, and only contributes to the background of the PIV correlation function. On the other hand, when the velocity field does not depend on Z, ghost bubbles in one snapshot remain so in the next one. Since the velocity field is now two-dimensional, and since the measurement volume is far away from the camera, so that parallax is small, the velocity error is small. Thus, the correlation step in PIV may help suppress false velocity information produced by ghost bubbles (Elsinga et al. 2011;Elsinga and Tokgoz 2014).
Let us now refine this idea. Clearly, the PIV error caused by ghost bubbles depends on the magnitude of the velocity gradients u∕ z and v∕ z , with u, v the velocity components in the plane. When the magnitude of these gradients is large enough, the separation of the bubbles in Fig. 6 may grow larger than their diameter, and they no longer make a ghost bubble. The typical Z-separation l * z of these bubbles, such that they cease to make a ghost bubble depends on the velocity gradient and the time t between two snapshots. When only a gradient u∕ z is present, l * z = r| u∕ z| t. Accordingly, when the gradient stretches across all of L z , the volume of potential ghost bubbles is reduced by a factor l * z ∕L z . Elsinga et al. (2011) provide a statistical analysis of the influence of ghosts on PIV, but it is for large-scale gradients.
The influence of ghost bubbles on velocity structures is illustrated in Fig. 9. By increasing the refocusing intensity threshold t SA = (m − 1)∕N cam , which is equivalent to increasing the number of cameras m that simultaneously view a ghost bubble, the severity of the influence of ghosts on measured quantities can be accessed. In PIV, the particle displacement is = u t , with associated displacement gradient ∕ z = t u∕ z . In the simulations, we can only prescribe a displacement gradient ∕ z . Figure 9 shows that for N = 65536 bubbles, a displacement gradient ∕ z can still be discerned for m = 2 , but not a velocity step in an otherwise stationary velocity field. Only for a larger threshold ( m = 3 ) can both velocity structures be retrieved faithfully. This result shows that the requirement on SA-PIV to detect small-scale velocity structures is more severe than for detecting large-scale gradients. As higher intensity thresholds must be used in the first case, too few refocused bubbles may remain for velocimetry.

Experimental results
Following the simulations with simulated bubble density and velocity fields, we now present the results of experiments. As in case of the simulations, particle image velocimetry on the refocused images was done using 32 × 32 pixels interrogation window sizes with 50% overlap. Outliers were detected and replaced through bilinear interpolation (Westerweel and Scarano 2005).

Density profile of a bubble curtain
At first, a single bubble curtain at Z = 0 is activated. Refocused images are created on planes throughout the whole measurement volume with a separation of 10 mm. The intensity in each refocused image is summed as a function of depth Z. To test the effect of the intensity threshold t SA on the refocused images, this analysis is repeated for different threshold values. Unlike in the simulation, where we expressed t SA in "fraction of cameras viewing a bubble", we now express the threshold value in counts, with 4096 counts saturating the CMOS image sensor.
The results are shown in Fig. 10. Without thresholding, no bubble density variation can be detected. The slow increase with increasing distance to the camera array (decreasing Z) is related to the slightly increased field of view when shifting towards planes at the back of the measurement volume. When the threshold is increased, the peak around Z = 0 becomes more pronounced. With thresholds of about 85 counts and higher, a clear distinction can be made between the peak corresponding to the bubble curtain and the "noise" in front of and behind the bubble curtain, with almost one order of magnitude difference. The width of the peak is related to large-scale and slow collective motion of the bubble curtain, combined with the depth resolution of SA-PIV. The density profile is approximately Gaussian, thus proportional to exp(−Z 2 ∕ Z 2 ) , with Z ≈ 93 mm (full-width at half-maximum 150 mm). From scanning through single bubbles (a single slice is shown in Fig. 10b), it is found that a bubble remains in focus over approximately 60 mm . This corresponds to the estimate of the depth resolution Z in Sect. 2.1.
The integrated intensity of a refocused bubble, averaged over 20 bubbles is approximately 2 × 10 3 counts at t SA = 170 counts (Fig. 10b). Thus, the integrated number of  Fig. 10a is 4.1 × 10 4 . However, the measured density profile of Fig. 10a is close to the refocusing depth resolution Z of SA-PIV. In this way, bubbles are counted several times, and this number cannot be compared to the number of bubbles estimated from the air flow rate ( 1.7 × 10 4 ). Incidentally, at an intensity threshold t SA = 170 counts , Fig. 10a shows that the apparent bubble density is reduced by a factor ≈ 0.15 . This would correspond to m ≈ 6 in Eq. 6. It must be realized, however, that the bubbles in the experiment have a broad distribution of intensities, unlike those in the simulation. When turning on multiple bubble curtains, the planes of interest are expected to experience more noise, caused by the presence of a large amount of bubbles in front of and behind a particular plane. Similar to the case of a single bubble curtain, a record is obtained with the five bubble curtains active. In Fig. 11, the results of the summed intensities are shown for a threshold of t SA = 250 and 92 counts. For both thresholds, the five bubble curtains can be distinguished. As expected the background level for a threshold of t SA = 92 is higher, because more ghost bubbles contribute. The peak heights corresponding to the two bubble curtains in the back are reduced. This may partly be explained by the slightly lower air flow rate through the tubes in the back. V (m/s) Fig. 12 Rising velocity of bubbles in the bubble field. a Shows the experiment with Z dependence measured through refocusing; averages are done over the entire refocused plane. The full line is the mean, the dashed lines indicate the standard deviation. b Shows the mean rising velocity in the central Z = 0 refocused plane. c Shows the simulation with set bubble displacement between two PIV frames of 6 pixels, as indicated by the gray line). We used N = 65536 bubbles and a refocusing intensity threshold m = 2 . In frames a and c the apparent bubble velocity increases when moving away from the camera array. Notice that the slope of the lines in a and c is comparable, which proves that this is a geometric effect 1 Page 12 of 14

Velocimetry
The bubbles have a life of their own: they rise, while a bubble swarm entrains fluid. Their rising velocity causes a bias of flow velocity measurements; it was removed from the measurements that are shown in Figs. 13 and 14. When a bubble is in steady motion, its terminal rising velocity V t can be calculated using the force balance between buoyancy and drag: ( l − g ) g 4 3 r 3 = 1 2 l V 2 t A C D , with l the density of water, g the density of air, g the gravitational acceleration, and A the area of the bubble projected on a plane perpendicular to the direction of motion and C D its drag. With C D ≈ 0.85 for Re ≈ 10 3 , where the bubble diameter is the characteristic length, this yields a rising velocity V t ≈ 0.25 m/s (see Kulkarni and Joshi 2005;Sommerfeld et al. 2003).
In Fig. 12a, the mean rising velocity of bubbles as a function of Z is plotted. The root-mean-square fluctuations are caused by PIV errors due to the limited number of bubbles in a refocused plane, combined with the zig-zag motion of bubbles (van Houwelingen 2018). From this figure, we conclude a mean rising velocity of 0.3 m/s, which is significantly larger than the one predicted. The mean Fig. 13 a Sketch of the 0.7 m-long vortex generator, with radii R o = 0.125 m, R i = 0.095 m . Water is expelled by pushing the piston. The induced circulation is ≈ 1 2 U 2 i T , with T the duration of the expulsion phase. b Visualization of a vortex ring, the contrast is enhanced following the preprocessing steps for refocussing (2.3). The vortex ring is moving from left to right. It is clearly visible, because bubbles are entrained into the vortex core.  Fig. 12b. It was averaged over 2 s . This field forms the background of the vorticity measurements in Figs. 13 and 14.
As one can observe, the apparent rising speed increases when moving towards planes in the back of the measurement volume. This apparent increase is due to out-of-plane bubbles in front of the plane of interest. They appear to move faster when showing up in planes towards the back of the measurement volume. The same phenomenon can be observed in the simulation, as shown in Fig. 12c. Both experiment and simulation show that the effect becomes smaller when increasing the intensity threshold t SA . For the simulation, we imposed a displacement of the entire bubble field between two PIV frames, while the Z dependence was again done through refocusing. In this case, the fluctuations are caused by PIV errors only.
The apparent increase of the rising velocity is the same in experiment and simulation, showing that this is a geometric effect. The simulation shows that the true rising velocity is measured in the plane closest to the camera array.
Another effect of the bubble field is its forcing on the water, resulting in an additional vertical velocity and a large-scale circulation in the pool. An estimate of the induced velocity V w follows from a simple and elegant integral momentum flux argument by Abraham (1974). For a constant vertical induced velocity, it is important that the plume widens. For this, an empirical geometrical factor k was established. The induced water velocity V w is then The effect is larger for smaller terminal rise velocities, as more bubbles provide buoyancy forces. With k = 32 (Abraham 1974), we find V w = 0.08 m∕s . Added to the rising velocity predicted from the bubble size, V t ≈ 0.25 m/s , this gives approximately 0.33 m/s , which can be compared well to the observed rising velocity of 0.3 m/s.

Vortex rings
A vortex ring is a well-defined and confined three-dimensional coherent flow structure and is, therefore, ideally suited to test the flow visualization. In a still medium, a vortex ring is stable and steadily propagates in a direction perpendicular to the plane of the ring. For finite Reynolds numbers, there is some decay related to viscous diffusion (Batchelor 2007).
The vortex rings are created with a vortex ring generator which is operated in the swimming pool. The vortex generator, shown in Fig. 13a, consists of a 0.7 m-long perspex tube with a radius R o = 0.125 m . At the frontal side of the tube, a ring with a inner radius R i = 0.095 m is attached. To help separation, the ring has a sharp edge. When the plug is pushed out with constant velocity U p , mass conservation dictates a flow velocity through the orifice U i = U p (R o ∕R i ) 2 . In a time T, the duration of the expulsion phase, the thin boundary layers at the edge of the inner ring roll-up into a circular vortex, with circulation ≈ 1 2 U 2 i T and radius R ring . This vortex ring propagates with induced velocity U ring = 4 R ring log R ring d ring , with d ring the thickness of the ring.
In Fig. 13b, a typical raw record of a vortex ring is shown. The vortex rings are created at the left side, just outside the field of view and the vortex rings have a self-induced motion towards the right. The circulation in three depth planes is shown in Fig. 13c. It was calculated over circles with diameter 0.13 m , which is approximately the extent of the vortex blobs. The maximum circulation found at Z = 0.10 m is ≈ 0.06 m 2 ∕s . Finally, the associated velocity vector fields are shown in Fig. 13d. Taking a crude estimate for the plug velocity U i ≈ 0.2 m∕s , and a vortex roll-up time of T = 2 s , we expect a circulation ≈ 0.1 m 2 ∕s . However, the vortex ring is far from ideal, as it is distorted by the rising flow field. From the circulation plots at different planes, it can be concluded that depth selectivity is important to reveal the dipolar structure of the vortex ring cross section.
We finally shown in Fig. 14 the circulation field left behind a human swimmer performing butterfly kicks (van Houwelingen 2018). The circulation is computed over circular loops with diameter 0.13 m . Two clear dipolar structures can be seen, corresponding to the intersection of vortex rings with the refocused plane at Z = 0.

Conclusion
Our purpose is to measure flow in a large-scale setup (a swimming pool) using bubbles as tracer and ambient lighting. Depth selectivity was achieved with synthetic aperture refocusing. The design of the setup, including magnification, bubble size, number of bubbles just matched the requirements of PIV. However, in future work, the bubble density should be increased substantially, despite the adverse effect of occlusion. Failing depth resolution and unsuccessful detection of velocity structures are caused by the emergence of ghost bubbles. The discrimination of ghost bubbles crucially depends on bubble illumination. A large illumination threshold kills ghost bubbles, but too few refocused bubbles may remain to analyze the velocity field. At present, this discrimination cannot be done effectively in the experiment. Depth discrimination of the bubble motion would greatly benefit from the addition of an illumination system. Another way to kill ghost bubbles is to track the path of bubbles in time. What is a ghost bubble at one instant is no longer a ghost bubble a while later (Elsinga and Tokgoz 2014;Elsinga et al. 2011). This mechanism also applies to PIV, which correlates refocused bubble planes at two instants. However, in PIV, it is only effective for largescale velocity gradients, and not for velocity structures buried in a bubble field.