1. Introduction

With the development of aerospace and industrial design, measurement of the flow field has become increasingly important. Thus, various measurement methods have been developed in recent decades. One of these is the optical-based measurement method, which can capture the instantaneous flow field properties of the entire field. It has the advantage of being a noncontact and nonintrusive technique. Therefore, it has been widely used for the measurement of flow fields. Interferometry [1–3], shadowgraph [4,5], and schlieren [6–8] methods are among the most common approaches. Interferometry can provide quantitative results; however, it has an extremely complex setup and requires strict conditions in the surrounding environment. Although the shadowgraph method has the simplest experimental setup, it is only suitable for flow fields with large gradient variations and only provides qualitative analyses [9].

The traditional schlieren setup generally consists of a light source, lens, camera, and sharp knife edge. Since its development around the 17^th century [10], the schlieren method has been adopted for the quantitative measurement of the refractive index, density, and temperature of the flow field [7,11,12]. Furthermore, schlieren tomography [13–15], in which schlieren images are obtained from multiple directions, has been developed. In 2000, the background-oriented schlieren (BOS) technique was proposed by Dalziel et al. [16] and Raffel et al. [17] at almost the same time. Compared with the traditional schlieren method, BOS has modest technical requirements, which only consist of a light source, background, and camera. The BOS technique has the advantage of a very large field of view without large, expensive, and precision optics [18]. Therefore, although BOS has some disadvantages such as vibration sensitivity and limited real-time visualization [18], it still can provide an excellent strategy for full-field measurements. And BOS can provide unique solutions to some problems [18].

The BOST method utilizes the deflection of rays passing through an inhomogeneous field for measurement. Therefore, BOST differs from general CT reconstruction. Its direct reconstruction result is the distribution of the refractive index gradient. An additional step is required to obtain the distribution of the refractive index — for example, direct line integration or solving Poisson equation.

Since the BOS method was proposed, many 2D experiments have demonstrated its applicability in the measurement of the flow field. Raffel et al. [17] first used the BOS method to measure the vortex of helicopter blades. In 2001, Richard and Raffel [19] measured a compressible vortex in the wake of a cylinder. They compared the results of the BOS technique and particle image velocimetry, and similar results were obtained. The separation of the boundary layer and shear layer, as well as vortex shedding, was clearly identified in both results. Direct line integration was used to calculate the density distribution in both studies. Richard and Raffel et al. [20] applied the BOS method to a supersonic jet. These experiments prove that the BOS method can be used for flow field measurements. However, these are 2D BOS measurements and cannot be used to obtain the 3D structure of a complex flow field.

If the flow is axisymmetric, the 3D flow can be reconstructed by inverse Abel transformation from one viewing direction, which is perpendicular to the axis of symmetry of the flow. For example, Sourgen et al. [21] compared numerical simulations and BOS results using the inverse Abel transformation. In 2014, van Hinsberg and Rösgen [9] measured an under-expanded supersonic free jet by using the inverse Abel transformation. However, the Abel transformation requires the assumption of axisymmetric flow, and it is susceptible to noise. Therefore, the filtered backprojection (FBP) algorithm was employed for the reconstruction of the BOS. In 2004, Venkatakrishnan and Meier [22] measured the density field for an axisymmetric supersonic flow over a cone cylinder model. In 2005, Venkatakrishnan [23] measured an axisymmetric under-expanded jet. The density field was obtained by solving the Poisson equation and using the FBP algorithm [22,23]. In addition, the FBP algorithm can be more universally used for nonsymmetric flows.

For a non-axisymmetric flow, BOS projections should be obtained from multiple directions for BOST reconstruction. When a stationary flow is considered, the BOS projections can be obtained by rotating the position of the camera and background relative to the flow [24], in which the multidirectional projections are obtained from 19 projection angles from 0° to 90° at 5° intervals. When the projection angle is limited or there are few projection data, the FBP algorithm is prone to artifacts [25]. Experiments performed by Ota et al. [26–28] used the same rotating strategy as the previously mentioned study [24], while the algebraic reconstruction technique (ART) was employed for the reconstruction. ART is more practical for reconstruction with limited data than the FBP algorithm. The density field was reconstructed from 30 projections by Ota et al. [29]. However, for the measurement of complex flow fields, this rotation strategy is only an expedient.

