Assessment and application of optical flow in background-oriented schlieren for compressible flows

Optical flow provides an opportunity to elevate the resolution and sensitivity of deflection sensing in background-oriented schlieren (BOS). Despite extensive relevant literature within the field of computer vision, there is a lack of proper quantification of its abilities and limitations with regard to the state-of-the-art BOS experiments. Thus, this study performs an assessment of accuracy and resolution limits in different flow field scenarios utilizing background patterns generated with random dot and wavelet noise distributions. Accordingly, a synthetic assessment over a theoretically generated Prandtl–Meyer expansion fan is conducted with variations introduced in the background patterns and operational parameters of optical flow. A clear superiority of accuracy and resolvable range of density gradient amplitudes over cross-correlation is demonstrated. Moreover, an experimental assessment of supersonic flow features over multiple wind tunnel models is performed. The influence of experimental constraints, limitations and uncertainties related to the application of optical flow in BOS and its comparative performance against the block-matching counterpart is characterized.


Introduction
Background-oriented schlieren (BOS) imaging is of great interest to fluid dynamics researchers enabling quantitative characterization of density varying flows by means of an inexpensive and simple setup that can be used for both two and three-dimensional applications. It is performed by acquiring a pair of images with or without the presence of density variations occurring in a flow field. These images contain randomly distributed illumination intensities whose positions on the image plane are modified by the presence of a density varying medium between the recorded scene and the image acquisition system (Settles 2001). As the pair of images is processed, the observed variations in the background pattern are traced by means of various digital image processing approaches (Hargather et al. 2011). Thus, the requirements of a BOS setup can be as simple as a single camera (Hargather and Settles 2010), which makes it a unique measurement technique to obtain a direct quantification of flow features with high versatility and applicability to a wide range of conditions (Raffel et al. 2000).
However, BOS is associated with a major sensitivity limitation originated by an inherent focusing problem (Raffel 2015). Considering a generic optical configuration for a BOS setup, increasing the sensitivity demands an amplification of the distance between the object and the background while minimizing the distance between the object and the camera. However this is constrained due to two main reasons. First, having the camera close to the object causes divergence of the light rays to dominate the measured displacements related to the refractive index variations within the test section which requires corrections and calibration (Elsinga et al. 2004). Secondly, since the background needs to be in focus to accurately measure the displacements, the object and the associated flow features of interest remain out of focus. Hence, increasing the distance between the background and the flow field amplifies the geometric blur. In combination with the diffraction limited minimum image diameter, it reduces the sensitivity and the physical resolution of the system, severely affecting the accuracy of the measurements (Gojani et al. 2013). Thus, an inevitable trade-off between the sensitivity and resolution of the flow structures that can be resolved with BOS exists from an optical stand point.
The sensitivity and resolution of a BOS system are also governed by the image processing approach utilized to characterize the deviations of the light rays from their original paths. The most common method employed for this purpose is the cross-correlation-based block matching algorithm ( Fig. 1, left) developed for particle image velocimetry (Westerweel 1997). The underlying theory of spatial cross-correlation computes a measure of similarity of two intensity (illumination) maps as a function of displacement with respect to the original position of the selected blocks of pixels referred to as interrogation windows (IW) (Loeb et al. 1983). Hence, the corresponding displacement vectors represent an average refractive index variation within the IW (Keane and Adrian 1992). However, as the minimum size of the IW governed by the necessary signal-to-noise ratios of intensity maps for higher correlation coefficients to be established, resolution of the resulting vector field becomes limited. This limitation is not only reduces the resolution, but also yields the gradient information being prone to oversmoothening due to the averaging effects (Hargather and Settles 2012). Even though multi-pass interrogation schemes exist to enable a larger dynamic range of pixel displacements to be covered (Scarano and Riethmuller 2000), these approaches are developed for tracing fluid particle motions within a bulk flow which makes them ineffective against very localized flow features of density variations. Thus, the dependence on neighboring pixel displacements and the documented degradation of accuracy in presence of large displacement gradients (Kähler et al. 2012) further alleviates the accuracy and resolution of BOS measurements when processed with block-matching approaches such as cross-correlation.
Nevertheless, the problem of inferior sensitivity and resolution can be addressed through the application of optical flow (OF) image processing technique (Fig. 1, right) which defines the motion of objects within the image plane by means of brightness patterns (Davies 2012). The biggest advantage of using OF for BOS can be referred to as its capability to capture the displacement of each pixel enabling one vector-per-pixel resolution. This is proven to enable OF go beyond the resolution and accuracy characteristics of PIV algorithms for reconstructing density field information (Atcheson et al. 2009). Thus, leveraging from its elevated sensitivity and resolution characteristics, various researchers have opted for optical flow to quantify the displacements due to varying refractive index medium using BOS. Accordingly, Nicolas et al. (2016) performed 3D reconstruction of a supersonic jet flow employing the Lucas-Kanade optical flow approach which relies on a window-based neighborhood dependency. Moreover, Hayasaka et al. (2016) used a variational formulation of optical flow to characterize the laser induced underwater shock waves where the residuals of brightness constancy assumption (non-zero divergence) are utilized to detect the shockwave shape and location. Furthermore, optical flow as the processing tool of choice has been applied to combustion studies for the tomographic reconstruction of flames (Grauer et al. 2018).
Furthermore, as the BOS measurements are increasingly employed for the development and validation of computational fluid dynamics models (Ramanah et al. 2007;Kirmse et al. 2011), the cross-correlation-based BOS becomes vulnerable to fall short of the corresponding resolution and accuracy requirements, whereas OF is showcased to offer a viable solution. Despite its utilization in various occasions, there exists a lack of quantitative characterization of its capabilities and application details especially in high-speed flows. Although Heineck et al. (2021) performed natural BOS through in-flight measurements using OF, its performance is demonstrated only qualitatively against cross-correlation. Considering the complexity of the models exposed to high speed flows, intensity of the localized density gradients and the varying scales of flow structures, this study showed that OF can measure high fidelity density field data from quantitative visualization of supersonic/hypersonic flows. Yet, its utilization requires proper quantification of the resolution, sensitivity and accuracy of OF-based BOS as well as the optimization of the experimental and image processing setups.
This study aims to quantify the accuracy characteristics of optical flow as the deflection sensing stage of BOS compared to cross-correlation. It starts with a brief explanation of the OF algorithms in Sect. 2.1 and gives a description of suitable background pattern candidates for generic BOS applications that can be tailored towards the use of OF in Sect. 2.2. Then, a synthetic assessment of various deflection sensing schemes for BOS is performed utilizing a theoretical Prandtl-Meyer (PM) expansion fan in Sect. 3. Finally, the experimental applicability of the deflection sensing and backgrounds generation methods under investigation is demonstrated in Sect. 4 by conducting BOS measurements of Mach 3.5 supersonic flows over a von Karman Ogive and a 10° wedge. The corresponding experimental setup and the data processing procedures are detailed in Sect. 4.1 and Sect. 4.2, respectively, and results are discussed in Sect. 4.3.

