Optical-flow-based background-oriented schlieren technique for measuring a laser-induced underwater shock wave

The background-oriented schlieren (BOS) technique with the physics-based optical flow method (OF-BOS) is developed for measuring the pressure field of a laser-induced underwater shock wave. Compared to BOS with the conventional cross-correlation method in PIV (called PIV-BOS), by using the OF-BOS, the displacement field generated by the small density gradient in water can be obtained at the spatial resolution of one vector per pixel. The corresponding density and pressure fields can be further extracted. It is particularly demonstrated that the sufficiently high spatial resolution of the extracted displacement vector field is required in the tomographic reconstruction to correctly infer the pressure field of the spherical underwater shock wave. The capability of the OF-BOS is critically evaluated based on synchronized hydrophone measurements. Special emphasis is placed on direct comparison between the OF-BOS with the PIV-BOS.

displacement field induced by the density change in fluid, which is used to infer the density and pressure fields. The BOS technique has been used in various measurements particularly in compressible air flows and gas-mixing flows (Raffel et al (2000); Meier (2002); Venkatakrishnan and Meier (2004); Murphy and Adrian (2011)). A review of the BOS technique is given by Raffel (2015). The BOS technique is an extension of the classical schlieren technique traditionally used in high-speed wind tunnel for aerodynamics measurements. Due to the significant advance of digital cameras and image processing, the experimental setup of the BOS system is much simpler than its classical counterpart with complex arrangements of optics. In general, a background pattern plate placed behind a measurement domain of fluid is imaged by a camera in the cases with and without the fluid density gradient. The background pattern image is disturbed by the density gradient that deflects the light rays radiated from the background pattern plate. The disturbed image appears shifted relative to the undisturbed image taken in the case with the homogenous fluid density. The displacement field in the disturbed image is related to the path integral of the density gradient. When the displacement field is measured, the density and pressure fields could be inferred. To determine the displacement field from the undisturbed and disturbed images, the cross-correlation method in particle image velocimetry (PIV) can be applied.
However, the application of the BOS technique in liquid is more difficult compared to its application in gas flows since the light-ray deflection due to the density change in liquids is much smaller. In particular, the measurement of an underwater shock wave by using the BOS technique is scanty. Recently, using the BOS technique with the cross-correlation method, we obtained the displacement field induced by the local density gradient of the underwater shock wave (Yamamoto et al (2015)). Based on these results, we tried to further reconstruct the pressure field from the displacement field. However, it is found that the pressure field could not be correctly reconstructed due to the insufficient number of the displacement vectors extracted by the typical correlation algorithm in PIV (Hayasaka et al (2016)). In order to overcome this problem, we propose a high-resolution BOS technique by incorporating the optical flow method, which is simply referred to as OF-BOS. Atcheson et al (2009) suggested that the typical optical flow methods in computer vision, such as the Horn-Schunck method and the Locas-Kanade method, could be used for the BOS technique. They found that these optical flow methods particularly the Horn-Schunck method performed better than the correlation method in terms of the accuracy and the spatial resolution.
In this paper, we adopt the physics-based optical flow method in BOS, which is developed by Liu and Shen (2008) for various flow visualizations. This method could achieve the spatial resolution of one vector per pixel and the better accuracy (Liu and Shen (2008)), which is important for the application of the BOS technique in water.
Following this instruction, we will describe the optical-flow-based BOS technique, tomographic reconstruction of the field of the divergence of the displacement vector, and determination of the density field by solving the Poisson's equation. Then, the experimental setup will be described for BOS measurements of the laserinduced shock wave. Next, the results will be discussed, focusing on the effects of the spatial resolution of the extracted displacement field on reconstruction of the density and pressure fields. In particular, we will quantitatively compare the pressure distributions obtained from both the OF-BOS and the PIV-BOS with hydrophone data.