When a complex flow is considered, a synchronous multicamera setup is employed to record the BOS projections from different directions. The first multicamera BOS measurement of nonstationary flow was demonstrated by Atcheson et al. [30] in 2008. There were 16 cameras surrounding the hot air flow above a gas burner in a 180° arc, and the refractive index field was calculated by Poisson integration. Zeb et al. [31] employed three cameras to measure the convectional flow field. However, the reconstruction results had strong artifacts because of the small number of cameras. In 2015, Nicolas et al. [32] designed a 3D BOS experimental facility with 12 noncoplanar cameras and measured three different hot air flows. A multicamera BOS approach was also used in the measurement of compressible flow [33] and a supersonic wind tunnel [34]. Since 2018, 3D BOS has also been used in the measurement of the flame field [35–38].

The deviation angle of the BOS is proportional to the integral of the refractive index gradient; thus, the direct reconstruction result obtains the distribution of the refractive index gradient. Then, an integration step is performed to calculate the refractive index field. The first method is direct line integration [17,19,25] which is a simple method, while this method produces line noise [19]. It is easy to produce noise accumulation during integration. The second extensively used method is to solve the Poisson equation [17,20,22,23,27,30,39]. However, this method is complicated and has the disadvantage of strongly smoothing the results [17]. The key to obtaining the refractive index distribution directly is to address the derivation process in the discretization. In 2015, Nicolas et al. [32] proposed a direct reconstruction method. To eliminate the integration step, they combined a finite-difference matrix (FDM) in the discretization process. Recently, this method has been adopted in many studies [35,37,38,40]. It is an excellent method; however, there are many kinds of differences, such as forward, backward, central, and Richardson differences. The order of the differences also influences the accuracy.

In this article, a calculation method of projection matrix is proposed based on the RBF for BOST. Unlike the “mechanistic” approach of Nicholas et al. [32], the projection matrix is calculated using the differentiable property of RBF. Thus, the weights in the projection matrix are the integration of partial derivatives of basis functions. The integration step is eliminated, and no additional finite-difference matrix is used. By using the 3D RBF, the derivation process can be combined with the calculation of the projection matrix. In Section 2, the BOS theory is described. In Section 3, the BOS projection model based on the 3D RBF is presented. Section 4 presents the results of the numerical simulation and the experimental reconstruction. Finally, conclusions are given in Section 5.

2. BOS theory

The BOS theory can be found in many articles, but, for the convenience of readers, it is briefly introduced here. The flow field is reconstructed by the BOS method through the deviation of the ray when propagating in an inhomogeneous field. A schematic diagram of the BOS theory is illustrated in Fig. 1. The background pattern is distorted when there is an inhomogeneous field between the background and the camera. The yellow and green lines are the ray paths with and without flow, respectively. In addition, $\Delta u$ and $\Delta v$ are the displacements of the ray in the imaging plane when there is a flow, ${l_A}$, ${l_B}$, and ${l_C}$ are the distance from the background to the probe volume, distance from the probe volume to the lens, and distance from the lens to the imaging plane (i.e., the focal length $f$), respectively. The deviation angle $\varepsilon$ is the integration of the refractive index gradient along the ray path:

Here, ${\varepsilon ^{(\alpha )}}$ is the deviation angle in the $\alpha$ direction. ${n_0}$ is the environmental refractive index and $ray$ is the ray path. Under the assumption of paraxial recording and a small deviation angle, the relationship between the deviation angle and displacement in the imaging plane is given by [41]

 

Fig. 1. Schematic diagram of the BOS theory: the yellow and green lines are the ray path with and without the field, and the ray path indicated by the green dashed line is actually nonexistent.

Download Full Size | PDF

Because the experimental setup is relatively simple, the BOS setup can be conveniently arranged as a multidirectional measurement device. As shown in Fig. 2(a), multiple synchronized cameras pass through the probe area to capture images of the background. The displacements can be extracted by comparing the distorted and undistorted images of the background. Thus, the refractive index can be reconstructed using displacements. In some cases, it may be necessary to obtain other parameters of the flow field, such as the density and temperature. The refractive index is related to the density using the Gladstone–Dale equation [20].