Optical flow
Various approaches to OF analysis can be broadly distinguished to have a common gradient-based scheme that is associated with two main assumptions of brightness constancy and spatial smoothness of the intensity distributions (Sun et al. 2008). The combination leads to the governing equation of motion of illumination intensities within the image plane (defined as the optical flow) in relation to the variations of illumination distributions between the two image sequences. Hence, a constancy of illumination intensity can be established which suggests the intensity of pixels moving within the image plane is constant (Kearney et al. 1987) . On the other hand, the assumption of small variation of optical flow vectors suggests that a 1st order Taylor series expansion (TSE) can be used since 2nd order gradients of illumination intensity variations are considered to be negligibly small. Accordingly, brightness constancy suggests the variation of local pixel intensity is directly related to the motion of illumination gradients within the image plane (Eq.1).
Hence, the displacement ( x, y ) of normalized illumination intensities (I) can be represented as functions of the temporal ( I∕ t ) and spatial pixel intensity gradients ( I∕ x , I∕ y ). For BOS applications, the temporal derivative of illumination intensities corresponds to the difference between the disturbed and reference image ( I∕ t=(I n − I 0 )∕Δt as shown in Fig. 1), while the nominal time interval is unitary since the difference between I n and I 0 is independent of time ( Δt =1). Nevertheless, the underdetermined structure of Eq.1 requires a problem closure to obtain a unique solution which leads to the well-known aperture problem (Galvin et al. 1998). Hence, gradient-based approaches contain additional energy functions (Eq.2) to close the system which include a penalization ( D , S ) for the displacement vector residuals, enforcing brightness constancy. Accordingly, Horn and Schunck (1981) proposed the use of a quadratic penalty function ( D ) for the displacement vectors corresponding to the deviations from the brightness constancy assumption. Moreover, Eq.2 included a quadratic penalization ( S ) of the displacement vector gradients to increases the robustness of the algorithm against spurious displacement vectors ( Barron et al. 1994). The weighting between the contribution of the brightness constancy violations and the vector gradients are controlled by a user defined regularization parameter, . However, this global penalization approach enforces the displacement vectors to vary smoothly over the reconstructed domain. Thus, precise selection of its value is crucial for an accurate displacement field estimation which would enable suppression of spurious vectors while preserving physical displacement magnitudes throughout the image plane (Nagel 1983). (1) Page 4 of 20 The biggest challenge of achieving high-resolution displacement fields with optical flow originates from the fact that practical image sequences contain multiple objects in motion independently. This causes occurrence of a global incoherence in terms of the acquired image intensity variations that also contain discontinuities (Deriche et al. 1996). Hence, as the TSE is employed to describe the nonlinear intensity variations based on a brightness constancy assumption, the linearized problem is only valid within the description of TSE. In other words, the derived set of governing equations is only valid for small displacements in absence of any discontinuities. Therefore, additional methods are required to relax the constraints on the TSE-based formulation of the optical flow problem that would allow large displacements within the image plane to be reconstructed accurately (Adelson and Bergen 1985). Accordingly, hierarchical algorithms are proposed to separate the computation of displacement vectors based on the scales of optical flow contained within the processed images (Anandan 1989). The acquired images are assumed to contain a wide spatial-frequency spectrum where the high amplitude (low spatial-frequency) intensity variations are represented by the large scale motion within the image plane, while low amplitude (high spatial-frequency) intensity variation provides the detailed structures over a smaller range. Nevertheless, identification of these scales and separating them accordingly requires a proper control strategy since spatial-frequency decomposition based on the sampling rates of the lower-frequency information would lead to a loss of resolution (Burt 1981). Such a proper decomposition of the image intensities leads to a pyramid scheme that is composed of different spatialfrequency channels (Tanimoto and Pavlidis 1975). The pyramid is generated by consecutive downsampling of the original images by a scaling factor (constant throughout the pyramid) until a certain resolution level based on the anticipated scale of largest motion is reached (Fig. 2, left panel). Then, the displacement detection process begins by applying the optical flow at the smallest resolution image (lowest spatial-frequency component) and the displacement field is upsampled to the resolution of the next pyramid level to be subtracted from the image intensity variation by unwarping the image (Fig. 2, right panel). The process is repeated until the pyramid level with the original acquired image resolution is reached.

