Background oriented schlieren technique with fast Fourier demodulation for measuring large density-gradient fields of fluids

In order to measure the large density-gradient fields of fluids such as underwater shock waves, we employed a fast Fourier demodulation called Fast Checkerboard Demodulation (FCD) method in Background Oriented Schlieren (BOS) technique. BOS is a simple image-based measurement technique that detects the apparent displacement (local distortion) of a background image caused by the density gradient of the fluid in front of the background. The cross-correlation particle image velocimetry (PIV) method, which is commonly employed in the BOS technique for detecting the apparent displacement, uses a random-dot background, whereas FCD uses a periodic pattern as the background (e.g., checkerboard pattern, lattice grid pattern (grid scale)) to detect the apparent displacement as the phase change of the pattern in the Fourier space. In this study, we measured the apparent displacement, which is proportional to density gradient of fluid, of an underwater shock wave using FCD and PIV. The results showed that FCD can measure a displacement gradient of up to 2.5 times larger than that which PIV can measure. Furthermore, we systematically investigated the measurement limit of FCD-BOS by changing several parameters of the periodic patterns, such as the grid size. Also, we explored the related parameters in the Fourier space to understand the limitations. It is worth noting that FCD-BOS with lattice grid pattern (grid scale) can measure the apparent displacement as accurately as that with a custom-made pattern, indicating that large density-gradient fields of fluids can be measured using a simple setup with a commercially available (inexpensive) pattern.


Introduction
Background Oriented Schlieren (BOS) technique is a non-contact method for measuring the density field in a fluid [1,2,3,4,5]. Compared to conventional Schlieren imagin and optical interferometry, the measurement system is quite simple: it consists of only a background, camera, and light source. BOS technique has been applied to various density fields, such as gas flows [6,7,8,9,10], shock waves [4,5,11,12], flames [13,14,15,16,17], and liquid-gas interfaces [18,19]. Furthermore, it has also been applied to various length scales, from micrometers [2,11,12,20] to meters [5,21]. Moreover, complex 3D density fields have been also measured using multi-directional BOS imaging [22,15,23]. Figure 1 shows the illustration of BOS technique. BOS technique quantifies the density-gradient field in front of a background by comparing two background images with and without the density gradient which are captured using a camera. The background image with density gradient looks distorted as compared to the background image without the den-sity gradient because the refractive index of a fluid changes with density. Such distortion is quantified as "apparent displacement". It is extremely important for the BOS technique to accurately measure the apparent displacement between the background images [24,25,26]. A random-dot background is commonly used to obtain apparent displacement of images through analysis by digital image correlation (DIC), which is known as particle image velocimetry (PIV) in the fluid dynamics research community [27,28,29,30] or optical flow [31,32,33]. However, neither PIV nor optical flow is capable of detecting displacements of large density gradients because they assume parallel displacement of the background image, while in reality, random dots experience a shape distortion due to large density gradients as illustrated in Figure 2(a).
Let us consider a simple situation in which a square ∆x × ∆x on a background is distorted in one direction (Figure 2(b)). Because the magnitudes of the displacements u 1 on the left side of the square and u 2 on the right side of the square are different, the size of the square after distortion is expressed as ∆x = ∆x+u 2 −u 1 , and the rate of shape distortion is (∆x − ∆x)/∆x = (u 1 − u 2 )/∆x ≈ ∂u/∂x. Thus, the magnitude of the distortion depends on the displacement gradient ∂u/∂x, i.e., the second-order spatial differential of the density field of the fluid. When the displacement gradient is too large, PIV cannot properly acquire the apparent displacement [34]. Hence, there is a limitation in the existing PIV-based BOS measurements for underwater shock waves and flames [35] due to their large displacement gradients.
In this study, to solve this problem fundamentally, we employed the Fast Checkerboard Demodulation (FCD) method which was proposed by Wildeman [36], because unlike PIV, FCD can measure the displacement even when the background pattern is distorted. FCD uses a background with a twodimensional periodic pattern to detect the displacement field from the change of wavenumber in the Fourier space. By using FCD, Wildeman [36] has successfully detected waves on a water surface. This study takes advantage on the robustness of FCD in detecting large apparent displacements. It is the first time that FCD method is used in a BOS technique. In our measurement system, the background patterns of horizontal and vertical lines (lattice grid pattern) and of a black-and-white square pattern (checker pattern) are used. For comparison, a density-gradient field was measured by both PIV and FCD. Also, the data measured using hydrophone was used as the reference (correct) data. We measured underwater shock waves since they are good examples of flow field with large displacement gradients. We also performed numerical validation by adding the displacement to a synthetic background image to gain a better understanding on the experimental values and the measurement limits of FCD-BOS. This paper is structured as follows. Principles of BOS technique and FCD method are described in Section 2. The experimental setup and results are ex- Figure 2: (a) (Left) A background image of dot pattern without apparent displacement (Reference image). (Right) A background image with apparent displacement (Distorted image). Black squares are not only translated but also deformed, causing problems to the displacement detection by PIV. (b) Simple situation under a lateral apparent displacements u 1 and u 2 . The magnitude of the dot deformation is represented by the displacement gradient ∂u/∂x. plained in Section 3. We compared the apparent displacement and displacement gradient obtained using FCD and PIV with those estimated from the pressure measured by the hydrophone. Finally, we also analyzed the synthetic image under ideal conditions for discussion on the measurement limit of FCD-BOS in Section 4.