Here, n is the refractive index, and $\rho$ is the density. Furthermore, ${K_{GD}}$ is the Gladstone–Dale constant, which is a function of the light wavelength. Combined with the ideal gas law, the temperature T can be obtained from the density, as follows:

 

Fig. 2. (a). Schematic diagram of multidirectional BOST. (b). Probe area discretization and ray path: the orange grid box is an RBF area in the green grid with a coverage range of $3d$, the blue line is the ray that can be used in the reconstruction, and the red line is the ray that cannot be used in the reconstruction.

Download Full Size | PDF

Here, P represents the atmospheric pressure, M represents the molar mass of the gas, and R represents the universal gas constant. Therefore, the temperature can be calculated using the refractive index.

3. BOST

3.1 BOS projection model based on RBF

To obtain the distribution of the refractive index directly, a BOS projection model based on the 3D RBF is proposed. The first basis-function-based CT was presented by Hanson and Wecksung [42]. This method is mainly used in 2D image reconstruction [42–46], but its application in 3D field reconstruction has rarely been reported. By combining the basis function with the BOS projection model, the distribution of the refractive index can be obtained directly instead of the distribution of the refractive index gradient.

As shown in Fig. 2(b), the probe area is discretized into grids of $L \times W \times H$, where each grid is a cube with side length d. For each grid, there is a 3D RBF with a certain coverage range, and the centroid of the RBF is at the center of the grid. An RBF area with a coverage range of $3d$ is shown in Fig. 2(b). In this approach, a function is represented as a linear expansion of the basis functions. It is assumed that the refractive index distribution in the 3D space is $n({\textbf x})$, and the basis function is $g({\textbf x})$. Therefore, $n({\textbf x})$ can be approximated with a set of N basis functions,

Here, $N = L \times W \times H$ is the number of basis functions, and, ${\textbf x} = {[x,y,z]^\textrm{T}}$ is a point in the probe area. ${g_i}$ is the ith basis function, and ${c_i}$ is the coefficient of ${g_i}$. Note that typically the coverage of the basis function is larger than one grid, therefore, the basis function will overlap with the surroundings. This work employs Gaussian RBFs which are defined as follows:

Here, ${{\textbf x}_i} = {[{x_i},{y_i},{z_i}]^\textrm{T}}$ is the centroid of ${g_i}$, $\sigma$ is the standard deviation of ${g_i}$, ${d_{th}}$ is the threshold distance (e.g., ${d_{th}} = 3\sigma$), and ${||\cdot ||_2}$ is the Euclidean norm. The α-direction gradient of $n({\textbf x})$ for $\alpha \in \{ x,y,z\}$, is

One can express the standard BOST measurement model using this basis. To begin, the α-direction deflection of the jth ray is given by a line-of-sight (LoS) integral over ${\nabla _\alpha }n({\textbf x})$ alone that ray,

All the coefficient of basis function and deflection data are collected into the vector, ${\textbf c} = \{ {c_i}\} _{i = 1}^N$ and ${{\mathbf \varepsilon }^{(\alpha )}} = \{ \varepsilon _j^{(\alpha )}\} _{j = 1}^K$, respectively. And ${{\textbf T}^{(\alpha )}}$ (defined in Eq. (10)) is formed for $i \in 1,2,\ldots ,N$, $j \in 1,2,\ldots ,K$, and $\alpha \in \{ x,y,z\}$. The result is a $3K \times N$ linear equation:

3.2 Calculation of the BOS projection matrix

To calculate the BOS projection matrix ${\textbf T}$, one must first know the ray trajectory equation corresponding to the sampling point on the image. Therefore, the pinhole model is employed to describe the camera imaging process. In the pinhole model, a 3D point ${(x,y,z,1)^T}$ in the world coordinate system is imaged as ${(u,v,1)^T}$ in the imaging plane by

