# Hybrid Adaptive Wavelet-Based Optical Flow Algorithm for Background Oriented Schlieren (BOS) Experiments

Authors: Xin-yu Zhang; Li-Ming Wang; Bin Liu; Yue Luo; Xing-cheng Han
Journal: Mathematical Problems in Engineering
 (2020)

---

## Abstract

Quantitative analysis of the flow field is an effective method to study hydrodynamics. As a flow field measurement technology, the Background Oriented Schlieren (BOS) is widely used. However, it is difficult to measure the complex transparent flow field (flow field with large refractive index gradient) using the BOS experiment. In order to overcome this difficulty and improve the accuracy of the BOS experiment, this paper presents a hybrid adaptive wavelet-based optical flow algorithm for the BOS. The current algorithm is a combination of the traditional optical flow algorithm and the wavelet-based optical flow algorithm. By adding the initial value constraints, the adaptive scale constraints, and the adaptive regularization constraints, the algorithm can effectively overcome the above-mentioned difficulty and also improve its accuracy. To further illustrate the feasibility of the proposed method, this paper uses the simulation data, the data of the DNS datasets, and the data of the BOS experiment to verify the performance of the algorithm. The experiment of comparing the proposed algorithm with the existing ones is carried out on the DNS datasets and the data of the BOS experiment. Finally, the proposed method is verified by a practical BOS experiment. The results show that the proposed algorithm can effectively improve the measurement accuracy of displacements.


---

## Body

## 1. Introduction

Background Oriented Schlieren (BOS) is widely used to measure the transparent flow field owing to the simplicity of its experimental device and its high measurement accuracy [1, 2]. BOS is an optical flow diagnostic technique used to measure the refractive gradients in a flow field by tracking the distortion of a target dot pattern. Due to the uneven refractive index distribution in the transparent flow field, the light deflects after passing through the flow field; therefore, an object viewed through the transparent flow field will appear distorted due to the refraction of light rays traversing the flow field [3]. BOS relies on this principle to measure the transparent flow field.The accuracy of the displacement estimation algorithm directly determines the measurement accuracy of the BOS. The displacement field obtained by the traditional cross-correlation algorithm is not dense, and especially in the region of large refractive index gradient, the algorithm is not accurate. The traditional optical flow algorithm and the improved multiscale optical flow algorithm overcome the shortcomings of the cross-correlation algorithm, especially in the region of large refractive index gradient. However, the optical flow algorithm is limited in practical application, as it is based on the assumption of small displacement and constant brightness. In addition, compared with the cross-correlation algorithm, the optical flow algorithm is more complex and computationally expensive. Jagadesh-Ramaiah et al. propose the use of GPU acceleration to improve the computational efficiency of the algorithm [4]. Because the wavelet transform itself has the characteristics of data compression and multiscale analysis, the wavelet-based optical flow algorithm has attracted people’s attention [5–8]. This paper presents a hybrid adaptive wavelet-based optical flow algorithm for the Background Oriented Schlieren (BOS) based on previous studies [2, 5–7, 9, 10]. The wavelet-based optical flow algorithm can effectively reduce the calculation cost of the algorithm by using the data compression characteristics of the wavelet. At the same time, the algorithm can effectively deal with the problem of large-scale displacement in the BOS experimental data by using the multiscale analysis characteristics of the wavelet. The rest of this paper is organized as follows: Section 2 introduces the related works at home and abroad; Section 3 describes the proposed algorithm; Section 4 presents the simulation method of the BOS experimental device based on COMSOL and the description of the truncation error; Section 5 presents the experimental results and discussion; and Section 6 summarizes the research.


## 2. Related Works

The performance of the data processing algorithm directly affects the measurement accuracy of the BOS. In order to improve the measurement accuracy of the BOS, more and more researchers are working on optimizing the data processing algorithm. Initially, an algorithm based on cross-correlation was used to process the BOS data. However, this algorithm has some shortcomings, which can be tackled by PIVlab, a useful MATLAB toolbox, proposed by Thielicke et al. [11]. The MATLAB toolbox uses cross-correlation algorithm to calculate pixel displacement. However, the calculated displacement is minimal and imprecise. Especially in the region of high refractive index gradient, the calculation error of the algorithm is large. Therefore, the optical flow method is introduced into BOS data processing. The optical flow method does not rely on the size of the interrogation window. It estimates the displacement of each pixel according to the overall correlation between two adjacent images, rather than depending on the relationship between the pixels of a certain region or a certain pixel. Optical flow algorithm is based on global correlation estimation of the displacement field, so that it can obtain dense, high-resolution displacement fields [5], such as Horn-Schunck algorithm [12] and Lucas-Kanade algorithm [13]. However, the optical flow algorithm also has its own potential drawbacks. It is based on the assumption of small displacement. This also limits the use of the optical flow algorithm. In order to overcome the limitation of displacement, the theory of multiresolution analysis is introduced into the optical flow algorithm [9, 14]. Although the improved optical flow algorithm can overcome the limitation of displacement, the optical flow algorithm cannot break through the assumption of constant brightness. This makes it difficult to apply the optical flow algorithm to the actual experiment. It may exhibit difficulties with experimental noise or events, such as particles, entering or leaving the image plane that commonly occur in experiments [6]. In order to enhance the adaptability of the optical flow algorithm, it needs some kind of regularization constraints, which are adjusted based on experience in most cases, or according to the characteristics of the flow field. There are several optical flow methods in the literature that have been designed and optimized for divergence-free motions [15]. However, some optical flows use periodic boundary conditions that are not well suited for experimental data [16]. The multiscales optical flow methods based on coarse-to-fine blurred images sequentially also have drawbacks, such as large errors occurring in coarser estimates cannot be subsequently corrected in the finer estimates [17]. However, the wavelet-based optical flow method could cope with this difficulty. The wavelet-based optical flow method representation hierarchy is different from the typical image coarse-to-fine hierarchy. The hierarchy is built along motion resolution level, not image resolution level. The wavelet-based optical flow algorithm was first proposed by Wu et al. [18]. Based on this, Kadri-Harouna et al. [10] improved the wavelet-based optical flow algorithm and mainly explored the application of regularization constraints in the optical flow algorithm based on the wavelet. Schmidt et al. [6] improved the wavelet transform on the basis of the Dérian algorithm. The transform in the WOF-X is performed using biorthogonal 9–7 wavelets and a lifting transform, unlike in previous wavelet-based optical flow methods. However, when dealing with the complex flow field (the flow field with great change of refractive index gradient), the calculation error of the above-mentioned wavelet-based optical flow algorithm is still large [6].