BOS
The basic equations in the BOS technique are briefly recapitulated for convenience of reading. The working principle of the BOS technique is illustrated in Fig. 1 (Venkatakrishnan and Meier (2004); Yamamoto et al (2015)). A typical BOS system consists of a camera, a background pattern plate (such as a random dot pattern plate), a light source, and a computer for image/data processing. For simplicity, it is assumed that the background reference plane (x, z) is parallel to the image plane. Thus, the image coordinates are equal to λ(x, z), where λ is a proportional constant in the orthographical projection. The light ray radiated from the background pattern plate is deflected from the original straight path through a fluid domain with the density gradient due to optical refraction, and therefore the image of the background pattern has an apparent displacement (shift) field.
The displacement on the image plane is proportional to the path integral of the small density gradient through a fluid domain. For BOS setting, van Hinsberg and Rösgen (2014) gave the following relation for the in-plane displacement vector w in the background plate where c (100 µm in this case) is the thickness of the density gradient domain, l 0 (8 mm in this case) is the distance from the density gradient domain to the background plate, n 0 is the refractive index of water at the condition without a shock wave, n is the refractive index of water in the test condition, and ∇ is the gradient operator on the coordinate plane (x, z). The relation between the refractive index and the fluid density is given by the Gladstone-Dale equation (Merzkirch (2012); Raffel (2015))  where K is the Gladstone-Dale constant (3.34×10 −4 m 3 /kg for water), and ρ is the density of the fluid. Substitution of Eq.
(2) into Eq. (1) yields where ρ 0 is the fluid density under hydrostatic pressure. It is emphasized that the displacement vector w obtained by the BOS technique is the path-integrated (or projected) quantity across the measurement domain with the density gradient.
In general, to extract the 3D field from the path-integrated quantity requires tomographic reconstruction from data at multiple viewing directions. Further, applying the dot product of ∇ to Eq. (3), we have the Poisson's equation for the density where the source term S is proportional to the divergence of the displacement vector w by In principle, when the displacement vector w is measured, the density field can be determined by solving Eq. (4). The Gauss-Seidel method is used to solve Eq. (4) numerically here. Then, the pressure field can be further determined by applying the Tait equation where p 0 is the hydrostatic pressure, B is a constant of 314 MPa, and the exponent α is 7 (Brujan (2010); Yamamoto et al (2015)). In this study, we measure a spherical shock wave in water with the peak overpressure of about 1 MPa. In this case, the density change due to the overpressure is less than 1 kg/m 3 . This means that the change in the refractive index is very small [O(10 −4 )].
The BOS system acquires an undisturbed background pattern image and the corresponding disturbed image. Figure 2 shows the typical images obtained in the experimental setup in this work. Figure 2(i) shows a shadowgraph snapshot of a shock wave propagating spherically around a laser-induced bubble, where the high-pressure region is expected at a shock front. Figure 2(ii) and 2(iii) show an undisturbed background reference image and the corresponding disturbed image, respectively. The zoomed-in image of the particle pattern is shown as well. Without image processing, it is difficult to see directly the shock wave in Fig. 2(iii) since the change of the refractive index of water is so small.