BOS background generation
The conventional background patterns in BOS are adjusted to suit the needs of the cross-correlation approach, have a high SNR in correlation maps, which yielded random dot patterns based on the optical configuration of the measurement system (Fig. 3, left). Accordingly, concentration, size, intensity profile and interparticle distance of the random dots are specified to maximize contrast, minimize peak-locking and maximize subpixel interpolation accuracy (Adrian and Westerweel 2011). Alternatively, in line with the sequential downsampling of images following a coarse-to-fine resolution layout, multi-scale procedural noise patterns can be used as the background image. Thus, the short comings of random noise functions can be avoided which often requires a match between the camera resolution and the background pattern for optimum performance (Yang et al. 2014).
Hence, the wavelet noise (WN) pattern proposed by Cook and DeRose (2005) is utilized to construct a background image where the scaling factors of the hierarchical OF algorithm is matched to capture both large and small scale structures with high accuracy. The content stored within these WN patterns is controlled by means of two parameters, 2 n I (initial resolution) and 2 n F (final resolution). The lower spatial frequency content is modulated by the lower resolution levels where large scale displacements are resolved. This also acts as an anti-aliasing filter when the hierarchical approaches of the OF downsamples the high-resolution image to a lower resolution by preserving orthogonality. Combined with the high frequency content governed by the final resolution specification, each resolution level of noise represents a bandwidth limited procedural texture that can be tailored to the specific spatial resolution needs of the experimental setup. Accordingly, Fig. 3 shows three different configurations of WN patterns. Absence of high-frequency content, when the final resolution level ( n F ) is decreased, is shown in Fig. 3 (mid-left) which disturbs the orthogonality of the noise pattern while removing the low-frequency content. Contrarily, increasing the initial resolution level ( n I ), amplifies the standard deviation of the intensity distribution where vulnerability against aliasing is increased (Fig. 3, right).