## 3. Hybrid Adaptive Wavelet-Based Optical Flow Algorithm

### 3.1. Multiscale Optical Flow Representation Based on the Wavelet

I0 and I1 are images obtained by imaging the same background pattern under the same BOS experimental device; I0 is the imaging result without the flow field, and I1 is the image with the flow field. The optical flow algorithm is used to estimate the pixel displacement ux between the two images. Based on the conservation assumption of the image brightness, images I0 and I1 fulfill the following formula:(1)I1x+ux−I0x=0.In order to calculate the displacement fieldux, it is necessary to change formula (1) to the standard quadratic function, as shown in the following formula:(2)u∗=argminFd=argmin∫x∈ΩI1x+ux−I0x2dx.Here,Ω represents the set of effective pixels of image I0x.By calculating the minimum value of formula (2), the pixel displacement ux is obtained. For a point x in two-dimensional space, formula (2) is a quadratic equation of two-dimensional vector ux. However, the number of unknowns in the equation is more than that in the equation, which is an underdetermined equation. In order to better solve the above optimization problems, a convex regularization term FRu needs to be added [19]. The objective function is given by the following formula:(3)u∗=argminFu=argminFDu+λFRu.The wavelet framework introduces the decomposition of a signal into several sets of details at various scales and a remaining coarse approximation. The projection of the two components,ux, uy, of the velocity field u onto such multiresolution approximation spaces proves to be efficient within the context of optical flow estimation, enabling the handling of large displacements. Here, let’s recall the basic principle of the two-dimensional wavelet transform. For any function gx,y∈L2R2, it can be expressed linearly by Vj2j∈z,Wj2j∈z, as shown in the following formula:(4)gx,y=∑k,m=02C−1g,φC,k,mφC,k,mx,y+∑t=13∑j=CF−1∑k,m=02j−1g,ϕj,k,mtϕj,k,mtx,y.The closed subspace sequence Vj2j∈z generated by the expansion and translation of the scale function φx,y is called the approximate space of L2R2; the closed subspace sequence Wj2j∈z generated by the expansion and translation of the wavelet function ϕ1x,y,ϕ2x,y,ϕ3x,y is called the detail space of L2R2.The projection of each component of the optical flow field on wavelet basis can be obtained from formula (4), as shown in the following formulas:(5)uxx,y=∑k,m=02C−1dC,k,m1φC,k,mx,y+∑t=13∑j=CF−1∑k,m=02j−1kC,t,k,m1ϕj,k,mtx,y,(6)uyx,y=∑k,m=02C−1dC,k,m2φC,k,mx,y+∑t=13∑j=CF−1∑k,m=02j−1kC,t,k,m2ϕj,k,mtx,y.Here,dC,k,m1=ux,φC,k,m, kC,t,k,m1=ux,ϕc,k,mt, dC,k,m2=uy,φC,k,m, kC,t,k,m2=uy,ϕc,k,mt, and the subscripts k, m represent the horizontal and the vertical displacements of the wavelet basis function; J, C represent the resolution scale; and the superscripts 1, 2, and 3 represent the horizontal direction, the vertical direction, and the diagonal direction, respectively.It is known from formulas (5) and (6) that the reconstruction of the optical flow field from the wavelet basis function is a linear transformation. Therefore, the relationship between the optical flow and the wavelet basis function can be expressed by the following formula:(7)ux,y=uxx,y,uyx,y=Φx,yΘ1T,Θ2T=Φx,yΘ,(8)whereΘ=Θ1T,Θ2T,(9)Θ1=dC,k,m1,kC,t,k,m1,(10)Θ2=dC,k,m2,kC,t,k,m2,(11)Φx,y=φC,k,m,ϕj,k,m1,ϕj,k,m2,ϕj,k,m3.Here,0≤C<L,C≤j≤J,0≤k.m≤2j−1,t=1,2,3, Θ1 and Θ2 represent the wavelet transform coefficients.Using formula (7), the solution of the optical flow problem can be converted into the wavelet coefficients of the optical flow field, and the objective function can be changed into the following formula:(12)Θ∗=argminFO=argminFDI0,I1,Θ+λFRΘ,(13)ux,y=Φx,yΘ∗.Here,Θ is the unknown wavelet coefficient used to estimate ux,y.


### 3.2. Implementation