Background Oriented Schlieren (BOS) technique
BOS is a technique which visualizes and quantifies the density gradient by detecting variations in the refractive index of the fluid. As shown by the illustration in Figure 1, the measurement target is placed between a camera and a background where the camera is focused onto the background to capture images with and without density gradients. FCD or PIV is then used to determine the apparent displacement by comparing both images. The apparent displacement in y-direction, v, (Figure 1) can be related to the gradient of the refractive index as shown in the following equation [5] : where Z D is the distance from the background to the center of the density gradient, n is the refractive index of the measurement target, and n 0 is the refractive index of the ambient fluid. The relationship between the refractive index n and the fluid density ρ is given by the Gladstone-Dale equation: where K is the Gladstone-Dale constant, which is 3.34 × 10 −4 m 3 /kg for water [5]. Substituting Eq.
where ρ 0 is the density of the ambient fluid. Similarly, the apparent displacement in x-direction, u, is In this study, the apparent displacement measured by the BOS technique and that estimated from direct measurement of pressure using a hydrophone are compared to validate the ability of FCD-BOS in detecting apparent displacement. The procedure for estimating the apparent displacement from the pressure measured by the hydrophone is explained in § 2.3

Fast Fourier demodulation of a periodic background
A spatial Fourier demodulation technique of a periodic background, which is known as Fast Checkerboard Demodulation (FCD) method and was proposed by Wildeman [36], allows the quantitative extraction of the apparent displacement field when a periodic checkered background is used. By photographing a background with a periodic pattern (a lattice grid pattern or a checkered pattern in this study) and performing Fourier transform, the apparent displacement on the background is obtained as the phase change. Reference image I 0 (ξ) is expressed as a general 2D periodic pattern as where ξ denotes the pixel coordinates (x, y) in the image, k 1 and k 2 denote the orthogonal wavenumber vectors extracted from the reference image (called carrier peaks), and a mn denotes the Fourier coefficient. The distorted image I(ξ) is the result when the reference image I 0 is modulated by the apparent displacement u(ξ), as follows: By performing Fourier transforming on I 0 and I, we obtain the respective power spectras P 0 and P with peaks at mk 1 + nk 2 , which is expressed using integers m, n, as shown in Figure 3. By extracting the signals within a radius of k rad centered from k 1 and k 2 in P 0 and P , respectively, and performing the inverse Fourier transform, we obtain where a c1 , a c2 , a c1 , and a c2 are coefficients derived from a mn . The phase fields φ 1 (ξ), φ 2 (ξ) extracted from two linearly independent carrier peaks are or where k 1x is the x-component of vector k 1 , k 1y is the y-component of vector k 1 , k 2x is the x-component of vector k 2 , k 2y is the y-component of vector k 2 , and u(ξ) = (u, v). By solving Eq. 13, we obtain the displacement u, v directly as To avoid aliasing and phase wrapping, following three criteria must be satisfied for the displacement measurement [36]: where k s is the wavenumber of the displacement, k c ( k c = |k 1 | = |k 2 | ) is the norm of the carrier peak of the reference image, and u is the displacement gradient. In this study, we focus on the value on xaxis. For this reason, although Wildeman [36] uses k c in the criteria, we utilized the norm of carrier peak A distorted image I (1500 × 750 pixels) was created by displacing the reference image I 0 . P 0 and P are the power spectra in the Fourier space (k-space) of I 0 and I, respectively. k 1 and k 2 are the carrier peaks. The radius of a circle in P 0 is k rad . k cx is the norm of carrier peak k c along x-axis. The synthetic displacement (true) and the displacement detected by the FCD are shown. along x-axis, i.e., k cx . Thus, the criteria correspond to Eq. 16 and Eq. 17 in this study are respectively. Eq. 15 is a criterion that varies with the displacement and the background pattern, while Eq. 18 and Eq. 19 are criteria which vary only with the background pattern. For example, in the k-space at P 0 in Figure 3, k cx = k 1x = k rad for the checkered pattern and k cx = k 2x = 2k rad for the lattice grid pattern. If the criterion in Eq.15 is not satisfied, aliasing artifacts will occur; if the criterion in Eq.18 is not satisfied, overlap will occur; and if the criterion in Eq.19 is not satisfied, phase wrapping will occur [36].