Optical flow method
The key element in the BOS technique is to determine the displacement vector field in the undisturbed background reference image and the corresponding disturbed image. This computation is usually carried out by using the crosscorrelation method in PIV. An effort has been made to adopt the classical optical flow methods in computer vision science in BOS (Atcheson et al (2009)). The optical flow method can be further developed based on rational physical foundations for various flow visualizations (Liu and Shen (2008)). The optical flow is described by the generic equation in the image plane, i.e., where u is the velocity in the image plane referred to as the optical flow, g is the normalized image intensity, ∇ is the spatial gradient in the image plane, and f is a term related to the diffusion and boundary fluxes which are negligibly small in most cases. In velocimetry, the optical flow in Eq. (7) has a clear physical meaning, that is, the optical flow is proportional to the light-ray-path-averaged velocity of fluid or particles in flows. We call this method as the physics-based optical flow technique. In a special case where ∇·u = 0 and f = 0, Eq. (7) is reduced to the Horn-Schunck brightness constraint equation which is the foundation of the classical optical flow method developed by Horn and Schunck (1981). It is noted that in a general case the optical flow is not divergence-free and ∇·u = 0 is not physically true. For BOS applications, the difference ∂g/∂t ≈ ∆g/∆t is used in Eq. (7), where ∆g = g -g ref is the difference between the disturbed image and the undisturbed reference image, and the nominal time interval is unitary (∆t = 1). Therefore, the optical flow in Eq. (7) is interpreted as the displacement vector field in the image plane (i.e. u = w) that is generated by the deflected light ray through the density gradient domain. The displacement vector in the background plane is related to that in the image plane by w = λw .
To determine the displacement vector field in the image plane, a variational formulation with a smoothness constraint is typically used (Liu and Shen (2008)). By minimizing the functional, the Euler-Lagrange equation is given for the optical flow. The standard finite difference method is used to solve the Euler-Lagrange equation with the Neumann condition on the image domain boundary for the optical flow. The optical flow algorithm used in this work has the routines: the Horn-Schunck estimator for an initial solution (Horn and Schunck (1981)) and the Liu-Shen estimator for a refined solution of Eq. (7) (Liu and Shen (2008)). The relevant parameters in pre-processing and optical flow computation should be suitably selected. The main parameters are the Lagrange multipliers for the Horn-Schunck and Liu-Shen estimators. Other parameters are the number of iterations in successive improvement of optical flow computation by using a coarse-to-fine iterative scheme, and the sizes of the Gaussian filters for correcting the effect of a local illumination intensity change and removing small random noise in images. A mathematical analysis of the physics-based optical flow and an iterative numerical algorithm are given by Wang et al (2015). Quantitative comparison between the optical flow and cross-correlation methods for PIV images has been evaluated by Liu et al (2015). The variational optical flow method based on Eq. (7) is a differential approach that can achieve the theoretical spatial resolution of one vector per pixel. Figure 3 illustrates the difference between the optical flow method as a differential approach and the cross-correlation method as a region-based integral approach. The optical flow method is particularly suitable to detect small displacements in narrow regions (such as shock wave).