Formula (12) gives the mathematical expression of the optical flow problem based on the wavelet transform, which transforms the optical flow problem into the wavelet coefficient of the optical flow. The L-BFGS algorithm is also selected to solve the above-mentioned optimization problem. This algorithm can effectively save memory space and improve the efficiency of the operation, and is suitable for solving large-scale optimization problems. L-BFGS algorithm needs to calculate the gradient of the objective function when solving the optimization problem. The gradient of the objective function to the wavelet coefficients Θ can be expressed as(14)∇FO=∇FDI0,I1,Θ+λ∇FRΘ.In order to solve∇FO, this paper first solves ∇FD, which is expressed as(15)∇FD=∬I1x,y,Θ−I0x,y∇I1x,y,ΘΦx,ydxdy.For the solution of∇FR, Kadri-Harouna et al. [10] made a detailed derivation and explanation. In this paper, the method of discrete approximation in [10] is used to approximate FR, and the expression of ∇FR is(16)FR∼ΘΛjΘT,∇FR∼ΛjΘT,where Λj is a diagonal matrix and the diagonal element is 42j,j=C,C+1,…,L.In the existing multiscale optical flow method based on the wavelet, the optimal scale is set by prior information of the flow field or experience, or is adjusted according to the calculation results. This setting method has certain subjectivity, aiming at a certain flow field or a certain kind of flow field. In order to enhance the adaptability of the algorithm and reduce its dependence on prior information, an adaptive scale selection method is added to the algorithm. Firstly, calculate the initial displacement fieldux. Secondly, calculate the variance S1 of ux, the gradient Gx of ux, and the variance S2 of Gx. Thirdly, the threshold judgment is used to determine the optimal scale. When S1≥D1,S2≥D2, it indicates that the overall amplitude of ux changes greatly and the distribution of the flow field characteristics is complex, which requires the use of the largest scale L−1 as the optimal scale; when S1≥D1,S2<D2, it indicates that the overall amplitude of ux changes greatly, and that the distribution of the flow field characteristics is relatively simple. It is possible that ux is a gradual flow field, which requires a smaller scale L−2 as the optimal scale; when S1<D1,S2≥D2, it indicates that the overall amplitude of ux change is small, and the distribution of the flow field characteristics is complex, which requires L−1 as the optimal scale; when S1<D1,S2<D2, it indicates that the overall amplitude of ux change is small and that the distribution of the flow field characteristics is simple, which requires the use of scale L−3 as the optimal scale. As shown in Figure 1, this paper analyzes the feasibility of adaptive scale selection using 1000 images in the DNS dataset. Since the image resolution is 256∗256, L=8. For example, the threshold is set to D1=500 and D2=10.(i)
S1<D1,S2<D2, scale=L−3=5. According to Figure 1, PSNR≥30, SSIM≥0.8.(ii)
S1<D1,S2≥D2, scale=L−1=7. According to Figure 1, PSNR≥45, SSIM≥0.95.(iii)
S1≥D1,S2<D2, scale=L−2=6. According to Figure 1, PSNR≥35, SSIM≥0.90.(iv)
S1≥D1,S2≥D2, scale=L−1=7. According to Figure 1, PSNR≥45,SSIM≥0.95.Figure 1
Description of adaptive scale selection. (a) The change ofS1 (variance of ux) with image sequence. (b) The change of S2 (variance of the gradient of ux) with image sequence. (c) The corresponding relationship between the PSNR of the images restored at the different scales and the image sequence. (d) The corresponding relationship between the SSIM of images restored at different scales and the image sequence.
(a)(b)(c)(d)PSNR and SSIM represent the similarity between the image at different scales and the original image. The larger the value of PSNR and SSIM, the higher the similarity between the image at this scale and the original image. Through the above analysis, the adaptive selection of the scale can be realized by setting an appropriate threshold.In order to avoid over-smoothing of the velocity fieldux and the impact on the accuracy of the velocity field ux reconstruction, we improve the regular term parameter λ, so that the objective function can adaptively adjust the smoothness of the velocity field according to the gradient variation of ux. The expression of λ is shown in the following formula:(17)λi+1=K∗0.5−ti+11+ti+1,i≥1,(18)whereti+1=e−Δi/Δi−1,i≥1,(19)Δi=Δui22,i≥1,(20)Δ0=Δu022,where Δui is the variance of the gradient of the velocity field ux, Δu0 is the variation of the gradient of the velocity field u0x, and λ1=1,K controls the value range of λ, λ,λ∈0,0.5∗K. According to formula (17), λ is controlled by the ratio of Δui to Δui−1. When Δui is less than Δui−1, it means that the distribution characteristics of the flow field tend to be concentrated, which indicates that the phenomenon of over-smoothing occurs, and the parameter λ of the regularization term should be reduced; when Δui is greater than Δui−1, it means that the distribution characteristics of the flow field tend to be dispersed, which indicates that there is an under constraint phenomenon, and the regularization parameter λ should be increased. By adjusting λ adaptively, the objective function can reach the convergence condition smoothly.


## 3.1. Multiscale Optical Flow Representation Based on the Wavelet