Estimation of the apparent displacement from the pressure measured by the hydrophone
The apparent displacement was estimated from the pressure measured by the hydrophone for comparison with FCD and PIV. First, the density ρ was obtained from the pressure using Tait equation: where the hydrostatic density ρ 0 is 998 kg/m 3 , B and α are 314 MPa and 7.15, respectively. Then, the refractive index gradient ∂n/∂r was calculated by differentiating the refractive index n, which was obtained from Eq. 2 using the density ρ. Because the shape and pressure of a laser-induced underwater shock wave are axisymmetric in the direction of the laser beam [37,38], the axisymmetric field of the refractive index gradient ∂n/∂x, as shown in Figure  4, was reconstructed. Finally based on Eq. 1, the apparent displacement u was obtained by integrating ∂n/∂r along z-axis as shown in Fig. 4. The shock wave is axisymmetric about the direction of the laser beam (dotted straight line). (b) Constructed refractive index gradient field ∂n/∂x. To obtain u, integration of the ∂n/∂x field was performed toward the z-direction in the pink area of (a).

Experimental set-up
A schematic of the experimental setup is shown in Figure 5.   To determine the apparent displacement on the background, an open-source PIV code (PIVlab, MATLAB, Thielicke and Stamhuis [39]) and an opensource FCD code [36] were used. A dot occupies about 5.5 -11.3 pixel on the image. The PIV setting has a 50% overlap with the fast Fourier transform (FFT) multipass interrogation (from 64 to 8 pixels).
To compare with the displacements estimated by the hydrophone measurement, the displacements obtained by the BOS in the range of ±0.125 mm in the y-direction from the center-line of the hydrophone were averaged along the x-axis (see Figure 6). Note that the positions of the shock wave measured by the BOS and hydrophone in x-direction are slightly different, as shown in Figure 6. The pressure is inversely proportional to the distance of the shock wave propagation in this experimental range [40]. Thus, for direct comparison with BOS results, the shock wave pressure measured by hydrophone was corrected us- Figure 6: A typical background image with a shock wave (left). x-components of the displacement field detected by FCD (right). The analyzed area was determined by considering the effective diameter of the hydrophone.
ing following equation: where p image is the pressure at the time of imaging, p hyd is the pressure measured by the hydrophone, r image is the radius of the shock wave at the time of imaging, and r hyd is the distance between the center of the shock wave and the tip of the hydrophone. In Figure 7, the results of the FCD and the displacement calculated from the pressure measured by hydrophone are shown. In following sections, we focus on the largest apparent displacement at the distance near 2.7 mm (first peak). Note that a second peak (distance around 1.2 mm) appears in the estimated displacement because a part of the shock wave was reflected on the background. The second peak in the image analysis does not appear because the integral of ∂n/∂x is very small.

Experimental results
The maximum displacement detected by PIV or FCD, u max (PIV/FCD) and maximum displacements calculated from the pressure, u max (hydrophone), are shown in Figure 8. For u max (hydrophone) > 4 µm, the maximum displacement detected by PIV is about 4 µm, which is significantly lower than that estimated from the measurements of the hydrophone. Such underestimation is due to a significant distortion of the shape of the dots when u max (hydrophone) > 4 µm, making it difficult for PIV to accurately detect the displacement. In contrast, the displacements obtained by the FCD reasonably agrees with the estimated displacement.
The error of the PIV and FCD methods is defined as the difference between the maximum displacement detected using each method and that estimated from the pressure measured using hydrophone. In Figure  9, the errors are plotted against the maximum displacement gradient (∂u/∂x) max , which represents the magnitude of the distortion. As shown in the figure, PIV underestimates the displacement when the maximum displacement gradient is greater than 0.20. This is because the displacement detected by the crosscorrelation method does not match the original displacement owing to the large distortion of the dots as discussed in Section 1. On the other hand, FCD can reasonably detect the maximum displacement for the displacement gradient ∂u/∂x < 0.5, which is approximately 2.5 times larger than 0.20, the measurement range of PIV-BOS. Among the FCD methods, the error of 20 µm checkers is less than the error of 40 µm checkers and 20 µm lattice grid. Besides, it is also observed that 20 µm lattice grid, the error shifts from overestimation to underestimation at the displacement gradient between 0.3 and 0.4. To analyze the reasons for such trends, we utilized synthetic images to investigate the applicable range for displacement detection of FCD with various background patterns. Synthetic images were used because the images captured from the experiment might be affected by random errors such as image blurring and slight misalignment of the experimental setup. The analysis are discussed in Chapter 4. 4 Numerical approach: measurement limit of FCD-BOS