Here, ${\textbf P}$ is the camera projection matrix, ${\textbf K}$ is the intrinsic matrix of the camera, ${\textbf R}$ is a rotation matrix, ${\textbf t}$ is the translation vector, and $\varsigma$ is a nonzero scale factor. After one obtains the camera projection matrix by camera calibration, the ray trajectory equation of the sampling point ${(u,v)^T}$ can be calculated by

Therefore, the points at which the ray enters and exits the probe area can be calculated. The grids that are crossed by the ray can also be calculated sequentially. The grids crossed by the rays are illustrated in Fig. 3. Then, the weights of these grids can be calculated one by one using Eq. (10). However, a ray only passes through a small part of the grids along the ray path; thus, most of the elements in the matrix ${\textbf T}$ are zero. Therefore, the matrix ${\textbf T}$ is stored as a sparse matrix. The programs are written in the c++ programming language. And each non-zero element in ${\textbf T}$ is represented by a triple, in which two variables of int type represent the position of the element and a variable of float type represents the weight. Therefore, the memory consumption can be significantly reduced.

 

Fig. 3. The grids crossed by the rays. For the convenience of viewing, draw a grid every 4 grids, and the grid is enlarged 3 times.

Download Full Size | PDF

The ray is discretized in the calculation. The discrete step can be set to $d/10$ [47]. A 3D RBF overlaps with the surrounding RBF. The overlap range depends on the coverage of the RBF. Therefore, the value of the RBF at a certain point should be added to the value of the surrounding RBF. To facilitate the calculation, the coverage of the RBF is generally set to an odd number. The rays passing through the top or bottom of the probe area should not be used in the reconstruction. As shown in Fig. 2(b), because there may still be a flow field above or below the probe area, the path used to calculate the integral is incomplete. Thus, rays that pass through the top or bottom of the probe area should be excluded.

The displacement in the image is a 2D vector $(\Delta u,\Delta v)$, whereas the deviation angle used in Eq. (11) is a 3D vector $({\varepsilon _x},{\varepsilon _y},{\varepsilon _z})$ — that is, the image displacement is decomposed into three directions of the world coordinate system, as given by Eq. (14).

Here, ${r_{ij}}$ is the element in the rotation matrix of the camera obtained by camera calibration.

3.3 Reconstruction

The Eq. (11) is commonly ill-posed in BOS theory. The projection matrix ${\textbf T}$ is typically large in BOS theory; thus, the matrix-theory-based method is computationally unpractical. Iterative reconstruction algorithms are widely used for such problems. Liu et al. [37] showed that ART outperforms the simultaneous iterative reconstruction technique and the Landweber algorithm in BOS reconstruction. Therefore, the ART algorithm was used in the experiments. The iterative process of ART is expressed as

Here, ${{\textbf c}^{(i,k)}}$ is the result updated by the i-th equation after the k-th iteration, and the initial guess can be filled with zero. In addition, ${\varepsilon _i}$ is the i-th component of the vector ${\mathbf \varepsilon }$, ${{\textbf T}_i}$ represents the i-th row of the projection matrix ${\textbf T}$, and $\lambda$ is the relaxation parameter.

Using ART, it is relatively easy to incorporate prior knowledge into the reconstruction process. Total variation (TV) regulation [48] was originally presented for image denoising, which can suppress noise while retaining the edge and sharpness of the image [49]. The TV regularization was used in BOST for the first time by Grauer et al [35], and achieved excellent results. There are often a lot of noise in the results of ART; therefore, TV regulation is introduced into the ART reconstruction process. The TV value is expressed as

Here, $TV ({\textbf c})$ is the TV of ${\textbf c}$ in the 3D discrete form. The minimization of $TV ({\textbf c})$ can be written as

The result of ART is updated by TV regularization based on the split Bregman method [50] at each iteration, and the updated result is used as the initial value of the next ART iteration.

The overall process of the 3D BOST is shown in Fig. 4. The detailed steps are listed below.

 

Fig. 4. Flow chart of 3D BOST: the blue rectangles represent some kind of computation, and the others represent the data or result.

Download Full Size | PDF

4. Experiments and results

4.1 Numerical simulation