I0 and I1 are images obtained by imaging the same background pattern under the same BOS experimental device; I0 is the imaging result without the flow field, and I1 is the image with the flow field. The optical flow algorithm is used to estimate the pixel displacement ux between the two images. Based on the conservation assumption of the image brightness, images I0 and I1 fulfill the following formula:(1)I1x+ux−I0x=0.In order to calculate the displacement fieldux, it is necessary to change formula (1) to the standard quadratic function, as shown in the following formula:(2)u∗=argminFd=argmin∫x∈ΩI1x+ux−I0x2dx.Here,Ω represents the set of effective pixels of image I0x.By calculating the minimum value of formula (2), the pixel displacement ux is obtained. For a point x in two-dimensional space, formula (2) is a quadratic equation of two-dimensional vector ux. However, the number of unknowns in the equation is more than that in the equation, which is an underdetermined equation. In order to better solve the above optimization problems, a convex regularization term FRu needs to be added [19]. The objective function is given by the following formula:(3)u∗=argminFu=argminFDu+λFRu.The wavelet framework introduces the decomposition of a signal into several sets of details at various scales and a remaining coarse approximation. The projection of the two components,ux, uy, of the velocity field u onto such multiresolution approximation spaces proves to be efficient within the context of optical flow estimation, enabling the handling of large displacements. Here, let’s recall the basic principle of the two-dimensional wavelet transform. For any function gx,y∈L2R2, it can be expressed linearly by Vj2j∈z,Wj2j∈z, as shown in the following formula:(4)gx,y=∑k,m=02C−1g,φC,k,mφC,k,mx,y+∑t=13∑j=CF−1∑k,m=02j−1g,ϕj,k,mtϕj,k,mtx,y.The closed subspace sequence Vj2j∈z generated by the expansion and translation of the scale function φx,y is called the approximate space of L2R2; the closed subspace sequence Wj2j∈z generated by the expansion and translation of the wavelet function ϕ1x,y,ϕ2x,y,ϕ3x,y is called the detail space of L2R2.The projection of each component of the optical flow field on wavelet basis can be obtained from formula (4), as shown in the following formulas:(5)uxx,y=∑k,m=02C−1dC,k,m1φC,k,mx,y+∑t=13∑j=CF−1∑k,m=02j−1kC,t,k,m1ϕj,k,mtx,y,(6)uyx,y=∑k,m=02C−1dC,k,m2φC,k,mx,y+∑t=13∑j=CF−1∑k,m=02j−1kC,t,k,m2ϕj,k,mtx,y.Here,dC,k,m1=ux,φC,k,m, kC,t,k,m1=ux,ϕc,k,mt, dC,k,m2=uy,φC,k,m, kC,t,k,m2=uy,ϕc,k,mt, and the subscripts k, m represent the horizontal and the vertical displacements of the wavelet basis function; J, C represent the resolution scale; and the superscripts 1, 2, and 3 represent the horizontal direction, the vertical direction, and the diagonal direction, respectively.It is known from formulas (5) and (6) that the reconstruction of the optical flow field from the wavelet basis function is a linear transformation. Therefore, the relationship between the optical flow and the wavelet basis function can be expressed by the following formula:(7)ux,y=uxx,y,uyx,y=Φx,yΘ1T,Θ2T=Φx,yΘ,(8)whereΘ=Θ1T,Θ2T,(9)Θ1=dC,k,m1,kC,t,k,m1,(10)Θ2=dC,k,m2,kC,t,k,m2,(11)Φx,y=φC,k,m,ϕj,k,m1,ϕj,k,m2,ϕj,k,m3.Here,0≤C<L,C≤j≤J,0≤k.m≤2j−1,t=1,2,3, Θ1 and Θ2 represent the wavelet transform coefficients.Using formula (7), the solution of the optical flow problem can be converted into the wavelet coefficients of the optical flow field, and the objective function can be changed into the following formula:(12)Θ∗=argminFO=argminFDI0,I1,Θ+λFRΘ,(13)ux,y=Φx,yΘ∗.Here,Θ is the unknown wavelet coefficient used to estimate ux,y.


## 3.2. Implementation

Formula (12) gives the mathematical expression of the optical flow problem based on the wavelet transform, which transforms the optical flow problem into the wavelet coefficient of the optical flow. The L-BFGS algorithm is also selected to solve the above-mentioned optimization problem. This algorithm can effectively save memory space and improve the efficiency of the operation, and is suitable for solving large-scale optimization problems. L-BFGS algorithm needs to calculate the gradient of the objective function when solving the optimization problem. The gradient of the objective function to the wavelet coefficients Θ can be expressed as(14)∇FO=∇FDI0,I1,Θ+λ∇FRΘ.In order to solve∇FO, this paper first solves ∇FD, which is expressed as(15)∇FD=∬I1x,y,Θ−I0x,y∇I1x,y,ΘΦx,ydxdy.For the solution of∇FR, Kadri-Harouna et al. [10] made a detailed derivation and explanation. In this paper, the method of discrete approximation in [10] is used to approximate FR, and the expression of ∇FR is(16)FR∼ΘΛjΘT,∇FR∼ΛjΘT,where Λj is a diagonal matrix and the diagonal element is 42j,j=C,C+1,…,L.In the existing multiscale optical flow method based on the wavelet, the optimal scale is set by prior information of the flow field or experience, or is adjusted according to the calculation results. This setting method has certain subjectivity, aiming at a certain flow field or a certain kind of flow field. In order to enhance the adaptability of the algorithm and reduce its dependence on prior information, an adaptive scale selection method is added to the algorithm. Firstly, calculate the initial displacement fieldux. Secondly, calculate the variance S1 of ux, the gradient Gx of ux, and the variance S2 of Gx. Thirdly, the threshold judgment is used to determine the optimal scale. When S1≥D1,S2≥D2, it indicates that the overall amplitude of ux changes greatly and the distribution of the flow field characteristics is complex, which requires the use of the largest scale L−1 as the optimal scale; when S1≥D1,S2<D2, it indicates that the overall amplitude of ux changes greatly, and that the distribution of the flow field characteristics is relatively simple. It is possible that ux is a gradual flow field, which requires a smaller scale L−2 as the optimal scale; when S1<D1,S2≥D2, it indicates that the overall amplitude of ux change is small, and the distribution of the flow field characteristics is complex, which requires L−1 as the optimal scale; when S1<D1,S2<D2, it indicates that the overall amplitude of ux change is small and that the distribution of the flow field characteristics is simple, which requires the use of scale L−3 as the optimal scale. As shown in Figure 1, this paper analyzes the feasibility of adaptive scale selection using 1000 images in the DNS dataset. Since the image resolution is 256∗256, L=8. For example, the threshold is set to D1=500 and D2=10.(i)
S1<D1,S2<D2, scale=L−3=5. According to Figure 1, PSNR≥30, SSIM≥0.8.(ii)
S1<D1,S2≥D2, scale=L−1=7. According to Figure 1, PSNR≥45, SSIM≥0.95.(iii)
S1≥D1,S2<D2, scale=L−2=6. According to Figure 1, PSNR≥35, SSIM≥0.90.(iv)
S1≥D1,S2≥D2, scale=L−1=7. According to Figure 1, PSNR≥45,SSIM≥0.95.Figure 1
Description of adaptive scale selection. (a) The change ofS1 (variance of ux) with image sequence. (b) The change of S2 (variance of the gradient of ux) with image sequence. (c) The corresponding relationship between the PSNR of the images restored at the different scales and the image sequence. (d) The corresponding relationship between the SSIM of images restored at different scales and the image sequence.
(a)(b)(c)(d)PSNR and SSIM represent the similarity between the image at different scales and the original image. The larger the value of PSNR and SSIM, the higher the similarity between the image at this scale and the original image. Through the above analysis, the adaptive selection of the scale can be realized by setting an appropriate threshold.In order to avoid over-smoothing of the velocity fieldux and the impact on the accuracy of the velocity field ux reconstruction, we improve the regular term parameter λ, so that the objective function can adaptively adjust the smoothness of the velocity field according to the gradient variation of ux. The expression of λ is shown in the following formula:(17)λi+1=K∗0.5−ti+11+ti+1,i≥1,(18)whereti+1=e−Δi/Δi−1,i≥1,(19)Δi=Δui22,i≥1,(20)Δ0=Δu022,where Δui is the variance of the gradient of the velocity field ux, Δu0 is the variation of the gradient of the velocity field u0x, and λ1=1,K controls the value range of λ, λ,λ∈0,0.5∗K. According to formula (17), λ is controlled by the ratio of Δui to Δui−1. When Δui is less than Δui−1, it means that the distribution characteristics of the flow field tend to be concentrated, which indicates that the phenomenon of over-smoothing occurs, and the parameter λ of the regularization term should be reduced; when Δui is greater than Δui−1, it means that the distribution characteristics of the flow field tend to be dispersed, which indicates that there is an under constraint phenomenon, and the regularization parameter λ should be increased. By adjusting λ adaptively, the objective function can reach the convergence condition smoothly.