Numerical condition
In order to determine the applicable range of displacement detection of FCD, we used synthetic images (image size: 1500 × 750 pixels, spatial resolution: 1 µm/pixel) that imitate three types of background patterns: random dots, lattice grid, and checkered. The synthetic images were modulated with synthesized apparent displacement u(r), i.e., "true displacement", based on the Gaussian distribution as follows: where u max varies between 1 -11 µm at an interval of 1 µm, r is the distance from the origin indicated in Figure 10, r c is set to 400 µm, and β is set to 400×10 −12 m 2 . This u(r) resembles the displacement of a shock wave. For each interval of u max , the origin of u(r) was issued to 20 different positions relative to Figure 9: The error which is the difference between the maximum displacement estimated from the pressure u max (hydrophone) and the maximum displacement detected by PIV/FCD u max (PIV/FCD) against the maximum displacement gradient of the displacement estimated from the pressure measured by hydrophone (∂u/∂x) max . Error bars indicate the standard deviations.
the background pattern. These 20 data were used to obtain the averaged value and its error range. Figure 10 shows the typical result for u max = 5 µm. Overall, the results show similar trends as the experimental results. The displacement detected by PIV is noisy, while those detected by the FCD methods are similar to the true displacement. However, the measurement accuracy of FCD depends on the background pattern where the displacement fields of 20 µm checkers, 20 µm (0 • ) lattice grid, and 20 µm (45 • ) lattice grid are almost identical to the true data, but the displacement field of 40 µm checkers is partially different from true data. The reason for this is discussed in the later part of this section. For the lattice grid patterns without rotation (0 • ) and with a rotation of 45 • , k cx are 0.32 pixel −1 and 0.23 pixel −1 , respectively although k rad are 0.16 pixel −1 for both patterns. Nevertheless, since the lattice grid patterns of both 0 • and 45 • satisfy Eq. 19, there is no notable difference between their displacement fields. Relations between the maximum displacement gradient and error are shown in Figure 11. PIV underestimates the apparent displacement as the displace- Figure 10: Synthetic images (1500 × 750 pixels) and the x-direction displacement u (maximum displacement is 5 µm (= 5 pixels) in this figure). The radius of the circle in each P 0 is k rad . Symbols are explained in the caption of Figure 3. Figure 11: The maximum displacement gradient of the displacement given modulated to the reference image (∂u/∂x) max and the error which is the difference between the maximum displacement given modulated to the reference image u max (estimation) and the maximum displacement detected by PIV/FCD u max (PIV/FCD). Error bars indicate the standard deviations estimated by using 20 sets of images. ment gradient increases, similar to the experimental results shown in Figure 9. On the other hand, the measurement error of FCD (20 µm checkers) is sufficiently low for (∂u/∂x) max < 0.5 while that of FCD (20 µm lattice grid) has non-trivial values for (∂u/∂x) max > 0.45. Since the value of k cx for the lattice grid pattern (20 µm) is 0.32 pixel −1 , based on Eq. 19, phase wrapping occurs when the displacement is larger than 9.8 µm, which corresponds to a displacement gradient of 0.42. Thus, FCD (20 µm lattice grid) underestimates the displacement when (∂u/∂x) max > 0.45. The FCD (40 µm checkers) gradually underestimates the displacement as the displacement gradient increases because of the resulting aliasing artifacts as described by Eq. 15.

