Formulation and demonstrations of three-dimensional background-oriented schlieren using a mirror for near-wall density measurements

Three-dimensional background-oriented schlieren (3D-BOS) is an effective method for reconstructing 3D density fields from optically measured data, but it has limitations in measuring near-wall regions, where most of the light paths are blocked. This paper proposes a new extension, 3D-BOS using Mirror, which uses a wall as a mirror to provide sufficient light paths. In this paper, first, the conventional formulations are modified for the proposed method to handle the mirror reflections of the light paths. Subsequently, the proposed method is validated using artificially generated model data of an ideal axisymmetric distribution. The proposed method can reconstruct the distribution as accurately as the conventional method for all the number of cameras examined. Finally, the proposed method is experimentally demonstrated using a candle plume. The proposed method can capture cylindrical low-density regions near the wall surface.


Introduction
The schlieren method is a technique used to visualize light refraction caused by density gradients. It enables to observe the instantaneous flow geometry of a wide area without disturbing the flow and is an essential tool for studying compressible flows, such as jets, wakes, mixing layers, and shock waves (Van Dyke 1982;Settles 2001). A variant of the schlieren method, background-oriented schlieren (BOS) (Raffel 2015), has also been widely used. The patterns on the background behind the flows appeared to shift because of light refraction caused by density gradients. The magnitude of refraction can be obtained quantitatively by measuring the displacement of the patterns in an image captured with a camera. Although BOS has some disadvantages, such as blurred images of the flow owing to out-of-focus and inferior spatial resolution compared with the schlieren method, it has gained interest because of its simple apparatus.
Despite the progress of these studies, the application of the 3D-BOS method for measuring the density distributions near a wall is still challenging. This is because the wall blocks most of the light paths, resulting in a lack of light from multiple directions for obtaining the refraction data for reconstructing the density distributions. This problem may arise for near-wall flows, such as transonic buffets (Kouchi et al. 2016), shock wave-boundary layer interactions (Bolton et al. 2017), and supersonic impinging jets (Akamine et al. 2021). In the third case, a supersonic jet impinges on an inclined flat plate, forming a 3D shock structure on the wall, that interacts with a turbulent jet shear layer to produce intense acoustic waves. To further understand the noise source mechanism, it is necessary to measure the unsteady behavior of 3D shock structures and the large-scale turbulent fluctuations of shear layers. Such measurements have been difficult because, for example, all fluctuations on the wall are integrated into schlieren images and the shear layers slightly distant from the wall cannot be captured by measuring the wall pressures. These difficulties can be overcome if the 3D-BOS can measure near-wall density distributions.
For near-wall density measurements based on the 3D-BOS method, modifications are required to provide sufficient light paths near the wall. One possible approach is the use of transparent walls; however, the light direction is still restricted because of the total internal reflection (the critical angle is approximately 42° for glass with a typical refractive index of 1.5). To the best of the authors' knowledge, two extensions have been proposed (Hashimoto et al. 2017;Bühlmann 2020), although each has its limitations. The method of Hashimoto et al. is an application of near-field BOS (van Hinsberg and Rösgen 2014), in which a wall surface is used as a background of the BOS optical system. This has the advantage of enabling density measurements with only minor modifications from the conventional 3D-BOS. However, for example, Fig. 3 of van Hinsberg and Rösgen (2014) shows that the smallest detectable refraction angle of the near-field BOS rapidly deteriorates when the distance between the wall and density distribution approaches about 40 mm or less. The other method by Bühlmann uses the speckle BOS instead of conventional BOS. It uses a speckle pattern created by the reflection of the laser on the rough surface of the wall instead of the usual background pattern, enabling measurements up to the vicinity of the wall. However, Bühlmann also noted that a phenomenon called speckle decorrelation occurs, in which the speckle pattern not only shifts but also deforms when the refraction angle is large. A method for designing optical systems to avoid speckle decorrelation has not yet been established. A method that is more similar to the conventional BOS with higher sensitivity in the vicinity of the wall would be useful. This paper proposes a new extension, 3D-BOS using Mirror, to measure the density field near the wall surface to overcome the limitations. The proposed method uses the wall as a mirror to provide sufficient light paths for the reconstruction. Unlike the method of Hashimoto et al., the sensitivity in the near-wall region can be maintained by using the mirror image of the background away from the density distribution. Additionally, unlike Bühlmann's method, most of the optical settings are the same as in the conventional 3D-BOS (Nicolas et al. 2017a). A major difference from the conventional method is that the light paths are reflected by the mirror. The refraction data for the reconstruction becomes complicated because both the direct and mirrorreflected images are superposed. This paper first proposes a formulation for such complicated data. Next, the proposed method is validated using artificially generated model data of an ideal axisymmetric density distribution. The differences between the true and reconstructed distributions are compared with those of the conventional method. Finally, the proposed method was experimentally demonstrated by measuring density distributions of a candle plume.