Tomographic reconstruction on a spherical surface
The tomographic reconstruction technique is used to reconstruct the relevant physical quantity (such as the fluid density or the divergence of the displacement vector) on a spherical surface of a propagating shock wave from the path-integrated or projected quantity obtained from the BOS measurements. In volumetric flow visualizations, an image can be modeled as a set of 1D projections of a field of a quantity with the given projection angles through a 2D domain as illustrated in Fig. 4. The basic tomographic problem is how to reconstruct the 2D field from these projections. The projection of the quantity q(x, y) can be expressed as an integral along a ray, i.e., where l = x cos θ + y sin θ is the coordinate along a projection line, θ is the angle defining the projection lines, R(q) denotes the Radon transform, and δ denotes the Dirac-delta function. An inversion of the Radon transform is sought for q(x, y). This tomographic problem has been studied in 3D flow measurements (Feng et al (2002); Venkatakrishnan and Meier (2004)). We introduce the Fourier transform Then, the solution of the integral equation Eq. (8) can be formally given by Clearly, the tomographic reconstruction of q(x, y) requires many projections in θ ∈ [0, π]. In actual computation, the inverse Radon transform in Matlab is used in this work. In this process, we generate the sinogram digitally. Although the Abel transform is simpler for the reconstruction of an axial symmetrical field, we do not utilize it because of its high sensitivity to noise (Venkatakrishnan and Meier (2004)).
In the BOS measurement of a spherical shock wave, the following procedures are proposed as illustrated in Fig. 5. 1. Image registration based on the affine transformation is applied to the reference and disturbed images to correct any global misalignment between them caused by the possible movement of the camera and other factors during measurements. 2. The projected quantity is selected as the divergence of the projected displacement vector, i.e., Q(l, θ) = ∇ · w from the BOS measurement since it is the source term of the Poisson's equation for the fluid density. The quantity ∇ · w is mapped onto the image plane based on an assumption of the axial symmetrical structure of the laser-induced shock wave. This is illustrated in Step (1-2) in Fig. 5.  3. The Radon transform Q(l, θ) = ∇ · w is given as a sinogram at section at a given z-coordinate, as illustrated in Step (2-3). In the sinogram, we generate 180 projected data in θ ∈ [0, π]. 4. The distribution of q(x, y) at that section is reconstructed by using the inverse Radon transform (symbolically expressed as q(x, y) = R −1 (Q)), as illustrated in Step (3-4). In Matlab computation, we apply the spline interpolation to the back projection method and calculate the inverse Radon transform without filtering. 5. The field of ∇ · w on the spherical shock wave is obtained by superposition of the reconstructed fields over a set of cross sections in the z-coordinate, as illustrated in Step (4-5). This superposition procedure for a spherical shock wave is further illustrated in Fig. 6. 6. The field of the density is obtained by solving the axisymmetric Poisson's equation for a given source term on the plane of symmetry and the corresponding pressure field is calculated by using Eq. (6), as shown in Step (5-6).
It is necessary to comment the selection of the projected quantity for tomographic reconstruction. In the BOS measurements by Venkatakrishnan and Meier (2004), the projected fluid density was obtained first by solving the Poisson's equation, and then the local density was reconstructed by using the tomographic technique. In contrast, the present procedures conduct the tomographic reconstruction before solving the Poisson's equation. The projected quantity is the divergence of the projected displacement vector, i.e., Q(l, θ) = ∇ · w , and thus the divergence of the local displacement vector is first reconstructed by using the tomographic technique. Next, the local density is calculated by solving the Poisson's equation. In principle, since the Poisson's equation is valid for both the local and path-integrated (projected) quantities, the approach used by Venkatakrishnan and Meier (2004) should be equivalent to the present approach. Nevertheless, the rationale behind the present arrangement is to avoid the propagation of the error in solving the Poisson's equation into tomographic reconstruction that is more sensitive to the errors.  Ltd.) was focused through a 20×microscope objective lens to a point inside a distilled-water-filled glass tank (450×300×300 mm 3 ). The underwater shock wave was generated at a laser-focused point, propagating spherically. The background random-dot-pattern plate was placed behind the laser-focused point inside the tank. The shock wave was recorded by using a camera with 2048×2048 pixels (FASTCAM Mini WX50, Photron Ltd.). The spatial resolution was set at 5.18 µm per pixel. The light source utilized for illuminating the background pattern was a laser stroboscope (SI-LUX 640, Specialized Imaging Ltd.) with the pulse width of 20 ns. The pulsed laser, camera, and light source were synchronized by using a delay function generator (Model 575, BNC Co.). We utilized a PVDF pressure sensor (i.e. hydrophone) (Muller-Platte Needle Prove, Mueller Instruments) to validate the results obtained by using the OF-BOS and the PIV-BOS. The pressure sensor was set toward a center of the laser-induced shock. The distance from the center to the hydrophone was 5.0±0.1 mm. In the PIV-BOS, the open-source Matlab PIV code called PIVlab is used to determine the displacement vectors, which is described by Thielicke and Stamhuis (2014).