Size and orientation of the background patterns
To investigate further how the criterion in Eq.15 offsets the results in our condition, the power spectra P d , which is obtained by the Fourier transformation of the true displacement u(r), is analyzed. Since the analysis results are similar for every displacement between u max = 1 and 11 µm, here only the results for the synthetic image of u max = 1 µm are discussed. Figure 12(a) shows the P d for u max = 1 µm, where the red dashed circle in the zoomed-in figure shows the boundary of k rad = 0.26 pixel −1 . Since Eq.15 shows that the measurement accuracy correlates to a ratio of the sum of P d contained within k rad to the total sum of P d , the ratio is plotted as a function of k rad in Figure 12(b). In the same figure, dashed vertical lines are Figure 12: (a) Power spectra P d , which is obtained by the Fourier transformation of the true displacement u(r), for u max = 1 µm. The subfigure on the upper right shows the enlarged view. The radius of the circle, k rad , varies depending on the pattern. (b) The ratio of the sum of P d contained within k rad to the total sum of P d .
plotted for the k rad values of the patterns of different sizes, i.e., k rad = 0.08 pixel −1 (40 µm checkers and lattice grid), 0.16 pixel −1 (20 µm checkers and lattice grid), and 0.26 pixel −1 (12 µm checkers and lattice grid). The ratios for k rad = 0.08 pixel −1 , 0.16 pixel −1 , and 0.26 pixel −1 are 65%, 87%, and 95%, respectively. Since ratio for 40 µm checkers is lower that that for 20 µm checkers, the resulting aliasing artifacts in 40 µm checkers have caused more noise in the detected displacements than those of 20 µm checkers, as shown in Figure 10. Besides, they have also caused the underestimation of the displacement by 40 µm as shown in Figure 11. Such results show that the measurement accuracy can be further improved when the ratio is increased through increasing the values of k rad using the patterns of smaller size.
Here, we focus on k rad = 0.26 pixel −1 (12 µm checkers and lattice grid) which correspond to a ratio of 95%. For that, the measurement accuracy of 12 µm checkers and lattice grid were evaluated based on the respective measurement errors shown in Figure 13 and Figure 14.
In summary, the limiting factors for FCD can be understood from the parameters in Fourier space as shown in Eq. 15, Eq. 18, and Eq. 19. It is worth noting from a practical view point that the results showed that FCD of lattice grid patterns can measure displacement as accurate as that of checkered patterns, although it is slightly more restricted.

Line thickness of lattice grid pattern
For practical purpose, we investigated the influence of line thickness of the lattice grid pattern. The results of the different line thicknesses of 12 µm lattice grid patterns for the maximum displacement u max -3.14 0 3.14 -3.14 0 3.14 -3.14 0 3.14 -3.14 0 3.14 -3.14 0 3.14 -3.14 0 3.14 -3.14 0 3.14 -3.14 0 3.14 14 0 3.14 -3.14 0 3.14 -3.14 0 3.14 -3.14 0 3.14 π −π −π π 0 0 ky [1/pixel] kx [1/pixel] -3.14 0 3.14 -3.14 0 3.14 -3.14 0 3.14  = 5 µm are shown in Figure 15. In these cases, all FCD criteria (Eq. 15, Eq. 18, and Eq. 19) are satisfied. While k cx and k rad are the same for all cases, the wavenumber intensity P of the synthetic images depends on the line thickness. If the lines are significantly thicker than the white regions, as shown in Figure 15(d), the u field is different from the true displacement described in Eq. 22. Since u fields in Figure 15(a)(b) are identical to true displacements, it is recommended to use lattice grid patterns with lines that are thinner than the white regions for FCD-BOS measurements.

Conclusion
We applied Fast Checkerboard Demodulation (FCD) method for Background Oriented Schlieren (BOS) technique. In general, BOS tecnique uses a randomdot background to detect the displacement. However, to measure large density-gradient fluid fields such as shock waves, the dotted pattern distorts caus-ing the displacement detection by cross-correlation PIV method fails to work properly. From the experimental results, FCD is superior to PIV in terms of detecting large displacement gradients and the maximum measurable displacement gradient of FCD is 2.5 times larger than that of PIV. Through analysis using synthetic data, the limiting factors of the FCD can be understood from the parameters in Fourier space as shown in Eq.15, Eq.18, and Eq.19. The results of the analysis shows that line thickness influences the limiting factors of the FCD of lattice grid patterns. On a final note, FCD-BOS can be easily constructed and customized for various thermal fluids since it just consists of a commercially available grid scale, still camera and light source. Based on measurement limitation of FCD-BOS reported in this paper, users of FCD-BOS could select proper background with suitable parameters such as the line width of lattice grid patterns and region size of checker patterns.