## 4. Simulation of BOS Experiment Based on the Ray Tracing Method

### 4.1. Simulation Description

In BOS, the density gradients in a flow are measured by the apparent shift of a dot pattern viewed through a variable density medium, where the displacement is evaluated using methods similar to PIV [1, 2]. In practical application, it is necessary to use simulation data to verify the accuracy of the algorithm. To assess the accuracy of the proposed displacement estimation algorithms in the paper, synthetic BOS images are required [3], which are simulated based on the ray tracing method [20, 21]. The layout of the BOS experimental setup modeled in the simulation experiment is shown in Figure 2, and the parameters describing the placement of the elements are summarized in Table 1.Figure 2
Basic principles of BOS [3].Table 1
Parameters used to simulate the BOS experiment.
SymbolDescriptionQuantityZADistance from background pattern to flow center0.3 mZDDistance from flow center to camera0.3 mfFocal length25 mmnRefractive index distributionnx,y,z=1.000293−0.000293∗exp−x2+y2+z2RRadius of flow field0.1 mlPixel pitch5µmMImage size256 ∗ 256PDot density20 dots/32∗32 pix


### 4.2. Simulation Results

The spatial refractive index distribution function of the flow field used in the simulation experiment is shown in the following formula:(21)nx,y,z=1.000293−0.000293∗exp−x−x02+y−y02+z−z02.Here,x0,y0,z0 is the center of flow. The pixel displacements can be easily obtained by ray tracing. The simulation results are shown in Figure 3.Figure 3
Simulation results of a simple Gaussian gas flow. (a) The background pattern. (b) Distorted pattern through the flow field. (c) Optical flow results. (d) Vorticity field derived from our simulation experiment.
(a)(b)(c)(d)In order to further illustrate the effectiveness of the simulation method in this paper, PIVlab [11] and LK algorithm [22] are used to calculate the displacement distribution of the simulation image. The results are shown in Figure 4.Figure 4
Optical flow results. (a) Result of simulation. (b) Optical flow calculation results of PIVlab. (c) Optical flow calculation results of LK algorithm.
(a)(b)(c)


## 4.1. Simulation Description

In BOS, the density gradients in a flow are measured by the apparent shift of a dot pattern viewed through a variable density medium, where the displacement is evaluated using methods similar to PIV [1, 2]. In practical application, it is necessary to use simulation data to verify the accuracy of the algorithm. To assess the accuracy of the proposed displacement estimation algorithms in the paper, synthetic BOS images are required [3], which are simulated based on the ray tracing method [20, 21]. The layout of the BOS experimental setup modeled in the simulation experiment is shown in Figure 2, and the parameters describing the placement of the elements are summarized in Table 1.Figure 2
Basic principles of BOS [3].Table 1
Parameters used to simulate the BOS experiment.
SymbolDescriptionQuantityZADistance from background pattern to flow center0.3 mZDDistance from flow center to camera0.3 mfFocal length25 mmnRefractive index distributionnx,y,z=1.000293−0.000293∗exp−x2+y2+z2RRadius of flow field0.1 mlPixel pitch5µmMImage size256 ∗ 256PDot density20 dots/32∗32 pix


## 4.2. Simulation Results