Typical case
A typical case with the laser power of 1.2 mJ is considered to compare the results obtained by using the OF-BOS and the PIV-BOS. Figure 8 shows the displacement magnitude fields obtained by using the OF-BOS and the PIV-BOS. The OF-BOS gives 2008×1085 vectors while the PIV-BOS gives 250×138 vectors in the same domain in Fig. 8. In optical flow computations, the Lagrange multipliers for the Horn-Schunck estimator and the Liu-Shen estimator are 20 and 2000, respectively. In PIV computations, the window size is iteratively reduced from the initial size of 32×32 pixels to 16×16 pixels and then to 8×8 pixels. Both the OF-BOS and the PIV-BOS capture the sharp change induced by the shock wave in the displacement fields. As expected, the cross-correlation computations in windows in the PIV-BOS tend to smooth out sharp features like the shock waves. Figure 9 shows the profiles of the measured displacement along a ray aligned with the hydrophone, where the averaging operation over 50 lines for the OF-BOS and 20 lines for the PIV-BOS in the z-direction is applied to reduce the noise. The displacements were estimated from pressure data by using Eqs. (3)-(6) with an assumption of constant speed of propagation. It is indicated that both the OF-BOS and the PIV-BOS give results consistent with that given by the hydrophone at that location. In particular, the displacement induced by the shock wave is well captured by both the techniques. Nevertheless, the OF-BOS achieves a much higher spatial resolution. Figure 10 shows the reconstructed pressure fields obtained by using the OF-BOS and the PIV-BOS. The OF-BOS capture correctly the sharp pressure change across the shock wave, while the PIV-BOS has a larger error in the reconstructed pressure field particularly near the 'north pole' and 'south pole' due to its low spatial resolution. For quantitative comparison, as shown in Fig. 11, the pressure profiles reconstructed based on the displacement vector fields given by the OF-BOS (2008×1085 vectors) and PIV-BOS (250×138 vectors) are plotted against the data obtained by the hydrophone. The pressure profile given by the OF-BOS along a ray aligned with the hydrophone is consistent with that given by the hydrophone. In contrast, the pressure profile given by the PIV-BOS exhibits a considerably broader distribution extending to the inner region although its peak value at that location agrees with the data given by the hydrophone. Clearly, the OF-BOS provides the more accurate reconstruction of the pressure field of the shock wave. In contrast, the PIV-BOS has a much larger deviation from the Fig. 9 The profiles of the measured displacement along a ray aligned with the hydrophone, where R is the radius of the shock wave, where the averaging operation over 50 lines for the OF-BOS and 20 lines for the PIV-BOS in the z-direction is applied to reduce the noise. data given by the hydrophone, which is particularly contributed by the large errors near the 'north pole' and 'south pole' (as shown in Fig. 10). In this case, the major issue of the PIV-BOS is its much lower spatial resolution that tends to corrupt the tomographic reconstruction of the shock wave.