Synthetic assessment
Comparative performance characterizations of cross-correlation (XCOR) and OF algorithms with different background patterns is performed using a theoretical expansion fan over which the density variations are computed by means of the Prandtl-Meyer (PM) function (Anderson 2001). The flow scenario is comprised of a 10° wedge exposed to supersonic flow at M ∞ =3.5 (Fig. 4). Variation of flow properties over the leading edge shock are calculated using the oblique shock relations which yielded the flow conditions upstream of the expansion fan (Anderson 2001). The region of interest for the synthetic assessment is selected slightly away from the model surface in order to prevent the optical distortion effects due to strong density gradients compromising the integrity of cross-correlation maps (Elsinga et al. 2005). Accordingly, the density gradients ( Fig. 6, top-middle) within the ROI is computed from the theoretical density field (Fig. 6,. Then, the ray displacements due to the density variations are rendered via a ray-tracing-based image generation methodology (Rajendran et al. 2019). The synthetic optical setup is configured such that a maximum displacement value of (‖( x, y)‖ 2 ) max < 8 pixels is obtained in order to match the dynamic range of the XCOR (Fig. 6, top-right).
In consideration of the well-established rules of crosscorrelation, the random dot (RD) backgrounds are generated for a high-accuracy reconstruction of displacements using IW sizes of 16 × 16 pix 2 . Hence, the random dot backgrounds are generated with particle diameter of 2 pixels ( d p = 2) and an image particle density of N I = 0.027 ppp. The dots are arranged to have a Gaussian intensity distribution (Fig. 5, top-left) for high accuracy sub-pixel interpolations which is also in agreement with the spatio-filtering requirements of optical flow (Christmas 2000). Furthermore, in order to investigate the influence of various spatial frequency ranges contained in the WN patterns, different backgrounds with varying high-and low-frequency levels are generated as shown in Fig. 3. The influence of random distribution of dots and gray scale levels of the random dot (RD) and wavelet noise (WN) backgrounds respectively on the uncertainty level is quantified by sampling each background pattern 50 times. This is achieved by keeping the same size ( d p ) and image particle density ( N I ) for the random dot backgrounds and the same initial ( n I ) and final ( n F ) resolution levels for the wavelet noise background constant for each group of 50 samples. Moreover, the synthetic images are provided with an additive Gaussian noise (generated via a build-in noise generation function MATLAB) to stimulate experimental imperfections (Fig. 5, following image warping with the displacement fields computed from the PM theory (Cai et al. 2022). This enables the camera noise to be simulated which is varying between the reference and wind-on images.
While both types of backgrounds are utilized for assessing the accuracy of OF, only the RD backgrounds are used with XCOR due to the downgraded SNR of the correlation maps caused by the reduced contrast of the WN patterns (Atcheson et al. 2008). The Horn-Shunck OF (described in Sect. 2.1) is employed through an open-source computer vision (CV) code by Sun et al. (2014). The hierarchical layout is constructed using 6 pyramid levels to minimize the interpolation errors between different levels (Hofinger et al. 2020) and 10 intermediate warping steps based on a convergence criteria of 10 −3 of the intermediate energy function residuals. For each background configuration, the regularization parameter is varied over a range (0.01≤ ≤ 250) to demonstrate the influence of different weighting schemes. Moreover, XCOR is performed using a multi-pass approach where 3 passes with consecutive reduction of IW sizes (64 × 64, 32 × 32, 16 × 16 pix 2 ) is employed with 50% overlap (Thielicke and Sonntag 2021). Both methods are accompanied with a median filtering-based outlier detection and removal scheme to eliminate spurious vectors (Tan and Jiang 2019). The median filter identifies the median image intensity value of the pixels within a kernel of user defined size and assigns the central pixel of interest the intensity value of the median of the neighboring kernel. This is performed by scanning the complete image with the median filter windows to prevent outlying signal entries acquired by the cameras in terms of image illumination intensities to be removed.
The obtained results using various methods are compared both qualitatively and quantitatively over the average endpoint errors (AEPE) is defined as the Euclidean norm of the error vector computed by performing a vectorial subtraction of the reconstructed displacement vector from the ground truth displacement vector at each pixel (Baker et al. 2007). The statistical data sampling for arithmetic averaging and computation of standard deviation of the average endpoint error (AEPE) is done over the entire region where the displacements due to the expansion fan are imposed. This prevents the damping of the error values by the displacement vector with zero synthetic density gradients. The AEPE values are quantified within the region of interest over 50 image samples comparing the arithmetic mean and median. The arithmetic mean and median differ only by 0.01 pixels, indicating that smearing effects on the sharp gradients associated with the arithmetic averaging remain relatively small for the investigated datasets. The error quantification is therefore provided via the arithmetic mean values.
To start with the qualitative assessments, the resolution superiority (0.016 vector/pixel for XCOR vs. 1 vector/pixel for OF) of OF inherently allows higher sensitivity in terms of the gradients of displacement vectors. Regardless of the background type, OF is able to capture larger displacement regions with improved accuracy as XCOR suffers from the averaging effects over the locally confined presence of density gradients which causes the displacements to be dissipated over larger regions (Fig. 6, bottom). Furthermore, comparing the performance of OF with RD and WN backgrounds, a very moderate increased level of noise is noticeable with the use of RD backgrounds. This is related to the fact that, in absence of illumination distributions induced by dots, OF is left with the low level additive noise on the images to trace the displacements. Thus, WN backgrounds covering the entire image plane with a controlled intensity variation prevents the dependence of the reconstruction procedure on the random distribution of dots and increases the consistency of the reconstruction accuracy. Nevertheless, the smooth variation of density gradients away from the first and last Mach waves that initiates and terminates (respectively) the expansion fan, suppresses these effects and yields very similar results to that of OF-based deflection sensing in comparison to XCOR.
Characterizing the differences obtained with XCOR and OF (with both WN and RD backgrounds) quantitatively, Fig. 7 shows the displacement magnitudes over three lines (Y = 20, 70 and 120 pix) reconstructed by means of different approaches. The accuracy increase obtained with OF concentrates around the regions where displacement gradients ( 2 ∕ x 2 , 2 ∕ y 2 ) are highest in terms of magnitude. Although the neighborhood dependency of OF imposed by the non-zero regularization parameter ( ) prevents the discontinuous gradients to be captured, the improved resolution allows much better localization of the expansion fan boundaries compared XCOR. Moreover, especially closer to the model where the magnitude of density gradients increase drastically, the severely localized high intensity features cannot be computed by the XCOR due to the limitations imposed by the limited dynamic range and resolution characteristics (Fig. 7, bottom). On the other hand, OF with both backgrounds shows a superior capability of having the same level of accuracy throughout the entirety of the spatial frequency spectrum when a hierarchical reconstruction scheme is used. Thus, computation of AEPE over the identified lines further supports these results where accuracy superiority of OF over XCOR is manifested with an average of 55% reduction in error magnitudes with both background patterns (Table 1). Moreover, dominance of the processing algorithm in comparison to the background is further indicated by Table 1, where a strong match between the error levels obtained with OF using both RD and WN backgrounds is present.
Furthermore, a parametric study on the influence of background patterns is performed with varying initial ( n I ) and final ( n F ) resolution levels of the WN pattern. First, the initial spatial resolution level is kept constant at the minimum level ( n I =2) while the final resolution level is modified (5≤ n F ≤10). Due to the wide spectrum of displacement vector magnitudes computed from the PM theory, an increase of accuracy with increasing n F is noted. The minimum error for reconstruction is obtained at the level of n F =8 which matches the resolution of the reference displacement field (Fig. 8, left). As the final resolution level is increased beyond the resolution of the flow field, the contribution of the displacements caused by the flow field are suppressed. Hence, increasing the regularization parameter shifted the influence of displacement gradients to the lower end of the spatial frequency spectrum which yielded lower reconstruction accuracy for n F ≥ 9. Then the procedure for background generation is inverted by keeping the final resolution level at n F =8 and varying the initial resolution level over a range of 2 ≤ n I ≤ 5 (Fig. 8, right). The impact of modifying the lower end of the spatial frequency spectrum is observed to be significantly small compared to the n F which is attributed to the characteristics of the reference case where density gradients are considerably smooth apart from the extremities of the expansion fan. Thus, having the correct range of resolution ( n I =2 & n F =8) levels where the motion over image plane is dominated by the displacement vectors, allows OF to resolve the full spectrum of spatial frequency domain. Moreover, increasing the initial resolution level influences the uncertainty levels of the OF reconstructions by increasing the possibility of aliasing due to the removal of lower frequency content.
Analyzing the influence of regularization parameter on the OF reconstructions using RD backgrounds, higher levels of regularization parameters are required to achieve higher accuracy. This is also related to the same contrast issue over the regions that suffer from the scarcity of dots where the considerably low level of illumination intensities (associated with the additive noise) provides the only features that OF can use to trace the motion of the background pattern. Accordingly, the spatial accuracy consistency of WN backgrounds over RD backgrounds is denoted in Fig. 8 in terms of error distributions throughout the displacement field. The error levels are not only lower with WN backgrounds but they are also confined to a lower range with varying regularization parameter compared to RD backgrounds. Thus, a high confidence reconstruction of ray displacements over a wider range of density gradient amplitudes can be established with optical flow especially coupled with a procedural noise pattern that contains multiple levels of spatial-frequency scales. This is also a major improvement in terms of reducing the sensitivity of the OF algorithms to the regularization parameter whose proper selection poses a very common challenge (Nicolas et al. 2015;Lu et al. 2021).
Finally, it should be noted that the computational time required for the cross-correlation (XCOR) and optical flow (OF) algorithms to complete the processing of one image pair is 18 seconds and 42 seconds respectively. Nonetheless, the OF algorithm that provides higher accuracy reconstruction of displacement fields achieve statistical converge with lower number of images. The evolution of mean AEPE values over the sampling domains of increasing number of  images reveals that XCOR requires ∼ 30 images to converge to a constant mean value whereas OF takes ∼ 15 images both for random dot and wavelet noise backgrounds (Fig. 9).
Therefore, while the computational cost of processing a single image pair is higher for the OF algorithm, the overall time effort to produce statistically converged results is comparable for the two techniques.

Experimental setup
The experimental setup employed in this work is implemented on the S-4 supersonic wind tunnel of the von Karman Institute (VKI). The S-4 rig is a blow-down wind tunnel equipped with a 2D contoured nozzle providing Mach 3.5 supersonic flow into a test section of 8 × 10 cm 2 . The reservoir pressure can be adjusted within a range of 3 to 18 bars to achieve a typical unit Reynolds number of 5 × 10 7 /m. Configuration of the BOS setup is composed of a Hamamatsu ORCA-Flash 4.0 V3 Digital CMOS camera, a Nikkor entocentric lens of focal length f = 200 mm, four Chazon 100W 10 × 10 blue (465 nm) LED arrays as the light source and backgrounds laser printed on 240 g/m 2 paper sheets (Figs.10, 11). The Hamamatsu ORCA-Flash 4.0 V3 digital camera has a scientific CMOS image sensor of size 2048×2048 pixels with 33000:1 dynamic range and equipped with a 16-bit analog to digital converter. The LED arrays are mounted on a brass plate which is actively air-cooled throughout the experiments. The LEDs, supplied with a DC power source, are utilized to illuminate the backgrounds from behind to enhance the recorded contrast and eliminate reflection effects. The distance between the light source and the background is adjusted to optimize the tradeoff between light intensity and uniformity. The supersonic tests are conducted with two different models of a von Karman ogive and a 10° wedge (end-to-end) at zero angle of attack to have a wide range of density gradients as well as comparing two-and three-dimensional flow features. The optical configuration of the BOS setup is adjusted by setting the main distances of the background ( Z D ) and the camera ( Z A ) from the schlieren object (Fig. 13, 1st stage). The distance between the object and the background reveals an optimization problem on its own which requires larger Z D for higher sensitivity, whereas the defocusing problem dictates the background should be close enough to the object to minimize the overall image blur ( d Σ ) and the diffraction effects in close proximity of the model (Raffel 2015). Moreover, in order to preserve the validity of the paraxial approximation while relating the displacement vectors to the density gradients, the background should be separated sufficiently from the object such that the width of the schlieren object needs to be smaller than the distance between the object and the background (Moisy et al. 2009). Thus, the background is placed 0.26 m from the wind tunnel wall which yields Z D =0.3 m. Moreover, non-negligible divergence of light rays induces a dependence of the refraction amplitudes to the position of the displacement vector with respect to the optical axis. Hence, in order to mitigate this influence, the camera is positioned at Z A =2.5 m yielding a field of view (FOV) of 0.18 × 0.18 m 2 (Fig. 13, 2nd stage).
After the optical setup is configured, calibration images are acquired to verify the estimates of the magnification factor (M). In accordance to the measured M, the necessary parameters for the random dot and wavelet noise backgrounds are determined based on the conclusions drawn from the synthetic assessment (Fig. 13, 3rd stage). Hence, a random dot background pattern with particle density of N I =0.027 ppp is created which corresponds to 7 particles per interrogation window of 16× 16 pix 2 . The dots are generated as white squares over a black background. Then, the diffraction effects yielded Gaussian illumination distributions which provide the desired intensity gradients for the optical flow algorithm as well (Fig. 5, bottom-left). Moreover, four wavelet noise backgrounds are generated based on 3 different noise patterns with constant initial resolution level ( n I =2) while the final resolution level is varied between 10≤ n F ≤ 12 (Fig. 12). Since the camera resolution level cannot be matched exactly by the printers, the highest resolution that can be captured within the acquired images can only be estimated up to a certain range. Then, the corresponding range is aimed to be matched by the WN patterns without exceeding resolution of the schlieren features dictated by d Σ . Additionally, the practical printing limitations are taken into account for determining the pixel projection for each noise component as the influence of printing quality becomes more profound for WN patterns due to the increased range of gray color levels.
It is well documented for BOS that the requirement of focusing on the background not only yields diffraction effects to influence the measurement in close proximity of the objects, but also causes the flow structures of interest to be out of focus. As a consequence, various optical artifacts such as blur over shocks are induced (Sourgen et al. 2012;Nicolas et al. 2017). Hence, a lower aperture is always preferred to increase the depth of field. Yet, this is prohibited by Fig. 11 Experimental setup of the background-oriented schlieren system in the S-4 supersonic wind tunnel of VKI the corresponding loss of illumination intensity. It reduces the SNR of illumination intensities over the recorded background images which hinders the reconstruction procedure severely regardless of the algorithm of choice. Thus, in consideration of the steady nature of the major flow structures of interest and limitation of the high-resolution camera to reach the temporal range of the expected unsteadiness in a Mach 3.5 supersonic flow stream, the exposure time of the camera is set to 100 ms. Accordingly, 40 wind-on images for each model and background combination are recorded (Fig. 13, 4th stage) with an acquisition frequency of f aq =10Hz which yields a random uncertainty of u r = ∼0.02 pix with at 95% confidence level (Coleman et al. 1999). An aperture of f # = 22 is used to minimize the defocusing effects of blur and diffraction in the vicinity of the model surface (Table 2).

Data processing
As the reference and wind-on images are acquired, the deflection sensing stage utilizing XCOR and OF is initiated without any image pre-processing (Fig. 13, 6th stage). This also enabled the quantification of sensitivity of the setup to brightness variations between different instances of wind-on images induced by the modulations in the power output of the DC power supply. The XCOR is applied with a multi-pass approach where 3 passes with consecutive reduction of IW sizes (64 × 64, 32 × 32, 16 × 16 pix 2 ) are employed with 50% overlap which yields a resolution of 0.016 vectors/pixel. After the final displacement maps are obtained, a local outlier removal step based on a median filter of size 3 × 3 IW 2 is utilized. While deflection sensing using XCOR is only employed over the RD backgrounds, the OF algorithm is utilized for capturing the displacements over both the RD and WN backgrounds. Similar to the synthetic assessment, the hierarchical scheme is implemented for the experimental data analysis. The number of levels is determined based on the size of the ROI in order to minimize interpolations errors between different levels of the pyramid scheme. Moreover, the intermediate warping steps at each pyramid level are determined to ensure the convergence of the energy function with residuals below 10 −2 . Then, the median filtering heuristic is applied at each intermediate warping step with a filter size of 5 × 5 pix 2 ( Table 3).
As demonstrated throughout the synthetic assessments, the selection of regularization parameter plays a key role in accurately reconstructing the displacement fields caused by the density variations with optical flow (OF). However, due  to the absence of a ground truth to compare with in case of an experimental scenario, the determination of regularization parameter becomes a user defined process. It is generally based on the assessment of each result by the user to find a suitable balance between noisy reconstructions with low values which allows preservation of sharper gradients and the oversmoothing of displacement vectors at with high values. Thus, in order to establish a quantitative framework for the selection of the value for the current experimental investigation, four samples of displacement vector magnitude variations over the leading edge shockwaves of the ogive and wedge are extracted and averaged by aligning the peak values of each line which is considered to be the shock location (Fig. 14, top). Then, the displacement magnitudes obtained with various lambda values are compared and the regularization parameter ( ) corresponding to the highest peak pixel displacement in absence of upstream perturbations (which is known to be zero apriori in the experiments) is selected for each case (Fig. 14, bottom).
Once the displacement fields are obtained, constant shifts of displacement magnitudes are denoted between each acquired image induced by the vibrations and illumination variations on the background in comparison to the reference image. These shifts are observed to be introducing a relatively constant density gradient contribution across the flow field. Thus, in order to correct for these spurious background effects, the displacements reconstructed using each image processing method along the central axis of the wind tunnel are averaged and subtracted from the whole field of displacement vectors. Then, the optical configuration details are used to relate to the deflection angles (Eq.3) induced by the refractive index variations with the assumption of deflections being small ( ≈tan ).
Moreover, the relation between the refractive indices (n) and the density field ( ) is provided by the Gladstone-Dale relation (n=1+K ) where K refers to the Gladstone-Dale constant which has a weak dependency on the wavelength of light utilized to illuminate the background and for air it receives a value of 0.223×10 −3 m 3 /kg. For two-dimensional investigations the density gradients are assumed to be constant in the line of sight of the camera (S). Hence, the integration of the light ray trace path through the refractive index variations can be computed using Eq. (4).
Furthermore, since the object width (especially for the 10 o wedge) is not negligible in comparison to Z D , the width of the schlieren object in direction of the optical axis (W) is taken into account. Hence, using W and K, the density gradients are computed from the displacement fields after Fig. 13 Schematic of the workflow for the BOS measurement system which they are integrated via a Poisson relation with Neumann boundary conditions to obtain the density field with respect to a reference point within the ROI (Fig. 13, 7th  stage). Finally, the Dirichlet boundary condition at the reference point is computed using the isentropic relations for the expansion through the wind tunnel nozzle utilizing the monitored reservoir pressure and temperature, and the nozzle design Mach number of 3.5.

Results
Inıtial assessment of the experimental results is performed by first analyzing the reconstructed flow structures of the supersonic flow over the two models. A strong conical shockwave is attached at the nose of the ogive (Fig. 15, left). Although this shock is generated with a conical form at the tip, its angle reduces towards the downstream of the model. This is owing to the distinctive design of the von Karman ogive with a relaxing surface tangent which causes expansion waves to be generated along the surface of the model which interact with the conical shock to reduce its intensity thus decreasing its steepness (Dolling and Bogdonoff 1979).
On the other hand, the leading edge shock of the wedge is observed to be straight due to the planar structure of the model. Yet its strength and thickness is increased further away from the model (Fig. 15, right). This is caused by the weak compression waves generated due to the effective compression surface modulations induced by surface roughness.
Since the misalignment of these compression waves are not as severe as the expansion waves generated over the ogive, they do not introduce a noticeable difference in the shockwave angle. Yet, changing the post shock flow properties in the direction of increasing pressure and density yields a thicker and more intense density variation region as the compression waves merge with the shock. Downstream of the oblique shock, the supersonic flow proceeds towards the expansion shoulder of the wedge where it re-aligns with the incoming stream and a strong expansion fan is generated.
All three approaches are able to visualize the variation of the conical shock angle over the ogive. The increased resolution brought by OF enables a sharper identification of the shock location at the tip both with RD (Fig. 15, middleleft) and WN (Fig. 15, bottom-left) backgrounds compared to XCOR (Fig. 15, top-left). Nevertheless, the differences between the utilized reconstruction approaches become more severe when finer details in the weak expansion region over the ogive surface is investigated. The line of sight integration of the BOS method incorporates the compression induced by the surfaces of the conical shock apart from the symmetry plane. The non-planar contribution to the deviation of light rays overwhelms the weaker expansion waves reconstructed with XCOR which suffers from lower sensitivity and dominance of large pixel displacement values within the IW (Fig. 15, top-left). On the contrary, the capability of OF to capture lower intensity flow features enables the weak expansion waves to be reconstructed, which yields a compression region confined in the vicinity of the shock wave as expected (Fig. 15, middle-left & bottom-left). Moreover, the clarity of the flow features in the base region of the ogive (separated shear layer, the lip shock and the reattachment shock (van Gent et al. 2019) further denotes the increased resolution characteristics of OF against XCOR (Fig. 15,  left).
The improved resolution of small scale flow structures with OF-based deflection sensing persists for the wedge model. The weak compression waves generated over the compression ramp are captured with greater fidelity compared to the XCOR where only mild displacement magnitudes are reconstructed (Fig. 15, right). In addition to the resolution limit, the bandwidth limited characteristic of XCOR restricts the dynamic range of the ray displacements. Hence, in close proximity of expansion shoulder of the wedge, the displacement magnitudes reach beyond ∼ 8 pixels which is above the resolvable limit of 16× 16 pix 2 IWs with 50% overlap. Thus, in combination with the spatial smoothing effect, the displacements are severely underestimated (Fig. 15, middle). Nonetheless, in regions further away from the model where the density gradients are more outspread with lower amplitudes, the displacement vectors are reconstructed with similar accuracy using both XCOR-and OF-based approaches (Fig. 15, right). Furthermore, the vertical displacement magnitudes are extracted over selected lines (Y = −16mm, −22mm and −28mm) to provide a quantitative characterization of the reconstruction quality of various methods (Fig. 16). In the vicinity of the shock wave where the density gradients are confined to small regions, displacement magnitudes captured by XCOR are not only poorly localized to precisely pinpoint the shock location, but also reconstructed with significantly lower amplitudes compared to OF. Accordingly, WN backgrounds yield higher noise levels over the displacement maps which interfere with the sensitivity of the density gradients of interest. The increased noise of the reconstruction is related to the changing contrast range at different image pyramid levels (as shown in Fig. 12), while the contrast range of RD backgrounds does not vary when the acquired images are subsampled at different levels of resolution. Thus, even though the WN preserved its orthogonality by keeping a Gaussian distribution at every level and region of sampling, the degraded intensity of gradients reduces the confidence level of OF estimation. The loss of contrast is observed to be due to a combination of experimental factors including the quality of the paper on which the backgrounds are printed, the printer itself and the light source. Hence, employing a quadratic penalty definition, this situation shifts the influence on the energy function towards the brightness constancy residuals which causes the generation of spurious small-scale structures with relatively large amplitudes.
In order to quantify the physical characteristics of the reconstructed flow structures and the differences achieved with various methods, the physical alignment of the shock waves and expansion fans are computed by identifying locations of higher magnitudes of second order density gradients (Laguarda et al. 2021). For the ogive, a third-order polynomial is fitted over the identified shock location to capture the variations in shock angle (Fig. 17, left). For the wedge, the leading edge shock and the expansion fan are reconstructed by a linear regression (Fig. 17, right). First of all, the OFbased reconstructions revealed a closer agreement compared to the XCOR for capturing the conical shock angle at the tip of the ogive where the displacements are significantly localized (Table 4). Furthermore, the linearity of the shock owing to the straight contour of the wedge allowed all three approaches to capture the relevant shock angle accurately irrespective of the sensitivity of reconstructions in close proximity of the model. The slight deviation from the theory for all three methods is caused by a small modulation in the angle of attack of the wedge in the negative direction during the wind tunnel operation which caused the upper surface deflection angle to be greater than 10°. Then, comparing the angles of the first and last Mach waves of the PM expansion fan, the OF slightly outperforms the XCOR. Nonetheless, the angle of the last Mach wave is overestimated using all three approaches. This is due to the lower magnitude of the density gradient jump at the end of the PM expansion fan which is in theory discontinuous. While the resolution limits of the BOS method fail to clearly identify the exact termination point of the expansion fan, the reduced level of displacement magnitudes enabled XCOR to achieve similar accuracy characteristics to that of the OF (Table 4).
The differences between the XCOR and OF are further demonstrated as Eqs. (3) and (4) are utilized to compute the density gradients (Fig. 18, left). The oblique shockwave is clearly identified to be attached at the leading edge of the model with OF. Additionally, the capability to capture higher density gradients over the expansion fan denotes the increased dynamic range of displacement detection against XCOR. The displacement magnitudes are severely underestimated with XCOR in close proximity of the expansion shoulder observed in Fig. 15 (right) propagate through the density gradients. Furthermore, significant differences can also be noticed throughout the boundary layers not only between the XCOR and OF but also between the different backgrounds utilized for OF-based deflection sensing stage. Firstly, the diffraction effects in close proximity of the model degrade the illumination intensities which reduces the SNR of the correlation maps for the XCOR, while OF yields a more robust reconstruction of the thermal boundary layers with increased sensitivity. However, this increased robustness is observed to be associated with the RD backgrounds. Lower contrast over the background pattern of WN suffers from the reduced illumination strength which restricts the capabilities of OF to capture the density gradients in contact with the surface (Fig. 18, middle-left & bottom-left).
The procedure of density gradient integration is only applied to the wedge flow data to compute the density fields. Two main reasons are attributed for this choice. Firstly, BOS measurements enable line of sight integration of the refractive index variations to relate the displacements in the background pattern between the wind-on and reference images. Thus, a three-dimensional model such as the ogive violates the assumptions that allows the construction of Eq.4 since the density gradients along the direction of the optical axis cannot be assumed constant. Hence, computation of two-dimensional density gradients along the symmetry plane of the ogive can either be performed through tomographic reconstructions (Amjad et al. 2020) or by the use of mathematical relations, i.e., Abel inversion, under the strict assumption that the flow properties are axisymmetric (Tan et al. 2015).
The incapacity of BOS for capturing the density jump accurately across the shock wave is regardless of the deflection sensing approach used (Fig. 18, right). The physical thickness of the shockwave is in the order of 200 nm (Puckett and Stewart 1950;Robben and Talbot 1966), while the BOS system has a pixel projection of 90μm/pix. Hence, the density variation across the shockwave is not resolvable due to the optical limitations of the setup. Yet the intensity of the density increase is observed to be elevated further away from the model. This is attributed to the combined effect of the dissipation of the shockwave and the compression waves emanating from the ramp surface. Furthermore, higher density values at the pre-shock state are computed using XCOR as the expansion waves originating from the nozzle surface are captured with lower amplitudes compared to OF. As the integration proceeded downstream towards the expansion fan, the severity of underestimations with XCOR amplified especially in close proximity of the wedge due to the limited dynamic range. Nevertheless, OF yields lower density values at the post-expansion state without being overpowered by the compression region further away from the model. The increased accuracy using OF is also elaborated over the thermal boundary layers which contain lower density fluid due to the heat transfer takıng place between the relatively hot model and the cold supersonic flow. On the other hand, XCOR suffered from both lower resolution and deteriorated reliability on the formation of high correlation coefficients in vicinity of the boundary layers (Fig. 18, top-right).
Furthermore, the density profiles over three lines (Y = −17mm, −21mm and −25mm) are extracted to perform quantitative comparisons with the theoretical values of density variations across the shock wave and the expansion fan (Fig. 19). Firstly, the aforementioned underestimation of the density gradients are noted between XCOR and OF as the overall density variation is consistently lower with XCOR. Secondly, the standard deviation of density fields over the acquired image instances is observed to be significantly higher for the OF with WN backgrounds. This is associated to the increased vulnerability of WN backgrounds to the variations in illumination strength as well as to the random measurement and reconstruction noise. The lower contrast between different levels of the gray scale noise pattern, causes WN backgrounds to yield noisier reconstruction fields which is further amplified by the higher sensitivity of  the optical flow approach. Hence, this vulnerability propagates from the density gradients (Fig. 18, left) to the density fields (Fig. 18, right). As a result, large modulations in the flow field are captured between different measurement samples which yields considerably higher standard deviation levels. Moreover, for comparisons between the flow features and the theoretical values, the variation of density over the selected regions are normalized with the pre-shock ( ∞ ) and pre-expansion ( 1 ) flow states for the shockwave and the expansion fan respectively. Then, oblique shock relations (Anderson 2001) are employed to compute the density jump across the shock and represented with respect to the preshock density value ( ∕ ∞ ) ( Fig. 19A.1, B.1 and C.1). Even though the OF-based approaches are able to capture a higher increase of density, the results are observed to be deviated up to 70% of the theoretical values for all three methods. This is mainly governed by the aforementioned resolution problem of the density jump across the shockwave which is discontinuous with respect to the resolution of the BOS setup.
Finally, a similar comparison is performed for the density variations reconstructed over the expansion fan ( Fig. 19A.2, B.2 and C.2) with respect to the density value at the pre-expansion state ( ∕ 1 ) utilizing the Prandtl-Meyer theory (Anderson 2001). Closer to the model, a good agreement of the density variations using OF in comparison to the theoretical values is obtained while measurements with XCOR are significantly off from the theory (Fig. 19top-A.2). This is a clear quantification of the inaccuracies introduced by the combined effect of the bandwidth filtering and lower resolution of the XCOR method in comparison to the multi-scale reconstruction algorithm of OF-based deflection sensing. Yet, the reduction of effective resolution of the schlieren object in order to have the background in focus causes density gradients to be dissipated even for OF. Additionally, the smoothness constraints associated with OF prevents displacements to be discontinuous which is a severe limitation for capturing the edges of the theoretical expansion fan. Moreover, further away from the model the differences between the XCOR and OF reconstructions are reduced owing to the mitigated localization of the expansion fan and the corresponding reduction in the amplitude of the density gradients. Nonetheless, deviations from the theoretical values are increased for all three methods due to the reflected shock near the wind tunnel walls, causing the density in the post-expansion region to rise above the theoretical predictions.

Conclusions
Simplicity of application and minimal requirements for its experimental configuration makes background-oriented schlieren (BOS) a cost-effective solution for quantitative visualization of density varying flows. However, degraded resolution and sensitivity characteristics of BOS outweigh its benefits in comparison to parallel light schlieren Fig. 18 Time averaged density gradients and density fields for supersonic flow past 10° wedge reconstructed with XCOR (top), OF with random dot (middle) and wavelet noise (bottom) backgrounds systems. In this regard, optical flow (OF) provides a unique opportunity to recover high fidelity density field variations at the deflection detection stage of BOS. Thus, in order to have a full comprehension of its capabilities over the conventional cross-correlation (XCOR) approach, a comparative study including synthetic and experimental measurements of supersonic flows over various models is conducted. The OF-based reconstructions revealed superior characteristics of resolution and sensitivity against XCOR which allowed higher accuracy reconstructions of displacement vectors. The hierarchical schemes enabled OF to resolve the full spatial-frequency spectrum of the displacements contained within the image plane where the bandwidth limited nature of XCOR suffered not only from smoothing due to spatial averaging, but also from a limited dynamic range. The influence of varying background patterns and OF configurations is also investigated with wavelet noise that allows the spatial frequency resolution specification of the image pyramids for OF to be complemented. Even though synthetic assessments performed in a controlled environment attributed the greatest accuracy to OF application with procedural background patterns, unavoidable experimental uncertainties of illumination and printing artifacts make random dot patterns to be the most robust option as multiscale patterns are vulnerable against the noise. Availability of data and materials All data that support the findings of this study are included within the article and accessible via the authors upon request.

Conflict of interests
The authors have no competing interests as defined by Springer, or other interests that might be perceived to influence the results and/or discussion reported in this paper.

Ethical approval Not applicable.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.