The spatial refractive index distribution function of the flow field used in the simulation experiment is shown in the following formula:(21)nx,y,z=1.000293−0.000293∗exp−x−x02+y−y02+z−z02.Here,x0,y0,z0 is the center of flow. The pixel displacements can be easily obtained by ray tracing. The simulation results are shown in Figure 3.Figure 3
Simulation results of a simple Gaussian gas flow. (a) The background pattern. (b) Distorted pattern through the flow field. (c) Optical flow results. (d) Vorticity field derived from our simulation experiment.
(a)(b)(c)(d)In order to further illustrate the effectiveness of the simulation method in this paper, PIVlab [11] and LK algorithm [22] are used to calculate the displacement distribution of the simulation image. The results are shown in Figure 4.Figure 4
Optical flow results. (a) Result of simulation. (b) Optical flow calculation results of PIVlab. (c) Optical flow calculation results of LK algorithm.
(a)(b)(c)


## 5. Experimental Results and Analysis

In order to verify the effectiveness of the proposed algorithm, firstly, this paper uses the simulation data generated by the proposed simulation method to verify the effectiveness of the proposed algorithm. Secondly, the validity of the algorithm is verified by using the public DNS datasets [6, 10]. At the same time, it is compared with a variety of existing optical flow algorithms on the datasets. Finally, in order to better illustrate the effectiveness of the proposed algorithm, this paper adds a BOS experiment. Detailed analysis will be carried out below.Experiment 1.
Algorithm verification on simulation data.
The proposed method was first validated on the data generated by the simulation method in this paper. Using synthetic data does not address certain complexities of experimental data, including 3D effects, noise, and nonuniform illumination. However, these issues are intentionally avoided in the current manuscript to assess the accuracy of the method in an ideal environment where a ground truth velocity is known [6]. The results obtained by using the hybrid adaptive wavelet-based optical flow algorithm at different scales is shown in Figure 5, and Figure 6 shows the error distribution.Figure 5
Velocity fields of scales 7, 6, 5, and 4 using the hybrid adaptive wavelet-based optical flow algorithm.
(a)(b)(c)(d)Figure 6
Error distribution in different scales.
(a)(b)(c)(d)Experiment 2.
Algorithm comparison on DNS datasets.
DNS datasets are open and recognized datasets for optical flow algorithm verification. In order to further illustrate the effectiveness of the algorithm, this paper uses DNS datasets to verify the proposed algorithm and makes a comparative experiment. Firstly, it is compared with a single-scale optical flow algorithm, such as the optical flow algorithm based on cross-correlation proposed by Thielicke and Eize [11]. Secondly, it is compared with the traditional multiscale optical flow algorithm, such as the multiscale LK algorithm proposed by Qi et al. [22] or the multiscale HS algorithm proposed by Meinhardt-Llopis et al. [23]. Finally, it is compared with the wavelet-based optical flow algorithms, such as the WOF algorithm proposed by Dérian et al. [7] and the WOF-X algorithm proposed by Schmidt et al. [6]. Figure 7 shows the real optical flow estimation results and the velocity field distribution. Figure 8 is the result of the optical flow algorithm based on cross-correlation. Figures 9 and 10 are the flow estimation results of two multiscale optical flow algorithms. Figures 11–13 are the flow estimation results of three wavelet-based optical flow algorithms.
Figures11–13 show the calculation results of three wavelet-based optical flow algorithms. The coarse scale and the fine scale of these three optical flow algorithms are 3 and 6. The differences of the three algorithms are as follows: firstly, the WOF algorithm does not use the regular term constraint, the WOF-X algorithm uses the regular term constraint of fixed parameters, and the algorithm proposed in this paper uses the regular term constraint of adaptive parameters; secondly, the random values of the first two algorithms are used as the initial values of the velocity field, but the values calculated by the traditional optical flow algorithm in this paper are used as the initial values of the velocity field.Figure 7
Ground truth. (a) Optical flow field. (b) Vorticity fields.
(a)(b)Figure 8
Correlation result. (a) Optical flow field. (b) Vorticity fields. (c) Error distribution.
(a)(b)(c)Figure 9
LK result. (a) Optical flow field. (b) Vorticity fields. (c) Error distribution.
(a)(b)(c)Figure 10
HS result. (a) Optical flow field. (b) Vorticity fields. (c) Error distribution.
(a)(b)(c)Figure 11
WOF result. (a) Optical flow field. (b) Vorticity fields. (c) Error distribution.
(a)(b)(c)Figure 12
WOF-X result. (a) Optical flow field. (b) Vorticity fields. (c) Error distribution.
(a)(b)(c)Figure 13
Our method result. (a) Optical flow field. (b) Vorticity fields. (c) Error distribution.
(a)(b)(c)Experiment 3.
BOS experiment.
In order to illustrate the effectiveness of the proposed algorithm, a practical BOS experiment is carried out. The experimental equipment is shown in Figure14(a), including camera, hot air gun, and background pattern. Figure 14(b) is a reference image (background pattern). Figure 14(c) shows the experimental results (the image with the flow field).
Figure15 is the result of the optical flow algorithm based on cross-correlation. Figures 16 and 17 are the flow estimation results of two multiscale optical flow algorithms. Figures 18–20 are the flow estimation results of three wavelet-based optical flow algorithms.Figure 14
Explanation of actual BOS experiment. (a) Experimental setup. (b) Background pattern. (c) Experimental result.
(a)(b)(c)Figure 15
Correlation result. (a) Vorticity fields. (b) Motion compensated image. (c) Error distribution between the compensation image and the image with the flow field.
(a)(b)(c)Figure 16
HS result. (a) Vorticity fields. (b) Motion-compensated image. (c) Error distribution between the compensation image and the image with the flow field.
(a)(b)(c)Figure 17
LK result. (a) Vorticity fields. (b) Motion-compensated image. (c) Error distribution between the compensation image and the image with the flow field.
(a)(b)(c)Figure 18
WOF result. (a) Vorticity fields (b) Motion-compensated image (c) Error distribution between the compensation image and the image with the flow field.
(a)(b)(c)Figure 19
WOF-X result. (a) Vorticity fields. (b) Motion-compensated image. (c) Error distribution between the compensation image and the image with the flow field.
(a)(b)(c)Figure 20
Our method result. (a) Vorticity fields. (b) Motion-compensated image. (c) Error distribution between the compensation image and the image with the flow field.
(a)(b)(c)