Effect of laser energy
Measurements at different levels of the laser power were conducted to further compare the data obtained by using the OF-BOS, the PIV-BOS and the hydrophone. Figure 12 shows the typical displacement distributions induced by the laser-induced shock wave at three levels of the laser energy, where the OF-BOS and PIV-BOS data are obtained along the ray aligned with the hydrophone. Overall, the OF-BOS and PIV-BOS data are in good agreement with the hydrophone data in the low and medium levels of the laser energy (0.6 and 1.4 mJ). In the cases of the higher levels of the laser energy (2.1 and 2.9 mJ) with larger displacements, a coarse-to-fine iterative scheme is adopted to improve the accuracy of optical flow computation (Liu et al (2012). In this scheme, images are initially downsampled by 2 for a coarse-grained velocity field and then a refined velocity field with the full image resolution is obtained in iterations for correction of the large displacements. Two and three iterations are applied to the cases of 2.1 and 2.9 mJ, respectively. Figure 13 shows the corresponding pressure distributions of the laser-induced shock wave at three levels of the laser energy. The OF-BOS is able to detect the weak shock wave and the pressure peaks of the shock wave are consistent with those given by the hydrophone. Figure 14 shows the peak pressure value of the shock wave as a function of the laser energy, indicating the favorable comparison between the data obtained by using the OF-BOS and the hydrophone. The PIV-BOS also detects the peak pressure in the lower energy cases, but the pressure distribution inside the shock front does not match that of hydrophone.

Effect of spatial resolution
It is conjectured that the spatial resolution of the extracted displacement field would have a significant impact on the tomographic reconstruction of the pressure field of the shock wave. This is because the spatial derivatives of the displacement vector field should be calculated accurately for the tomographic reconstruction of the divergence of ∇ · w , which will have the direct influences on the extracted density and pressure fields. To examine this problem, the displacement field of 2008×1085 vectors in the half-circle domain obtained by the OF-BOS is downsampled into several fields with the lower resolutions. Then, the tomographic reconstruction procedures are applied to these fields to observe the effect of the spatial resolution. Figure 15 shows the reconstructed pressure fields based on the selectively downsampled displacement fields, where the field obtained by the PIV-BOS is also shown for reference. It can be observed that the quality of the reconstructed pressure field is degraded as the spatial resolution decreases. The degraded result given by the PIV-BOS is particularly obvious.
The quality of the pressure field particularly near the shock wave directly depends on the spatial resolution of the reconstructed field of the divergence of the displacement field, which can be illustrated in Fig. 16. Furthermore, as shown in Fig. 17, the effect of the spatial resolution on the reconstructed pressure field can be quantitatively observed in the pressure profiles along the ray aligned with the hydrophone. For the OF-BOS, the displacement field of 2008×1085 vectors is downsampled to that of 268×148 vectors. As shown in Fig. 17(a) the reconstructed pressure distribution exhibits that the pressure peak is broaden increasingly with the decrease of the spatial resolution and the peak value deviates from the data Fig. 11 Comparisons between the pressure profiles obtained by using the OF-BOS, the PIV-BOS and the hydrophone, where R is the radius of the shock wave. Fig. 12 The profiles of the measured displacement at three levels of the laser energy, where the OF-BOS and the PIV-BOS data are obtained along the ray aligned with the hydrophone. given by the hydrophone. In addition, the displacement field obtained by the PIV-BOS can be interpolated to generate a pseudo high resolution field of 2091×1000 vectors. As shown in Fig. 17(b), such interpolation could indeed improve the quality of the reconstructed pressure field even though the peak value becomes lower.

Conclusion
In BOS measurements of a shock wave induced by a laser in water, it is critical to obtain the displacement vector field with the high spatial resolution. In this work, the physics-based optical flow method is incorporated into the BOS technique, simply called the OF-BOS, which can obtain the displacement vector field from BOS images at the theoretical resolution of one vector per pixel. The tomographic reconstruction for a spherical laser-induced shock wave is proposed, which  is applied to the fields of the projected divergence of the displacement vector at cross sections in the vertical direction. Superposition of the tomographic solutions at these cross sections yields the local field of the divergence of the displacement vector in the domain of the spherical shock wave. Then, the local density field is obtained by solving the Poisson's equation and further the corresponding pressure field is calculated. Fig. 16 The effect of the spatial resolution of the field of the divergence of the displacement vector on the tomographic reconstruction on the pressure field of the shock wave. Fig. 17 The effect of the spatial resolution of the displacement field on the pressure profiles along the ray aligned with the hydrophone: (a) the OF-BOS, and (b) the PIV-BOS, where R is the radius of the shock wave.
The OF-BOS gives the pressure field consistent with the result obtained by the hydrophone. In contrast, the BOS with the conventional cross-correlation method in PIV (simply called the PIV-BOS) fails to correctly reveal the shock wave in the reconstructed pressure field due to its much lower spatial resolution. It is found that the quality of the tomographic reconstruction of the divergence of the displacement vector field significantly depends on the spatial resolution of the extracted displacement vector field. Therefore, the reconstructed density and pressure fields of the spherical shock wave are directly affected by the spatial resolution. Clearly, the OF-BOS has the advantages over the PIV-BOS in this aspect.