Methods
The proposed method, 3D-BOS using Mirror, shares most of the concepts of the measurement and analysis with the conventional 3D-BOS (Nicolas et al. 2016(Nicolas et al. , 2017a. In this section, the methods are described in four parts to illustrate the similarities and differences between the conventional and proposed methods: (i) measurement setup (Sect. 2.1), (ii) measured data (Sect. 2.2), (iii) reconstruction analysis from the measured data (Sect. 2.3), and (iv) camera calibration in experiments (Sect. 2.4). Figure 1a shows a typical measurement setup for a conventional 3D-BOS. Around the target density distribution, multiple cameras are placed within 180° on one side, and background boards with patterns are placed on the other side. In contrast, the proposed method considers a case in which a flat wall is located adjacent to the density distribution, as shown in Fig. 1b. Unlike conventional methods, the cameras cannot be placed behind walls. Thus, the cameras are placed within 90° on one side and the background boards on the other side. Using the wall as a mirror to reflect light, the proposed method provides light paths from directions blocked by the wall (Cam 1 to Cam 11 in Fig. 1b). It may be desirable to provide light paths that are almost parallel to the mirror to accurately measure the density gradient normal to the mirror plane. However, depending on location of the background boards and the sizes of the density distribution and mirror, there is a minimum camera direction relative to the mirror plane at which the reflected images of the background boards can be captured. In the experiment described in Sect. 3.2, an additional camera was placed almost parallel to the mirror (Cam 12 in Fig. 1b).

Measured data
The data measured in the setups described in Sect. 2.1 is explained using the schematic of the light paths shown in Fig. 2. Despite differences owing to the mirror reflection, both the conventional and proposed methods capture images with cameras and calculate the light direction vectors before and after passing through the reconstructed region including the target density distribution. The directional changes of those vectors are used to reconstruct the density distribution in Sect. 2.3. Specifically, the light direction vectors are calculated based on the following concepts. Figure 2a shows the light paths in the conventional method. Images are captured by cameras when the density distribution is present (disturbed case) and removed (quiescent case). The shifts of the background pattern between the two images due to light refraction are measured using methods such as optical flow. Figure 2a shows the case when the pattern appearing at D 1 in the disturbed case is shifted to Q 1 in the quiescent case. Under the camera model assuming that all the projecting light paths are through a pinhole O, the straight light paths Q 1 -O and D 1 -O are calculated from the pixel coordinates of Q 1 and D 1 , respectively. In the quiescent case, by extending Q 1 -O,  Figure 2b shows the light paths in the proposed method. The light paths are reflected at the mirror located at a surface of the reconstructed region. As in the conventional method, the lines D 1 -O and Q 1 -O are first extended to calculate points D 2 and D 3 in the disturbed case and Q 2 in the quiescent case. By reflecting and extending the line Q 1 -Q 2 , the location of point Q 3 is obtained as the intersection with the background board. The intersection point D 4 with the surface of the reconstructed region can also be determined by reflecting and extending the line D 1 -D 3 based on the approximation of the conventional method that the light paths in the reconstructed region can be considered straight lines. From the directions of the straight lines D 1 -D 2 and D 4 -Q 3 , the unit light direction vectors e in and e out are obtained, respectively. These vectors represent the light directions at points D 2 and D 4 where the light path enters and exits the reconstructed region, respectively.

Fundamentals
This section describes the method of reconstructing the density distribution using the changes between the light direction vectors e in and e out obtained in Sect. 2.2. In this study, the proposed method was formulated based on a direct approach (Nicolas et al. 2016;Grauer et al. 2018) among the various conventional 3D-BOS methods. In both the conventional and proposed methods, a measurement vector b defined in Sects. 2.3.2 and 2.3.3 is first calculated from the light direction vectors e in and e out corresponding to each pixel on the image plane of each camera. Next, a reconstructed region is defined that includes the target density distribution. A vector denotes the density vector, which is reshaped from the density values on a grid discretizing the reconstructed region, and a matrix A represents the transformation from to b . The vector is calculated such that the norm of the difference between b and A , ‖b − A ‖ 2 2 , is minimized: The second term is a regularization term to keep the solution smooth. The regularization parameter is determined using the L-curve method for every data point (Nicolas et al. 2016), which was = 10 −4 in Sect. 3.1 and = 10 −5 in Sect. 3.2. This paper uses an iterative method called LSQR (Paige and Saunders 1982) to solve the minimization (1). The obtained density distribution has an arbitrary constant corresponding to the constant of integration. The constants are determined to minimize the difference from the true distribution in Sect. 3.1, and to make the densities at the outer surfaces of the reconstructed region become the atmospheric density in Sect. 3.2. The matrix A and vector b in Eq.
(1) were derived from the following two basic equations. One is the ray equation in geometric optics (Eq. 2.27 in (Träger 2012)): where s is the distance along the light path, e is the unit light direction vector, and n is the refractive index. The other equation is the Gladstone-Dale law (Merzkirch 1987): where G = 2.26 × 10 −4 m 3 /kg is the Gladstone-Dale constant for dry air. In the succeeding subsections, the formulations of the conventional and proposed methods are explained based on these equations. For simplicity, the light beams were treated as the principal rays in this study.

Conventional 3D-BOS
Although the conventional and proposed methods share the basic Eqs. (2) and (3) where n 0 is the refractive index of the surrounding atmosphere, and s D2 and s D3 are the distances measured along the light path from D 1 to D 2 and D 3 , respectively. By selecting two vectors e and e such that e , e , e in is an orthonormal basis, and transforming the coordinates of Eq. (4) to this basis by determining the inner product, Only two of three equations in Eq. (5) are independent because e out is the unit vector, i.e., e ⋅ e out 2 + e ⋅ e out 2 + e in ⋅ e out 2 = 1 . We select the first and second equations in Eq. (5). By using the scalars b and b on the left-hand side of the first and second equations in Eq. (5), the measurement vector b in Eq. (1) can be constituted as follows. By writing e = e x e y e z T , the right-hand side of the first equation in Eq. (5) becomes where the x, y, and z axes are those of the grid of the reconstructed region. By discretizing Eq. (6) and summarizing those for all the pixels of all the cameras, The vector b consists of the scalars b of all the pixels of all the cameras. The vector is the density vector as in Eq. (1), which consists of the density values on the grid points. The numbers of elements of the vectors b and are denoted by N i and N j , respectively. The N j × N j matrices D x , D y , and D z are the operators that transform the density vector into vectors consisting of the density gradients x , y , and z on the grid points. In this paper, these are defined as the coefficient matrices for the second-order central difference or, at the end points, second order onesided differences. The N i × N j matrix P is the operator that transforms the values on the grid points to the integral between D 2 and D 3 . In this paper, this is defined as the coefficient matrix for the trapezoidal integration of the values on the light path linearly interpolated from the grid points. The N i × N i matrices E x , E y , and E z are the diagonal matrices consisting of the scalars e x , e y , and e z for all the pixels of all the cameras, respectively. For the scalar b in Eq. (5), we can express the same form as in Eq. (7); thus, where O is the N i × N j zero matrix, and the N i × N i matrices E x , E y , and E z are the diagonal matrices consisting of the scalars e x , e y , and e z for all the pixels of all the cameras, respectively, where e = e x e y e z T . The 2N i -vector on the left-hand side and 2N i × N j matrix on the right-hand side in Eq. (8) are the vector b and matrix A in Eq. (1), respectively, for the conventional method.

Proposed 3D-BOS using mirror
This subsection summarizes the derivation of the measurement vector b and transformation matrix A for the proposed method. The major difference from the conventional method is an additional refraction in the light path after the mirror reflection, D 3 -D 4 , as shown in Fig. 2b. Substituting Eq. (3) into Eq. (2) and integrating between D 2 -D 3 and D 3 -D 4 , respectively, where s D2 , s D3 , and s D4 are the distances measured along the light path from D 1 to D 2 , D 3 , and D 4 , respectively. The unit vectors e D3 and e ′ D3 represent the light directions before and after the mirror reflection, and they have the following relationship (Householder transform): where e n is the unit normal vector of the mirror plane. Using e and e defined as in the conventional methods, Eqs. (9a)-(9c) can be summarized into two equations as follows. By determining the inner product of Eq. (9a) and e , By determining the inner product of Eq. (9b) and the vector e � = e − 2 e ⋅ e n e n , which is the mirror reflection of e , By determining the inner product of Eq. (9c) and e ′ , and using that e n is a unit vector, i.e., ‖e n ‖ 2 2 = e n ⋅ e n = 1, ds.
Substituting Eqs. (10c) and (10a) into Eq. (10b), After performing Eq. (10d), and similarly for e and its mirror reflection e ′ , we obtain the following equations: This is the equation for the proposed method corresponding to Eq. (5) in the conventional method.
The right-hand side of the first equation in Eq. (11) can be expressed as follows, similar to Eq. (6) in the conventional method: where e � = e � x e � y e � z T , and the x, y, and z axes are those of the grid of the reconstructed region. By discretizing Eq. (12) and summarizing those for all the pixels of all the cameras, The matrices D x , D y , D z , P, E x , E y , and E z are the same as in Eq. (7). The N i × N j matrix P ′ is the operator that transforms the values on the grid points to the integral between D 3 and D 4 . (10d) ds.

Camera calibration
Before the analysis described in Sects. 2.2 and 2.3, camera calibration was performed to obtain the locations and attitudes of the cameras, background boards, and mirror. Similar to the camera calibration for the conventional method (Le Sant et al. 2014), a marker board of known size was placed near the target density distribution and several images were captured with the cameras while changing the attitude of the marker board, as shown in Case A of Fig. 3a. The pixel coordinates of the markers were detected from the images, providing the intrinsic parameters, locations, and attitudes of the cameras relative to the marker boards (Zhang 2000), and thus the relative locations and attitudes between the cameras. Furthermore, as shown in Case B in Fig. 3b, the locations and attitudes of the background board were obtained by attaching the markers to the background board and applying the same analysis. The proposed method additionally requires the location and attitude of the mirror. In this study, the location and attitude of the mirror were obtained by fixing the marker board and capturing several images while changing the mirror attitudes (Takahashi et al. 2012(Takahashi et al. , 2016, as shown in Case C in Fig. 3c. To further reduce the error in the location and attitude of the mirror, we captured the mirror-reflected images while changing the attitude of the marker board placed near the density distribution, as shown in Case D in Fig. 3d, to obtain information about the locations of the light paths reflected by the mirror. A parameter optimization called the bundle adjustment (Hartley and Zisserman 2004) was then performed to minimize the difference between the pixel coordinates of the measured markers and those calculated with the obtained parameters for all the data in Cases A-D.
In Cases A and D, a circle grid (3 rows × 7 columns asymmetric grid with 2 mm dot diameter and 5 mm column spacing) was used because the markers were not in focus, and in Cases B and C, ChArUco markers (Garrido-Jurado et al. 2014) with a chessboard square length of 10 mm and marker square size of 5 mm were used for robust detection. Both the markers were detected using the OpenCV library implementation. In Zhang's method, an ideal simple camera was assumed to stabilize the estimation of the intrinsic parameters: the scale factors of the pixel coordinate in the horizontal and vertical directions were = , the skew factor between the two axes on the image plane was = 0 , the pixel coordinates of the pinhole u 0 , v 0 was fixed to the image center, and the lens distortion was ignored.

Model data tests
The proposed method, 3D-BOS using Mirror, was validated by examining whether the reconstructed distribution agrees with the true distribution for an artificially generated model data. First, a model data of an ideal axisymmetric distribution was generated. By setting cameras virtually, the values b and b were calculated for the model data using the right-hand sides of the first and second equations of Eq. (5) for the conventional method and Eq. (11) for the proposed method. With the calculated b and b , we reconstructed the distribution using the conventional and proposed methods, respectively. The differences between the reconstructed and true distributions were compared for the conventional and proposed methods to examine whether the differences between these two methods affected the reconstructed results. Because the number of cameras was known to affect the reconstruction error (Nicolas et al. 2016), we varied the number of cameras in the range N cam = 3 − 24.
The model data were defined as an axisymmetric distribution around y-axis in the region [−0.5, 0.5] × [−0.5, 0.5] × [−0.5, 0.5] with a radial profile: where r = √ x 2 + z 2 , and the radius and boundary thickness of the cylinder were set at r e = 0.25 and e = 0.08 , respectively. The grid resolution was heuristically determined to be 31 × 31 × 31. The isosurfaces in Fig. 4a show the 3D distribution of the model data, and the gray surface ( z = 0.5 ) indicates the mirror plane when the proposed method was examined. Figure 4b shows the distribution on the plane perpendicular to the y-axis in Fig. 4a. The distributions on the horizontal and vertical dashed lines are plotted on the upper and right sides, respectively. Figure 5a and b show examples of virtual camera configurations for the conventional and proposed methods, respectively. The red and blue lines indicate the image plane and its normal vector, and the thin gray lines are the light paths indicating the angle of view for each camera. The orange lines show the reconstructed region, and the thick gray line in Fig. 5b indicates the mirror plane. In the conventional method, N cam cameras were equally spaced on a semicircle of a radius of 10 centered at the origin. In the proposed method, N cam cameras were equally spaced within the range of 10 to 75 degrees to the z-axis, on a circle of a radius of 10 centered at the point (x, y, z) = (0, 0, 0.5) located on the mirror plane. The resolution of each camera was 51 × 29 pixels, and the scale factor to the pixel coordinates (Sect. 2.4) was set at = 300. were computed with the right-hand side of the first and second equations in Eq. (5) for the conventional method and Eq. (11) for the proposed method. It should be noted that the principal rays (i.e., infinite depth of focus) were assumed for simplicity. The Gauss-Kronrod quadrature was used to calculate each integral, and the density gradient was obtained analytically from Eq. (15) as   shows red and blue vertical lines, which were caused by refraction at the edge of the cylindrical distribution shown in Fig. 4a. Because of the axial symmetry, the same distributions were obtained for all the cameras. In contrast, for the proposed method, Fig. 5d shows a more complicated distribution, where the direct and mirror-reflected images overlapped depending on the camera directions. Figure 6 shows the distribution reconstructed by solving Eq. (1) with the proposed method with N cam = 8 cameras. The 3D isosurfaces in Fig. 6a show an axisymmetric distribution as in Fig. 4a. Figure 6b shows the slice view as in Fig. 4b. In the upper and right plots, the orange lines coincide with the black lines, indicating that the reconstructed distribution agreed with the true distribution. The near-wall cylindrical distribution could be accurately reconstructed from the complicated data in Fig. 5d using the proposed method. Figure 7 shows that the proposed method can accurately reconstruct the distribution with a different number of cameras. Figure 7a shows the relative error ‖ − true ‖ 2 ∕‖ true ‖ 2 between the true distribution true and reconstructed distribution . The horizontal axis is the number of cameras. For both the conventional and proposed methods, the error decreased to approximately 2%-3% as the number of cameras increased. These two plots almost agreed, indicating that there is no significant error inherent to the proposed method. Figure 7b shows the reconstructed distributions for N cam = 3, 4, 6, 8, and 24 cameras. The left two columns show the reconstructed distributions of the conventional method, and the right two columns show those of the proposed method. The profiles on the vertical dashed lines in the first and third columns are plotted on the second and fourth columns, where the orange and black lines represent the reconstructed and true distributions, respectively. For N cam = 8 and N cam = 24 , the reconstructed distributions agreed well with the true distributions for both the conventional and proposed methods, indicating that the distributions were accurately reconstructed. Additionally, for N cam = 3, 4, or 6 , the reconstructed distributions of the proposed method did not exhibit significant errors compared with the conventional method. The proposed method could reconstruct the distribution as accurately as the conventional method for all the number of cameras examined. The remaining errors may be attributed to the discretization and approximation common to both the conventional and proposed methods.
In addition to the basic validation of the concept of simple cylindrical distribution described above, we investigated whether the proposed 3D-BOS using a mirror could reconstruct more complicated model data. Future applications of the proposed method, such as supersonic jets impinging on an inclined flat plate, should include wall-attached asymmetric density distributions. Additional model data were reconstructed to investigate the applicability of measuring these distributions. The model data were defined using the following profile: where is the distance on a line inclined at 30° relative to the mirror plane (the origin is the center of the mirror plane), and r and are the radial distance and azimuthal angle around the -axis, respectively. The mean radius of the cylinder was r 0 = 10 mm, the thickness was e = 3 mm, the azimuthal wavenumber was n = 8 , the amplitude was Δr = 1 mm, the length of the cylinder was 0 = 20 mm, and the boundary thickness in -direction was = 3 mm. As a realistic setting, the grid and camera resolutions were the same as those employed in the experiments in Sect. 3.2. The grid discretized the 60 mm × 60 mm × 30 mm region into 81 × 81 × 41 voxels (0.75 mm-side). N cam = 12 cameras were equally spaced within the range of 10-75° of the z-axis on a circle with a radius of 500 mm from the mirror center. The resolution of each camera was 180 × 135 pixels, and the scale factor of the pixel coordinates (Sect. 2.4) was set at = 925. Figure 8a shows the 3D isosurface of the true distribution. The slice view of the center plane is shown on the left side of Fig. 8b, and the profile of the dashed line is shown on the right side. Figure 8c shows the distributions of b and b (horizontal and vertical components, respectively) obtained through numerical integration as in the case of the simple model data described above. Using these data, the distributions shown in Fig. 8e and f were reconstructed using the proposed method. Features such as overall shape and azimuthal wavenumber were qualitatively captured. Although optimization of factors such as the camera layout is required, this result indicates the applicability of the proposed method to a wall-attached asymmetric density distribution.

Experimental demonstration for a candle plume
As an experimental demonstration of the proposed method, 3D-BOS using Mirror, the density distributions of a candle plume was measured. The candle plume was used for demonstration in a previous study of the conventional 3D-BOS (Nicolas et al. 2016), because the density gradient is large, making it easy to measure, and the density distribution is simple and almost cylindrical immediately above the flame.
The experimental setup is shown in Fig. 9a. A 200 mm × 200 mm, 2 mm-thick first-surface mirror was used as the wall, and a candle was placed at a distance of approximately 20 mm from the wall. Twelve cameras (The Imaging Source, DMK33GX273) and background boards were placed around the candle. The parameters of the BOS optical systems were determined by considering the sensitivity and diameter of the circle of confusion, as in the conventional method (Nicolas et al. 2017a), because mirror reflection does not affect these factors. The distance from the camera to the candle was m ≈ 500 mm and the distance from the candle to the background was l ≈ 300 mm. A focal length was f = 25 mm and the f-value was f # ≈ 8 . The diameter of (a) Relative error (b) Slice view (colormap is the same as in Fig. 6(b); black, true distribution; orange, reconstructed distribution) Although the camera had a resolution of 1440 × 1080 pixels, the resolution was reduced to 1/8 (i.e., 180 × 135 pixels at maximum) after the optical flow analysis described below was conducted, resulting in a scale factor of the pixel coordinates (Sect. 2.4) of ≈ 925 . The exposure time was 2 ms and the frame rate was 1 Hz. A total of 231 images were captured by synchronizing all the cameras with a pulsed signal. For the background, semi-randomly distributed dots proposed by Nicolas et al. (2017a) were used (dot diameter = 0.4 mm, average dot spacing = 0.6 mm, and standard deviation of dot displacement = 0.1 mm). A K-type thermocouple was placed approximately 80 mm above the candle to measure temperature.
Before the measurement, the camera calibration described in Sect. 2.4 was performed by capturing a total of 1754 images, as shown in the bottom row of Fig. 3. Figure 9b shows the obtained locations and attitudes of the cameras, background boards, and mirror: the thick gray lines are the mirror and background boards, the thin gray lines indicate the angle of view of the cameras, and the red and blue lines represent the image plane and its normal vector of each camera. The 60 mm × 45 mm × 45 mm region immediately above the candle, indicated by the orange lines, was set as the reconstructed region. The grid resolution was heuristically determined to be 81 × 61 × 61 (0.75 mm-side voxels). Figure 9c shows an example image captured with Cam 6 in the disturbed case. The direct and mirror-reflected images of the thermocouple and flame appear on the right and left sides, respectively. Behind them, the dots drawn on the background board appeared as shown in the enlarged part (100 × 100 pixels). The shifts in these dots were measured using the Farnebäck method, a type of optical flow implemented in the OpenCV library, with an averaging window size of 15 × 15 pixels. To reduce the extra resolution, the images were resampled by averaging each 8 × 8 pixel window. The displacement due to thermal deformation of the mirror was subtracted, resulting in the distribution as shown in Fig. 9d. Only the distributions of the horizontal displacements are shown; the vertical components, which were approximately zero, are not shown. Similar to Fig. 5d of the model data, the direct and mirror-reflected images were superposed, producing a complicated distribution that differed for each camera. Excluding the shaded areas in which the flame or thermocouple prevent the dot-shift measurements, we calculated the values b and b using the left-hand side of Eq. (11). Figure 10 shows examples of the instantaneous density distributions reconstructed from b and b by solving Eq. (1) using the proposed method. Figure 11 shows the time-averaged density distribution. Both show the value − 0 , where 0 is the atmospheric density. In the instantaneous distributions, the high-temperature, low-density region owing to the candle plume appeared to be circular ( Fig. 10a and b) or deformed ( Fig. 10c and d). The time-averaged distribution in Fig. 11 captures a nearly cylindrical distribution of the candle plume. These results showed that the proposed method can experimentally reconstruct the almost cylindrical density distribution near the wall.
The white circles in Figs. 10a, c, and 11a show the location of the thermocouple determined from the images. The temperature history measured using this thermocouple is shown in the upper part of Fig. 12, and density history converted by the equation of state for the ideal gas is shown in the lower part of Fig. 12 as a blue line. The orange line in Fig. 12 shows the density history immediately below the thermocouple reconstructed with the proposed method. Both plots show a fluctuating but almost steady densities. The time-averaged density was approximately −0.81 kg/m 3 for the thermocouple and −0.88 kg/m 3 for the proposed method, with a difference of approximately 9%. This may be attributed to the composition of the candle plume. These plots were obtained using the gas constant of air, R ≈ 287.1 J/kg-K, to convert the thermocouple temperature to density, and the Gladstone-Dale constant of air, G = 2.26 × 10 −4 m 3 /kg, to calculate Eq. (11), although the candle plume contained the combustion gas. If all the oxygen in the air were used to combustion of the paraffin wax of the candle, the gas constant would be R ≈ 288.4 J/kg-K, a change of approximately 0.5%, and the Gladstone-Dale constant would be G ≈ 2.41 × 10 −4 m 3 /kg, a change of approximately 7%, which was roughly estimated based on the studies of Merzkirch (1987) and Qin et al. (2002). Considering these uncertainties, the densities obtained from the thermocouple and reconstructed with the proposed method were consistent.

Conclusions
This paper proposes a new extension, 3D-BOS using Mirror, of the conventional 3D-BOS to measure three-dimensional near-wall density distribution by using the wall as a mirror to reflect light. The modifications to the formulation of the conventional method to incorporate the mirror-reflected light paths were presented. A validation using the artificially generated model data of the ideal axisymmetric distribution showed that the proposed method could reconstruct the density distribution as accurately as the conventional 3D-BOS for all the number of cameras examined. In an experimental demonstration using a candle plume, the proposed method reconstructed the almost cylindrical density distribution near the wall. Future work will include improving the results using newer optical flow algorithms, a more realistic treatment of light rays as conical beams, and more parametric studies on grid resolutions and camera configurations. It is also expected that the proposed method can be applied to unknown density distributions such as supersonic jets impinging on an inclined flat plate.