First, the feasibility of this method was validated by simulation data. Twelve simulated cameras were evenly distributed within a range of $165^\circ$. The distances from the camera to the probe area and background plate are 600 and 1300 mm, respectively. The temperature field used is flame-like:

There are four small flames surrounding a large flame evenly at a distance of 25 mm. The diameter and height of large flame are 20 mm and 80 mm, respectively. The diameter and height of the small flame are 10 mm and 40 mm, respectively. The max value of temperature is ∼1500 K. The forward projection process can be simulated according to Eqs. (1) and (2); thus, one can calculate the displacements in the imaging plane of each camera under the influence of the temperature field by ray tracing. For example, the displacements in the fourth camera are shown in Fig. 5(a).

 

Fig. 5. Simulated displacements in fourth camera: (a) displacements with 0% Gaussian noise, (b) displacements with 10% Gaussian noise, and (c) displacements with 20% Gaussian noise.

Download Full Size | PDF

In real experiments, the errors caused by the vibration, light fluctuation, displacement extraction, and image noise may influence the reconstruction. Therefore, to evaluate the influence of noise, we added 0%–20% Gaussian noise to the displacements. The mathematical expression for adding noise is as follows:

Here, $GN$ is the Gaussian noise distributed between 0 and 1, $NL$ is the noise level, ranging from 0% to 20%, and $MD$ is the maximum displacement of the image. The displacements with different Gaussian noise are shown in Fig. 5. The size of the probe area is $100 \times 100 \times 100$ mm and is divided into $100 \times 100 \times 100$ grids. In each imaging plane, $167 \times 167$ points were selected. Therefore, the size of the matrix ${\textbf T}$ is $1004004 \times 1000000$ and is much larger in the following real experiments. Therefore, the matrix ${\textbf T}$ is stored as a sparse matrix to reduce the memory consumption.

To investigate the effect of RBF coverage on reconstruction results, the reconstruction is performed under different RBF coverages. The coverage of the RBF is determined by the $\sigma$ in Eq. (6). Moreover, the mean relative error ($MRE$) can be calculated. The MRE is given by

 

Fig. 6. Slices of reconstruction results. (a). Reconstruction result of the method of Nicolas et al. (FDM). The reconstruction at the edges is unsatisfactory (i.e., left of the figure), therefore, the error is relatively large. (b)∼(f). Reconstruction results of our method under different RBF coverages.

Download Full Size | PDF

Slices of the reconstruction results are shown in Fig. 7. The figure shows that the BOS projection model can reconstruct the temperature field accurately. However, owing to the influence of Gaussian noise, there are also many noises in the reconstruction results of the ART algorithm. Moreover, in the BOS experiment, the displacement is relatively small; thus, the reconstruction is easily influenced by noise. TV regularization can effectively suppress these noises and make the results smoother. Figure 8 shows the profiles of the reconstruction results. The results of the TV regularization are smoother and have less noise. Especially at the edge of the reconstruction area, noise can be effectively reduced. The figure shows that the TV regularization can constrain the results effectively at different noise level. As shown in Fig. 9, the $MRE$ of the different noise levels are calculated from 0% to 20%. The $MRE$ under TV regularization is significantly reduced in both methods, however the $MRE$ of our method increases more slowly as the noise level increases.

 

Fig. 7. Slices of reconstruction result of our method in the plane $z = 0$: the first row is the result of ART, and the second row is the result of ART with TV regularization — each column is the reconstruction result at a different noise level (NL).

Download Full Size | PDF

 

Fig. 8. Profiles of slices presented in Fig. 7 along X = 0.

Download Full Size | PDF

 

Fig. 9. $MRE$ under different noise levels. The x-axis is the noise level, and the y-axis is the $MRE$. The $MRE$ under TV regularization is significantly reduced.

Download Full Size | PDF

4.2 Experimental setup and calibration

The experimental setup consisted of 15 AVT Guppy F-125B cameras mounted on an ${\sim} 160^\circ$ arc surrounding the hot air above the candle flame, as shown in Fig. 10. A Sony ICX445 1292- × 964-pixel sensor was included in each camera and attached to a lens with a focal length of 12 mm. All cameras were synchronized by an external trigger signal, and the frequency of the trigger signal was 30 Hz. Three background plates painted with a sine wave pattern were fixed in a semicircle approximately 750 mm away from the candle, and the distance from the candle to the camera was approximately 650 mm. The background plates were illuminated using six light-emitting diode light sources.

 