### 5.1. Analysis of Experimental Results

In this paper, root-mean-square endpoint error (RMSE) is used to quantitatively evaluate the experimental results [6], as shown in the following formula:(22)RMSE=1N∑urefx,y−ux,y2.The lower the values of RMSE, the better the displacement filed. The evaluation results of Experiment1, Experiment 2, and Experiment 3 are shown in Table 2, Table 3, and Table 4, respectively.Table 2
The evaluation results of Experiment1.
ErrorScale7654RMSE0.0127430.0154000.0333940.064844Table 3
The evaluation results of Experiment2.
ErrorMethodCorrelationMultiscale HSMultiscale LKWOFWOF-XOur methodRMSE0.284640.0681620.263130.156040.137970.038876Table 4
The evaluation results of Experiment3.
MethodCorrelationMultiscale HSMultiscale LKWOFWOF-XOur methodRMSE6.72666.18746.42887.07616.3746.1865PSNR19.620220.345920.013519.180220.087920.3472SSIM0.936140.945110.941070.925720.942110.94514Comparing the RMSE values of different scales, as shown in Table2, it is obvious that the larger the scale, the higher the accuracy of the hybrid adaptive wavelet-based optical flow algorithm. This also accords with the principle of the wavelet transform. The larger the scale, the more abundant are the details contained in the wavelet coefficients. To ensure the accuracy and improve the efficiency of the algorithm, it is necessary to choose the appropriate resolution scale. In order to illustrate the effect of the resolution scale on the accuracy of the algorithm, the results of four scales are presented in Table 2. The optimal resolution scale calculated by the adaptive scale selection method proposed in this paper is scale = 6. Table 2 shows that the accuracy of the algorithm in scale = 6 is similar to that in scale = 7, but the calculation amount of the algorithm is reduced by four times. However, the accuracy of the algorithm in scale = 6 is twice as high as that in scale = 5, and four times as high as that in scale = 4. This shows that the adaptive scale selection method is effective and feasible. This not only ensures the accuracy of the algorithm but also improves the efficiency of the algorithm. If the error allows, a lower resolution scale can be selected, such as scale = 5 or scale = 4.Comparing the RMSE values of different algorithms, as shown in Table2, the RMSE value of our method is the smallest. As can be seen from Figure 11(c), the error of the WOF-X algorithm in some regions with a large refractive index gradient is large, which may be caused by the inability of the fixed regularization constraint parameters to adapt to some regions. Compared with the WOF-X algorithm, the algorithm proposed in this paper has higher calculation accuracy and smaller calculation error in the whole image area, as shown in Figure 12(c). This shows that the adaptive regularization constraint method proposed in this paper is effective and feasible.In order to observe the accuracy of each algorithm more intuitively, this paper extracts the 50th, 100th, 150th, and 200th lines of the estimated velocity field of each algorithm and compares them with the corresponding lines of the true velocity field, as shown in Figure21. According to Table 5, the velocity field estimated by the HS algorithm, the WOF-X algorithm, and the algorithm proposed in this paper has a high coincidence with the true velocity field. Through the above analysis, it is easy to draw the following conclusion: the accuracy of our algorithm is higher than that of other algorithms mentioned in this paper.Figure 21
Results from the comparison algorithms compared to the true vorticity fields for one component of velocity taken over four horizontal slices through the domain.
(a)(b)(c)(d)Table 5
The evaluation results of the 50th, 100th, 150th, and 200th lines.
RMSE ErrorMethodCorrelationMultiscale HSMultiscale LKWOFWOF-XOur method50th0.013510.003380.009150.012470.005370.00079100th0.010060.003800.008040.010240.005190.00064150th0.012180.003730.008340.014330.006650.00233200th0.008910.003250.008780.018350.006590.00154In this paper, the effectiveness of the proposed algorithm is verified by the actual BOS experiment, and the motion-compensated image is reconstructed according to the calculated displacement fieldux, as shown in Figure 16(b). The motion compensation image was compared with the experimental image (image with flow field， as shown in Figure 14(c)) by RMSE, PSNR (Peak Signal to Noise Ratio), and SSIM (structural similarity). PSNR is the most common and widely used objective image evaluation index. However, it is based on the error between corresponding pixels, that is, the image quality evaluation based on error sensitivity. The larger the PSNR value, the smaller the distortion. SSIM, an image quality evaluation index, evaluates image quality based on structural similarity. It measures the image similarity from three aspects: brightness, contrast, and structure. The higher the SSIM value, the higher the image similarity. Comparing the RMSE values, the PSNR values, and the SSIM values of different algorithms, as shown in Table 4, the RMES value of our method is the smallest and the PSNR and SSIM values of our method are the highest. According to the RMSE, PSNR and SSIM values in Table 4, the proposed method is superior to other comparison methods.


## 5.1. Analysis of Experimental Results

In this paper, root-mean-square endpoint error (RMSE) is used to quantitatively evaluate the experimental results [6], as shown in the following formula:(22)RMSE=1N∑urefx,y−ux,y2.The lower the values of RMSE, the better the displacement filed. The evaluation results of Experiment1, Experiment 2, and Experiment 3 are shown in Table 2, Table 3, and Table 4, respectively.Table 2
The evaluation results of Experiment1.
ErrorScale7654RMSE0.0127430.0154000.0333940.064844Table 3
The evaluation results of Experiment2.
ErrorMethodCorrelationMultiscale HSMultiscale LKWOFWOF-XOur methodRMSE0.284640.0681620.263130.156040.137970.038876Table 4
The evaluation results of Experiment3.
MethodCorrelationMultiscale HSMultiscale LKWOFWOF-XOur methodRMSE6.72666.18746.42887.07616.3746.1865PSNR19.620220.345920.013519.180220.087920.3472SSIM0.936140.945110.941070.925720.942110.94514Comparing the RMSE values of different scales, as shown in Table2, it is obvious that the larger the scale, the higher the accuracy of the hybrid adaptive wavelet-based optical flow algorithm. This also accords with the principle of the wavelet transform. The larger the scale, the more abundant are the details contained in the wavelet coefficients. To ensure the accuracy and improve the efficiency of the algorithm, it is necessary to choose the appropriate resolution scale. In order to illustrate the effect of the resolution scale on the accuracy of the algorithm, the results of four scales are presented in Table 2. The optimal resolution scale calculated by the adaptive scale selection method proposed in this paper is scale = 6. Table 2 shows that the accuracy of the algorithm in scale = 6 is similar to that in scale = 7, but the calculation amount of the algorithm is reduced by four times. However, the accuracy of the algorithm in scale = 6 is twice as high as that in scale = 5, and four times as high as that in scale = 4. This shows that the adaptive scale selection method is effective and feasible. This not only ensures the accuracy of the algorithm but also improves the efficiency of the algorithm. If the error allows, a lower resolution scale can be selected, such as scale = 5 or scale = 4.Comparing the RMSE values of different algorithms, as shown in Table2, the RMSE value of our method is the smallest. As can be seen from Figure 11(c), the error of the WOF-X algorithm in some regions with a large refractive index gradient is large, which may be caused by the inability of the fixed regularization constraint parameters to adapt to some regions. Compared with the WOF-X algorithm, the algorithm proposed in this paper has higher calculation accuracy and smaller calculation error in the whole image area, as shown in Figure 12(c). This shows that the adaptive regularization constraint method proposed in this paper is effective and feasible.In order to observe the accuracy of each algorithm more intuitively, this paper extracts the 50th, 100th, 150th, and 200th lines of the estimated velocity field of each algorithm and compares them with the corresponding lines of the true velocity field, as shown in Figure21. According to Table 5, the velocity field estimated by the HS algorithm, the WOF-X algorithm, and the algorithm proposed in this paper has a high coincidence with the true velocity field. Through the above analysis, it is easy to draw the following conclusion: the accuracy of our algorithm is higher than that of other algorithms mentioned in this paper.Figure 21
Results from the comparison algorithms compared to the true vorticity fields for one component of velocity taken over four horizontal slices through the domain.
(a)(b)(c)(d)Table 5
The evaluation results of the 50th, 100th, 150th, and 200th lines.
RMSE ErrorMethodCorrelationMultiscale HSMultiscale LKWOFWOF-XOur method50th0.013510.003380.009150.012470.005370.00079100th0.010060.003800.008040.010240.005190.00064150th0.012180.003730.008340.014330.006650.00233200th0.008910.003250.008780.018350.006590.00154In this paper, the effectiveness of the proposed algorithm is verified by the actual BOS experiment, and the motion-compensated image is reconstructed according to the calculated displacement fieldux, as shown in Figure 16(b). The motion compensation image was compared with the experimental image (image with flow field， as shown in Figure 14(c)) by RMSE, PSNR (Peak Signal to Noise Ratio), and SSIM (structural similarity). PSNR is the most common and widely used objective image evaluation index. However, it is based on the error between corresponding pixels, that is, the image quality evaluation based on error sensitivity. The larger the PSNR value, the smaller the distortion. SSIM, an image quality evaluation index, evaluates image quality based on structural similarity. It measures the image similarity from three aspects: brightness, contrast, and structure. The higher the SSIM value, the higher the image similarity. Comparing the RMSE values, the PSNR values, and the SSIM values of different algorithms, as shown in Table 4, the RMES value of our method is the smallest and the PSNR and SSIM values of our method are the highest. According to the RMSE, PSNR and SSIM values in Table 4, the proposed method is superior to other comparison methods.


## 6. Conclusions

In order to improve the accuracy of the BOS device in transparent flow field measurement, this paper presents the hybrid adaptive wavelet-based optical flow algorithm for Background Oriented Schlieren (BOS) experiments. The algorithm can effectively estimate the displacement field. Firstly, this paper simulates the BOS experiment and uses the simulation data to verify the effectiveness of the proposed algorithm. In Experiment1, simulation data are used to verify the proposed algorithm, and the influence of the resolution scale on the accuracy and efficiency of the proposed algorithm is discussed. The results of Experiment 1 show that the algorithm proposed in this paper is effective. Secondly, in order to further prove the effectiveness of the algorithm proposed in this paper, in Experiment 2, this paper uses DNS datasets to verify the algorithm proposed in this paper and compares it with the existing related algorithms. The results of Experiment 2 show that the accuracy of this algorithm is better than the other algorithms mentioned in this paper. Especially, in the region of large refractive index gradient, the accuracy of the algorithm proposed in this paper is higher than that of other comparison algorithms. It also shows that the proposed algorithm is effective and feasible. Lastly, the method proposed in this paper is verified by the actual BOS experiment. The results of Experiment 3 show that the accuracy of this algorithm is better than other algorithms mentioned in this paper.



---

*Source: 5138153-2020-09-22.xml*