Fig. 10. Experimental setup. Fifteen cameras are arranged around the hot air flow above the candle flame. The cameras numbered with green are for reconstruction, and the cameras numbered with blue are for validation. The background plates were illuminated using light-emitting diode light sources.

Download Full Size | PDF

The main problem of this multicamera calibration is that it is impossible for all the cameras to view the calibration plate. Therefore, the multicamera calibration procedure was performed using the flexible multicamera calibration method with a rotating calibration plate [52]. The 15 cameras were divided into three groups, with five cameras in each group. A chessboard calibration plate was placed in the probe area, as shown in Fig. 11(a). When the calibration plate rotated, all cameras simultaneously captured the images of the calibration plate. Within each group, the 3D world coordinates of the rotating calibration points were calculated by multiview stereo vision. Then, by using the overlapping calibration points of the adjacent group, the three groups of world coordinate systems could be transferred to a unified world coordinate system. Finally, the positions and orientations of the cameras were calculated based on the pinhole model.

 

Fig. 11. (a) Calibration plate and rotation device and (b) 3D coordinates of rotating points of three groups in a unified world coordinate system.

Download Full Size | PDF

The 3D coordinates of the rotating points of the three groups are shown in Fig. 11(b). The calibration results are shown in Fig. 12. The reprojection errors of all cameras were within 0.3 pixels. The absolute mean error (AME) and standard deviation (Std) of the reprojection error along the u-axis were 0.02 and 0.0257 pixels, respectively, and the AME and Std of the reprojection error along the v-axis were 0.025 and 0.033 pixels, respectively. After the intrinsic and extrinsic parameters of the cameras have been obtained, the ray path of any point in the image can be calculated. The four calculated ray paths are shown in Fig. 12(a).

 

Fig. 12. Calibration results: (a) positions and orientations of 15 cameras (the red, green, and blue lines in the middle represent the world coordinate system, and the yellow lines are the ray path of four sample points) and (b) reprojection error of 15 cameras.

Download Full Size | PDF

4.3 Experimental results

By comparing the distorted and undistorted images, the displacements in the imaging plane can be calculated using the cross-correlation algorithm [51]. As shown in Fig. 13, $126 \times 185$ sample points were selected in each image. The size of the probe area was $120 \times 120 \times 192$ mm, and it was divided into $200 \times 200 \times 320$ grids.

 

Fig. 13. Displacements of all cameras: the cyan cuboid is the projection of the probe area on the image, and the orange, green, and blue lines in the middle represent the projection of the world coordinate system.

Download Full Size | PDF

The projection matrix was calculated according to the description given in Section 3 and Nicolas et al.’s method (FDM) [32]. The reconstruction process was performed using ART and TV regularization. To validate the correctness of the method, 12 cameras (green in Fig. 10) were selected to participate in BOS reconstruction, and the other three cameras (blue in Fig. 10) were used for validation. After the reconstruction result was obtained, the reprojection displacements of the other three cameras could be calculated according to Eqs. (1) and (2) by using ray tracing. The comparisons between the measurement displacements calculated by the cross-correlation algorithm and the reprojection displacements calculated by ray tracing are shown in Fig. 14(a). The correlation coefficient (CC) was used to indicate the similarity of each pair of images. The CC of images A and B is given by

Here, ${\mathop{\rm cov}} (A,B)$ is the covariance of images A and B, and ${\mathop{\rm var}} (A)$ is the variance of image A. The CCs of each pair of images were calculated to be 0.94, 0.91, and 0.90, which is larger than the result of FDM. The figure shows that the measurement displacements and reprojection displacements are similar to each other. Further, we conducted a more detailed analysis of measurement displacements and reprojection displacements. We use Euclidean distance to evaluate the error between two displacements. As shown in Fig. 14(b), the position with a larger error is clearly visible in the figure. The error of FDM is larger in the lower part of the plume. Compared to FDM, our method produces smaller errors. Besides, the distribution of the reprojection errors is shown in Fig. 16. The AME and Std of the reprojection error along the u-axis are 0.049 and 0.080 pixels, and the AME and Std of the reprojection error along the v-axis are 0.027 and 0.041 pixels, respectively. Only a few errors are relatively large, and 90% of errors are less than 0.05 pixels. The comparisons between our method and Nicolas et al.’s method in different frames are listed in Table 1. The slices of results are shown in Fig. 15. The figure reveals that a complicated structure can also be reconstructed. However, we notice that there are the petal-like artifacts in the surrounding area, which is similar to the results in Ref. [32]. In this article, though the “3D mask” can improve the accuracy, the artifacts still exist. The number of petal-like artifacts are related to the number of projection angle in the linear tomography [53]. In the results of FDM, the refractive index is larger than ours in the lower part of the flow. This problem can be improved by using a “3D mask” [32]. However, the “3D mask” is not employed in our method. Besides, the gradient of refractive index is shown in Fig. 17. In this experiment, 50 frames of BOS images were captured at a frequency of 30 Hz, and BOS reconstruction was performed using the proposed method. The reconstruction results for five different frames are shown in Fig. 18. The quantitative reprojection errors were calculated as described above. The details are listed in Table 1. The reprojection errors of different frames are approximately the same; therefore, the stability of the method is verified.

 

Fig. 14. (a). Comparisons between measurement displacements and reprojection displacements: the first row is the measurement displacements. The second and third rows are the reprojection displacements of our method and Nicolas et al.’s method (FDM), respectively. (the rendering scalar is consistent with Fig. 13). (b). The Euclidean distance error between measurement displacements and reprojection displacements. The first row is the error of our method, and the second row is the error of Nicolas et al.’s method (FDM).

Download Full Size | PDF

 

Fig. 15. Slices of one instantaneous reconstruction result: slices along the vertical direction from −17.7 mm to −2.7 mm, and slices along the horizontal direction from −47.7 mm to 6.3 mm. The first row is the result of our method, and the second row is the result of Nicolas et al.’s method (FDM).

Download Full Size | PDF

 

Fig. 16. Distribution of errors between measurement displacements and reprojection displacements: (a) scatter plot of errors, (b) frequency of different errors along the u-axis direction, and (c) frequency of different errors along the v-axis direction.

Download Full Size | PDF

 

Fig. 17. Gradient of refractive index corresponding to Fig. 15. The figure shows the gradients in three directions for the results of our method and that of Nicolas et al. (FDM).

Download Full Size | PDF

 

Fig. 18. Three-dimensional refractive index distribution of different frames. The more distributions of consecutive frames are shown in Visualization 1. And the display of distribution under different viewing directions are shown in Visualization 2.

Download Full Size | PDF

Table 1. Reprojection errors of five frames. The value in “()” is the results of Nicolas et al.’s method.

View Table

5. Conclusion

A discrete projection model based on a 3D RBF was proposed for BOST. Using this method, the 3D flow field distribution can be obtained directly without an additional integration step or additional finite-difference matrix. Firstly, the accuracy of this method was verified by simulation. The RBF coverage of $7d$ or $9d$ is recommended in the BOST by the simulation. Besides, the refractive index of the hot air flow above a candle flame was reconstructed using 12 cameras. The reprojection displacements in the other three cameras were calculated by ray tracing. Comparing the measurement displacements calculated by the cross-correlation algorithm demonstrated that the reprojection displacements were similar to the measurement displacements. This method can obtain accurate results and has excellent anti-noise ability. It can provide a feasible method for flow field measurements. However, there are many rays that need to be processed; thus, the calculation of the projection matrix is time consuming, requiring approximately 23 h in the experiments, although 15 CPU threads were used. Therefore, in future research, simplification of the BOS projection matrix to reduce the time consumption and memory usage is planned. In addition, in some experimental scenarios, there may be vibrations, water mist or light source fluctuation that influence the experiment. This will have a great impact on the extraction of the BOS image displacements. Displacements may be generated where there should not be displacements. In future research, we will focus on these problems.