# Design rules for Background Oriented Schlieren experiments with least-squares based displacement calculation

Bertrand Mercier, Tom Lacassagne

► To cite this version:

Bertrand Mercier, Tom Lacassagne. Design rules for Background Oriented Schlieren experiments with least-squares based displacement calculation. Flow Measurement and Instrumentation, 2023, 92, 10.1016/j.flowmeasinst.2023.102395. hal-04249140

HAL Id: hal-04249140

https://hal.science/hal-04249140v1

Submitted on 19 Oct 2023

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
<figure><img src="image_2.png" /></figure>

# Design rules for Background Oriented Schlieren experiments with least-squares based displacement calculation

Bertrand Mercier, Tom Lacassagne*

IMT Nord Europe, Institut Mines-Télécom, Univ. Lille, Centre for Energy and Environment, Lille, F-59000, France

## ARTICLE INFO

MSC:

0000

1111

Keywords:

Background oriented Schlieren

Design rules

Optical flow

Matlab

## ABSTRACT

Background Oriented Schlieren (BOS) is a non-intrusive optical method used to visualize density variations in fluids. A camera is used to measure the distortion of a background image displaying a random pattern due to the refraction of the light induced by the local gradients in the refractive index of the under study. In contrast to conventional schlieren, for which direct observation is possible, the images obtained with the BOS technique have to be processed. In this contribution, a computationally affordable approach based on linear least-squares minimization and adapted from optical flow techniques is presented in detail and tested in the context of BOS with the Matlab GUI application comBOS. A series of synthetic images are first analyzed to devise design rules for the definition of an optimal random Gaussian dot background pattern. Tests are performed for comBOS, the commercial BOS module included in DaVis 10, and the free PIV tool PIVlab. It is found that the diameter of the dot images on the camera sensor should be at least 3 pixels, but can be larger without affecting the accuracy of the measurement with comBOS and DaVis 10, while PIVlab severely degrades the quality of the results for dots larger than 5 pixels. The contrast and the density of dots must also be sufficient such that a minimum of 75 % of the surface of the background is shaded. If these conditions are met, the maximum systematic error in the measured distortions is 0.006 pixel for comBOS and 0.001 pixel with DaVis 10 when local gradients are small. Similar errors are found with PIVlab for dots smaller than 5 pixels in diameter, but they increase rapidly for larger dots. The random error is typically smaller than 0.004 pixel in optimal conditions, but it increases significantly for images affected by a random noise or local gradients of the distortion. Finally, it is concluded with a realistic test case that both comBOS and DaVis 10 are well suited to capture the main features of a complex compressible flow with Schlieren imaging. Thus, the main advantage of comBOS is its affordability and robustness for researchers and educators willing to implement BOS, for example in educational contexts.

## 1. Introduction

Schlieren techniques have been used under a wide variety of conditions and include many forms of implementation to describe the structure of transparent flows with inhomogeneous refractive index. The present study focuses on the Background Oriented Schlieren (BOS) method which is the most recent Schlieren method and that continues to be a source of innovation for experimentalists. The concept of this method indeed first appeared in a peer-reviewed journal in 2000 under the name "synthetic Schlieren" by Dalziel et al. [1] and as "Background Oriented Schlieren" by Raffel et al. [2] the same year. The latter is, today, the most common name in the fluid dynamics community to refer to this technique. BOS has many advantages in comparison with other Schlieren implementations, the most important being the simplicity of the setup which relies only on the use of a camera that captures images of a background pattern and the pattern

itself. The output of a BOS analysis is a measure of the deformation of the pattern as seen by the camera, which is related to the deflection of the light induced by the gradients of refractive index within the medium between the camera and the background. The random pattern is often displayed or projected onto a screen, but it can also be a natural random image such as a landscape [3], a plugged field [4], or even the fluctuation in intensity of the band-pass filtered light at the surface of the sun [5] which respectively raised the possibility of observing the flow around an airborne helicopter and fighter jet. The ability to observe large fields is therefore the second main advantage of BOS. The last major advantage is the possibility to reconstruct three-dimensional density fields by using multiple cameras and tomographic techniques as demonstrated in this non exhaustive selection from the long list of studies since the early 2000's: [6–11]. The drawback of BOS in comparison with conventional Schlieren techniques is that direct

<footnote>* Corresponding author.
E-mail address: tom.lacassagne@imt-nord-europe.fr (T. Lacassagne).</footnote>
observation of the measurement is not possible because it requires a post-processing step that takes a significant amount of computing time. There are various methods to perform this post-processing, some of them are specifically optimized for efficiency with a high level of accuracy such as the use of optical flow based methods [12–14] or other specific methods [11,15,16], but most of the studies using BOS in the literature rely on the classical cross-correlation and sub-pixel interpolation that are included in PIV tools [7,17–21]. Cross-correlation is probably the dominant approach because many PIV tools are available either in open-source [22–25], in-house, or commercial packages, and because the quality of the results obtained with this method is satisfactory. The counterpart is that these algorithms are optimized to accurately determine displacements of several pixels between images with a large amount of distortion. Therefore, they are not expected to perform optimally for displacements of 0.1 pixel and a low level of distortion, although these are the generally observed in pairs of BOS images in air in the absence of shockwaves, thus in most cases. The robustness of PIV tools can still be beneficial when these conditions are not met, but this may require computationally expensive processing.

The simplicity of the arrangement of BOS apparatus and the potential outcomes of this flow visualization technique for pedagogical application stressed the authors’ motivation to develop an open-source tool that can be operated by students on their personal laptop, thus requiring an optimized algorithm to measure accurately the small displacements with a limited computational cost. This need led the authors to develop a Matlab application called comBOS which has already demonstrated its ability to provide results of quality similar to commercial packages for BOS processing with smartphone cameras [26].

The objectives of this article are two folds. First, it intends to present the post-processing method implemented in comBOS. Second, and more importantly, this study is also motivated by the need for providing design rules and guidelines to set up a BOS experiment with a least-squares based method, which are not expected to be the same as with a cross-correlation method as in PIV. There are indeed recurring questions in the literature about how to choose the parameters to generate a random pattern, and the effect of these parameters on the overall uncertainty of the BOS measurements. These two concerns are in fact not answered by the few former studies [7,17–19,27,28] addressing the capabilities of BOS measurements. Only Hartberger [29] focused on the properties of the pattern, but with a limited number of adjustable parameters.

A brief description of the BOS principle is presented first, followed by a description of the post-processing method implemented in comBOS. Design rules are then devised from a series of benchmarks to design an optimal pattern with comBOS and the commercial BOS package included in DaVis 10 by LaVision. It should be noted that although there is no public description of the algorithm used by LaVision for the BOS module, it is thought to be an algorithm specifically designed for BOS and not a cross-correlation based algorithm. Based on these rules, the errors induced by both software are evaluated for realistic constraints including large displacements and gradients. The errors are also compared with displacements calculated by PIVlab using cross-correlations. A summary and concluding remarks complete this contribution.

## 2. Principles of Background Oriented Schlieren visualization

The BOS technique, like other Schlieren techniques, involves the deflection of light rays caused by the inhomogeneities in the refractive index $n$ of the medium being analyzed. The light path is governed by the following set of differential equations [30] written in the system of coordinates illustrated in Fig. 1:

$$ \begin{cases} \frac{d^2 x}{dz^2} = \frac{1}{n} \left( 1 + \left( \frac{dx}{dz} \right)^2 + \left( \frac{dy}{dz} \right)^2 \right) \left( \frac{\partial n}{\partial x} - \frac{dx}{dz} \frac{\partial n}{\partial z} \right) \\ \frac{d^2 y}{dz^2} = \frac{1}{n} \left( 1 + \left( \frac{dx}{dz} \right)^2 + \left( \frac{dy}{dz} \right)^2 \right) \left( \frac{\partial n}{\partial y} - \frac{dy}{dz} \frac{\partial n}{\partial z} \right). \end{cases} \quad (1) $$

These equations describe the local curvature of the light rays induced by the gradients of the refractive index. This problem is generally simplified by assuming paraxial conditions, hence small values for $dx/dz$ and $dy/dz$ with respect to 1, and by assuming that the gradients of refractive index are similar in all directions. This leads to the following new set of equations

$$ \left\{ \begin{aligned} \frac{d^2 x}{dz^2} &= \frac{1}{n} \frac{\partial n}{\partial x} \\ \frac{d^2 y}{dz^2} &= \frac{1}{n} \frac{\partial n}{\partial y}. \end{aligned} \right. \quad (2) $$

The validity of these simplified equations may be questioned at the edges of BOS images obtained with systems equipped with lenses of short to medium focal length because the light rays experience moderate to large angles with respect to the optical axis $\bar{z}$. However, the BOS is often used solely as a qualitative diagnostic technique, and as long as 3D reconstruction is not intended to be done, this set of equations is generally considered as satisfactory for describing the BOS principle in most experimental conditions. In addition, this simplification allows a direct and straightforward integration of the curvature along the light rays from $z = 0$ to the background at $z = L_B$ which results in the deflection angles $\epsilon_x(z)$ and $\epsilon_y(z)$ with respect to the undisturbed ray in the $(\bar{x}, \bar{z})$ and $(\bar{y}, \bar{z})$ planes with

$$ \left\{ \begin{aligned} \epsilon_x &\approx \tan(\epsilon_x) = \frac{dx}{dz} = \int_0^{L_B} \frac{1}{n} \frac{\partial n}{\partial x} dz \\ \epsilon_y &\approx \tan(\epsilon_y) = \frac{dy}{dz} = \int_0^{L_B} \frac{1}{n} \frac{\partial n}{\partial y} dz. \end{aligned} \right. \quad (3) $$

Noting that if BOS is employed to analyze fluids at low pressure relative to their critical pressure, $n$ is expected to be close to 1, and an average refractive index $n_0$ can be used to further simplify Eq. (3). Furthermore, the position of the point at which a given light ray reaches the background is shifted from the interception point of the undisturbed ray by $\Psi_x$ and $\Psi_y$ in respectively the $\bar{x}$ and $\bar{y}$ directions, respectively, according to the notations in Fig. 1. They are given by

$$ \begin{aligned} \Psi_x &= \frac{1}{n_0} \int_0^{L_B} \frac{\partial n}{\partial x} z_B dz \\ \Psi_y &= \frac{1}{n_0} \int_0^{L_B} \frac{\partial n}{\partial y} z_B dz \end{aligned} \quad (4) $$

where $z_B = L_B - z$ is the remaining distance to the background.

The BOS technique ultimately consists in measuring $\Psi_x$ and $\Psi_y$ which are related to the magnitude of the gradients of the refractive index between the lens and the background. However, unlike conventional Z-type Schlieren apparatus, the results also depend on the position in $z$ at which the gradients occur. The displacements measured from a phenomenon occurring close to the lens of the camera will be amplified in comparison to the same phenomenon happening near the background. This is a useful way to adjust the sensitivity of the system, but it makes it challenging to interpret the results if the axial extent of the phenomenon is significant compared to $L_B$.

## 3. Least-squares based algorithm for the estimation of the displacements

Images acquired for BOS visualization are often processed using correlation based algorithms developed for particle image velocimetry (PIV) applications. These algorithms are optimized to determine the displacement of groups of particles advected by a flow. Particles typically move by several pixels over the period of time between the two images, therefore the non-uniformities of the flow velocity introduce a significant amount of distortion in the pattern represented by the particle group. These displacements are typically at least an order of magnitude larger than with BOS. The relevance of algorithms based on linear least-squares fitting such as the one implemented in the toolbox
<figure><img src="image_3.png" /><figcaption>Fig. 1. Schematic of a BOS imaging system.</figcaption></figure>

comBOS for BOS lies in the fact that they are optimized to perform accurate and time-efficient estimations of the displacements in BOS images due to the low level of distortion across a pair of images. Once projected onto the camera sensor, the displacements $\Psi_x$ and $\Psi_y$ of the background are in fact smaller than a pixel in most applications.

As is the case with standard PIV algorithms, the present method requires that the image be first divided into subsets. A displacement is then measured by estimating the translation that should be applied to a subset of the deformed image to match the subset at the same position in a reference image. Here, the subsets are squares of arbitrary size.

In details, and following the notations defined in Fig. 2, let $f(x_i, y_i)$ and $g(x_i, y_i)$ be the grey levels of the pixel $i$ at position $(x_i, y_i)$ in the reference and the deformed subsets, respectively, the horizontal and vertical displacements $\Delta x$ and $\Delta y$ are defined such that

$$af(x_i, y_i) + b = g(x_i + \Delta x, y_i + \Delta y). \quad (5)$$

Coefficients $a$ and $b$ account for overall changes in the luminosity and the contrast between the reference and the deformed images. $\Delta x$ and $\Delta y$ are further decomposed into integer pixel displacements $D_{x,y}$ and sub-pixel displacements $\delta_{x,y}$ following

$$\Delta x = D_x + \delta_x \quad (6)$$

$$\Delta y = D_y + \delta_y \quad (7)$$

### 3.1. Integer-pixel displacements

If maximum displacement in a BOS image is smaller than one pixel, $D_x$ and $D_y$ are both equal to zero and do not need to be considered. However if this condition is not met, $D_x$ and $D_y$ must be computed prior to $\delta_x$ and $\delta_y$ which is enabled by an option in the processing settings in comBOS.

When needed, this first step is performed through the computation of the cross-correlation function $C(k, l)$ between $\tilde{f}(x_i, y_i)$ and $\tilde{g}(x_i, y_i)$ where $\tilde{f} = f - \langle f \rangle$ and $\tilde{g} = g - \langle g \rangle$, $\langle \cdot \rangle$ being the spatial average. $D_x$ and $D_y$ are the values of $k$ and $l$ where $R(k, l)$ is maximum. Note that if the peak of the correlation was estimated with a sub-pixel resolution, this step would be similar to a standard PIV processing. Also, given the fact that the displacements are expected to remain small, and in order to keep the computation time as low as possible, the cross-correlation function is computed on interrogation windows centered in the subset of interest and of size $M$ equal to half the size of the subset, thus considering a quarter of the pixels, with

$$R(k, l) = \sum_{i=0}^{M-1} \sum_{i=0}^{M-1} \tilde{f}(x_i, y_i) \tilde{g}(x_i - k, y_i - l). \quad (8)$$

### 3.2. Sub-pixel displacements

Sub-pixel displacements are computed using a technique adapted from Pan et al. [31] and similar to the Lucas–Kanade method [32] which has been recognized as a suitable algorithm for BOS [12], but with additional degrees of freedom for the variations of luminosity and contrast of the images. The method consists in estimating $\delta x$ and $\delta y$. To this end, the gray level in this deformed subset follows the function $h(x_i, y_i)$, $x_i$ and $y_i$ still being the pixel position in the subset reference, not in the full image.

With this new subset, Eq. (5) becomes

$$af(x_i, y_j) + b = h(x_i + \delta x, y_j + \delta y) \quad (9)$$

The Taylor expansion at order one of Eq. (9) $F_{i,j}$ is written for each pixel at position $x_i$ and $y_j$, which leads to

$$F_{i,j}(p) = h_{i,j} + \frac{\partial h_{i,j}}{\partial x} \delta x + \frac{\partial h_{i,j}}{\partial y} \delta y - af_{i,j} - b \approx 0. \quad (10)$$

The derivatives $\frac{\partial h_{i,j}}{\partial x}$ and $\frac{\partial h_{i,j}}{\partial y}$ are calculated with the centered finite difference scheme

$$\frac{\partial h_{i,j}}{\partial x} \approx \frac{1}{12} h(x_{i-2}, y_j) - \frac{2}{3} h(x_{i-1}, y_j) + \frac{2}{3} h(x_{i+1}, y_j) - \frac{1}{12} h(x_{i+2}, y_j) \quad (11)$$

$$\frac{\partial h_{i,j}}{\partial y} \approx \frac{1}{12} h(x_i, y_{j-2}) - \frac{2}{3} h(x_i, y_{j-1}) + \frac{2}{3} h(x_i, y_{j+1}) - \frac{1}{12} h(x_i, y_{j+2}). \quad (12)$$

No specific treatment is implemented to compute the derivative for the first two layers of pixels on the edges of the subset, therefore these pixels are not considered.

Eq. (10) can be written as a linear system of equations for the $n$ pixels considered in the subset.

$$Xp = b \quad (13)$$

where

$$X = \begin{pmatrix} \frac{\partial h_1}{\partial x} & \frac{\partial h_1}{\partial y} & -f_1 & -1 \\ \frac{\partial h_2}{\partial x} & \frac{\partial h_2}{\partial y} & -f_2 & -1 \\ \vdots & \vdots & \vdots & \vdots \\ \frac{\partial h_n}{\partial x} & \frac{\partial h_n}{\partial y} & -f_n & -1 \end{pmatrix}, \quad p = \begin{pmatrix} \delta x \\ \delta y \\ a \\ b \end{pmatrix}, \quad b = \begin{pmatrix} -h_1 \\ -h_2 \\ \vdots \\ -h_n \end{pmatrix}. \quad (14)$$

If $n \ge 4$, the vector $p$ can be determined by a linear least-squares method:

$$p = (X^T X)^{-1} X^T b \quad (15)$$

The use of this single pass linear method is however restricted to small displacements. Indeed, a 5% error is observed on the computed values $\delta x$ and $\delta y$ for displacements of 0.25 pixels. This bias is significantly reduced by repeating the same sub-pixel displacement calculation method on the subsets described by $f(x_i, y_i)$ and by $\hat{h}(x_i, y_i)$. $\hat{h}(x_i, y_i)$ is computed by mapping $h(x_i + \delta x, y_i + \delta y)$ onto the regular grid formed by $(x_i, y_i)$ using bilinear interpolation.

### 3.3. Least-squares sensitivity

The sensitivity of the unknown vector $p$ to the noise in the observed value of $b$ is directly related to the conditioning of the matrix $X$. Let $\|\delta b\|/\|b\|$ be the relative uncertainty of $b$, then, the relative error $\|\delta p\|/\|p\|$ for $p$ is, according to [33], upper bounded as

$$\frac{\|\delta p\|}{\|p\|} \le \kappa(X) \frac{\|\delta b\|}{\|b\|} \quad (16)$$

where $\kappa(X)$ is the conditioning number of $X$ that is equal to the ratio between the largest and the smallest singular values of $X$.

A small value of $\kappa(X)$ mitigates the propagation of the uncertainty in $b$ to $p$. Conversely, if the value of $\kappa(X)$ is infinite, which happens when the smallest singular value $X$ is equal to zero, then $p$ will not
<figure><img src="image_4.png" /><figcaption>Fig. 2. Schematic of the coordinate systems in the reference and deformed subsets.</figcaption></figure>

be consistent. In practice, a null singular value can be found in X if the considered subset of the BOS images contains patches of strictly constant gray levels. In that case, the derivatives of $h(x_i, y_i)$ are equal to zero, and the values of $f(x_i, y_i)$ are constant, so the corresponding columns in X are identical. This issue could be mitigated by regularizing Eq. (15), but unless the camera sensor is overexposed, and thus the gray level of the pixels involved is equal to the maximum output value of the camera, the noise of the sensor is sufficient to prevent $\kappa(X)$ from being infinite. However, it will remain high and the estimated displacement will be affected by a significant amount of noise. Regularizing would therefore make $X^T X$ always invertible, but it would not guarantee that the calculated displacements are consistent.

This brief analysis shows that it seems preferable to use speckles that feature a high density of gray level fluctuations with methods based on least-squares methods such as the one implemented in comBOS. This is also consistent with analysis proposed with PIV algorithm [28,29]. However, unlike the PIV algorithms, the present method does not rely on a subpixel interpolation, the quality which is affected by the shape of the peak in the cross-correlation function. Therefore, the constraints on the speckles are not expected to be the same with least-squares-based methods or cross-correlations. The following section aims to provide a quantitative set of requirements to optimize the speckle.

## 4. Sensitivity of comBOS and DaVis 10 algorithms to the speckle properties

A series of benchmarks are proposed in this section to evaluate the performance of comBOS and compare it to the BOS module of the commercial software DaVis 10. Comparisons with the free PIV tool PIVlab are also included as it is a commonly used tool for PIV and BOS. The result of this benchmarking is a list of design rules that can be used to generate optimal speckles.

### 4.1. Generation of synthetic images

The benchmarks consist of pairs of synthetic speckles, one being the reference, and the second being distorted by a known function. Speckles made of a random pattern of dots are considered because they are commonly used in experiments and can be generated with the speckle generator included in comBOS.

The random dot pattern for the reference speckle is generated as follows:

Step 1: Set the mean distance in pixels $\bar{d}_{\text{dot}}$ between the dots and initialize a regular grid with this distance as the grid spacing, and the size of the desired speckle.

Step 2: Alter the coordinates of each node of the grid by a random value, here bounded within $\pm 1/3 \bar{d}_{\text{dot}}$.

Step 3: Set the full width at half maximum of the dots $\phi_{\text{dot}}$ and add to the speckle image a black Gaussian patch centered on each node of the modified grid.

<caption>Table 1</caption>

Average (< . >) and rms (< .^2 >^{1/2}) horizontal and vertical displacements measured with DaVis 10 for four different methods to generate the deformed speckle.

<table><thead><tr><th></th><th>&lt; Δx &gt;</th><th>&lt; Δy &gt;</th><th>&lt; Δx^2 &gt;^{1/2}</th><th>&lt; Δy^2 &gt;^{1/2}</th></tr></thead><tbody><tr><td>Analytic</td><td>0.261</td><td>-8.99e-6</td><td>2.18e-4</td><td>3.05e-4</td></tr><tr><td>Linear Interp.</td><td>0.250</td><td>-4.76e-5</td><td>3.19e-4</td><td>4.63e-4</td></tr><tr><td>Cubic Interp.</td><td>0.250</td><td>-1.81e-5</td><td>2.15e-4</td><td>3.47e-4</td></tr><tr><td>Spline Interp.</td><td>0.262</td><td>-1.44e-5</td><td>2.21e-4</td><td>3.17e-4</td></tr></tbody></table>

Step 4: Trim the gray levels to ensure that they are all contained within 0 and 255, even in the case of overlapping dots, and save the image as a bitmap file.

This method ensures that the density of the speckle is globally homogeneous and avoids the generation of large areas with too few or too many dots that can be encountered if the coordinates of the dots were allowed to be completely randomized instead of constrained to a grid.

Four methods were tested to deform the reference speckle with an imposed displacement field. The first one consists in generating a new speckle following the same procedure as for the reference one, but by adding a local additional shift of known value to the positions of the nodes, hence to the dots. This will be called the analytic method. The three other methods consist in interpolating the original speckle on a grid modified in accordance with the desired deformation. Linear, cubic and spline interpolations are tested. For the comparison of the methods $\bar{d}_{\text{dot}} = 7$ pixel, $\phi_{\text{dot}} = 3.3$ pixel, and the deformed image is shifted by a uniform displacement of 0.25 pixel along x and 0 pixel along y. The average and standard deviation of the displacements computed with DaVis 10 from 17 × 17 pixel subsets are summarized for all the four methods in the Table 1. The use of the analytic method to shift the dots is found to introduce a significant error, as can be seen from the average value of $\Delta x$. This error appears because the gray levels of the pixels were calculated at discrete locations using a two-dimensional Gaussian function, when it should have been the average of a continuous Gaussian surface over the pixel area. Spline interpolation does not perform well either. The best results both in terms of both accuracy and noise level are obtained with the cubic interpolation, this method is therefore chosen to generate deformed speckles in the following benchmarks.

### 4.2. Effect of dot size and density

The first question to be addressed when designing a speckle for a BOS experiment is the size of the dot image projected onto the camera sensor. This is determined by the geometry of the experiment, the lens of the camera, the size of the pixels on the camera sensor, and the size of the dots on the speckle. This last parameter is often the easiest to control.

The effect of this parameter on the accuracy of the displacement calculation is tested by generating a speckle of size $2500 \times 2500$ pixel$^2$. The dot diameter $\phi_{\text{dot}}$ linearly increases from 1 pixel on the left of the
speckle to 10 pixels on the right and $d_{\text{dot}}$ is kept constant at 16 pixels. In addition, a constant shift $\Delta x = \delta_x = 0.25$ pixel is applied to the deformed image using the cubic interpolation methods as discussed previously. A sample of the speckle in presented in Fig. 3(a).

The images are then processed with subsets of 51 × 51 pixel². The results on the accuracy of the displacement calculation and the resulting noise are presented in Fig. 3(b) and (c). Mean displacements $\langle \delta x \rangle$, which indicate the accuracy, and rms displacements $\langle \delta x^2 \rangle^{1/2}$, indicating the noise in the estimation of the displacements, are calculated for each horizontal position, hence dot diameter, from the 2500/51 ≈ 49 vertically aligned vectors. A large value for $\langle \delta x^2 \rangle^{1/2}$ means that the algorithm is not stable for that dot diameter, so that the resulting displacements are affected by a significant amount of noise.

Both comBOS and DaVis 10 behave similarly with a noise level that increases drastically for dots smaller than 3 pixels (see Fig. 3(c)). For larger dots, it stays at a low level below $2 \times 10^{-3}$ pixel with a tendency to decrease slightly for increasing $\phi_{\text{dot}}$. Similarly, accuracy remains within $\pm 5 \times 10^{-4}$ pixel for $\phi_{\text{dot}}$ larger than 3 and is altered for smaller dots (see Fig. 3(b)). This observation is consistent with the discussion about the sensitivity of the method in Section 3.3; in fact, a large part of the subset is uniformly white for small dots, therefore, the X matrix is ill-conditioned.

This ability to handle dots with a wide range of diameters is in contrast to the 3 pixel diameter particle images typically required for PIV. The difference in sensitivity to the size of the dots of the two techniques, namely least-squares minimization and conventional cross-correlation with sub-pixel interpolation, is tested. Results are presented in Fig. 3(d) and (e) for a variety of configurations that involve comBOS, PIVlab, and additionally for this first benchmark, the PIV module of the commercial software Davis 8.4, as the laboratory does not own the PIV license for DaVis 10.

First, it appears that all tested configurations using a PIV algorithm produce results with a larger average error than the least-squares method implemented in comBOS when applied with two passes according to the procedure described at the end of Section 3.2. This error is also found to be highly dependent on the dot size, although increasing the number of passes seems to mitigate this problem with PIVlab. The noise in the PIVlab's results follow the same trend as the average error, reaching values comparable to comBOS for the optimal dot diameter and if the “robust correlation” option is set to “high”.

Davis 8.4 performs better than comBOS for dots less than 8 pixels in diameter, but results are similar to comBOS when the second pass is disabled. Fig. 3(d) and (e) demonstrate well the issues that may arise when using advanced algorithms with many parameters in conditions that would not require such features. In fact, the BOS modules of Davis 10 and comBOS require only the dimension of the window, and the overlap as input parameters, with performances that beat PIV algorithms in most cases. In a realistic situation, the user does not know the true displacement, and since the results are not stable with respect to the parameters, one may increase the number of passes in the hope of obtaining more accurate results. This would certainly be a wise procedure, considering the comparison between 2 passes and 3 passes in PIVlab, but it would lead to a significant increase in processing time. These processing times are summarized in Table 2 with interrogation windows of 48 × 48 pixels and 50% overlap for all configurations. In each case, the indicated duration is the average processing time per image of a sequence of 20 images. The duration are all in the order of one second, with a clear effect of processing complexity. If the diameter of the dots is between three and four pixels, the best free alternative to comBOS is PIVlab with 2 passes and the “correlation robustness” set to “high”. It indeed offers similar computation time and noise figure as the 2 passes processing of comBOS, but with slightly lower accuracy. This is the configuration that will be used in the following when it is relevant to compare least-squares based methods with cross-correlation and sub-pixel interpolation.

<table><caption>**Table 2**</caption><caption>Average computing time for the configurations considered in Section 4.2 with 48 × 48 pixel² subsets and 50% overlap. Processing are performed on a laptop equipped with an I7-12850HX cpu, and Matlab version is 2021b. For PIVlab, “robust” means that the “robust correlation” parameter is set to “High”, it is set to “Standard” otherwise. Unless a specific configuration is mentioned, the configuration in bold is used for all the test cases discussed in this article.</caption><tbody><tr><th colspan="2" scope="colgroup">comBOS V.1.1:</th></tr><tr><th scope="row">Single pass</th><td>Time per image<br>0.95 s</td></tr><tr><th colspan="2" scope="colgroup">Standard</th></tr><tr><th scope="row">Standard + large Δx</th><td>1.3 s<br>2.1 s</td></tr><tr><th colspan="2" scope="colgroup">Davis V.10.2 (BOS):</th></tr><tr><th scope="row">Standard</th><td>0.46 s</td></tr><tr><th colspan="2" scope="colgroup">PIVlab V.2.61:</th></tr><tr><th scope="row">2 passes</th><td>0.59 s</td></tr><tr><th scope="row">2 passes, robust</th><td>1.6 s</td></tr><tr><th scope="row">3 passes, robust</th><td>2.0 s</td></tr><tr><th colspan="2" scope="colgroup">Davis V.8.4 (PIV):</th></tr><tr><th scope="row">2 passes</th><td>1.0 s</td></tr></tbody></table>

The effect of the spatial extent of the area shaded by the dots on the accuracy of the displacement estimation is now evaluated. To this end, a series of new speckles with different dot diameters are generated with the dot distance $d_{\text{dot}}$ decreasing linearly from the left to the right as shown in the sample in Fig. 4(a).

In these speckles, whose original size is 4500 × 1200 pixel², a 25 × 25 pixel² large subset contains in average 3 dots if it is on the left edge and 66 dots if it is on the right edge. The number of dots per subset is referred to as the dot density in the following. By repeating the same analysis with identical dot densities but different dot diameters, it is possible to determine whether the number of dots is important or whether it is the area covered by the dot that is important. In fact, although it is reasonable to expect that increasing the number of dots in a subset should lead to a decrease in noise, it has just been shown that the diameter of the dots also affects the noise figures because of the varying area of uniform gray levels. These two effects are therefore compared on this series of speckles with dots in the same positions but different diameters. Two parameters are used in this analysis. First, the density, which follows the definition above. Second, the uniformity ratio, that is equal to the ratio of the number of pixels where both the horizontal and vertical derivatives of the gray level as defined in Eqs. (11) and (12) are equal to zero to the total number of pixels in the subset. A uniformity ratio of 1 corresponds to a uniform image whereas a dense speckle has a uniformity ratio of 0.

The analysis is conducted for $\phi_{\text{dot}} = 2.2$ pixel, $\phi_{\text{dot}} = 3.3$ pixel, and $\phi_{\text{dot}} = 5.0$ pixel with 25 pixels large subsets and a $\delta_x = 0.25$ pixel horizontal shift imposed to the deformed speckle. The uniformity ratio varies from 0 to 0.65 for $\phi_{\text{dot}} = 2.2$ pixel, 0 to 0.55 for $\phi_{\text{dot}} = 3.3$ pixel, and 0 to 0.49 for $\phi_{\text{dot}} = 5.0$ pixel. The reliability of the displacement estimation is assessed via the root mean square of the horizontal displacements $\langle \delta x^2 \rangle^{1/2}$ in subsets of dot density being within a given range. For instance $\langle \delta x^2 \rangle^{1/2}$ is computed from the displacements of all the subsets containing between 3 and 4 dots in the least dense region, or for subsets containing between 55 and 66 dots in the densest region.

Fig. 4(b) shows the noise level of the displacements computed with comBOS for the three dot sizes against the uniformity ratio. For uniformity ratio less than 0.25, all the three results collapse into a single curve that describes a noise which increases with the uniformity. This is the expected behavior from the prior discussion in Section 3.3 if the stability of the displacement calculation was driven by the uniformity ratio alone. Nevertheless, the results tend to diverge and become dependent upon the number of dots for more uniform speckles. For example, for the same uniformity ratio of 0.5, the noise level is significantly higher if there are dots with a diameter of $\phi_{\text{dot}} = 5$ pixel, hence few large dots, than for dots with a diameter of $\phi_{\text{dot}} = 2.2$ pixel,
<figure><img src="image_5.png" /><figcaption>Fig. 3. Average (b)–(d) and standard deviation (c)–(e) of the displacement computed for dot diameters varying from 1 pixel to 10 pixels while keeping the mean distance between the centers of the dots equal to 16 pixels. (b) and (c) focus on the comparison of the BOS algorithm, whereas (e) and (f) compare comBOS with PIV algorithms. (a): sample of the speckle.</figcaption></figure>

hence many small dots. To emphasize this result, fronts of equal dot density were manually added to Fig. 4(b) for 3.5, 6.5, and 13 dots per subset with dashed gray lines. The reason why covering the same area of the speckle with large dots is less effective than with smaller dots is understood as begin a consequence of the lower values of the derivatives of the gray level at the outer limits of the large dots than those of the thinner dots.

These same results are now analyzed from the perspective of the dot density in Fig. 4(c) instead of the uniformity ratio. It shows the variation of the noise $(\delta x'^2)^{1/2}$ with the number of dots per subset for the same speckles as in Fig. 4(b) with the additional results from DaVis 10 and PIVlab. For dot densities lower than 10, the noise tends to increase with smaller dots, for both comBOS and DaVis 10. This corresponds to the results observed in the region of large uniformity ratio in Fig. 4(b).

Besides, the noise becomes less sensitive to the dot diameter above the density of 12 dots per subset, this configuration being that on the low uniformity ratio side of in Fig. 4(b) where the noise is found to be weakly influenced by the uniformity ratio. The $\phi_{\text{dot}} = 2.2$ case is however associated with a slightly higher noise level consistent with the conclusion derived from Fig. 3 in which it was shown that dots smaller than $\phi_{\text{dot}} = 3$ pixel should be avoided. In contrast, PIVlab is only slightly affected by the number of dots per subset.

In summary, it can be concluded that an optimal speckle least-squares algorithms should be designed such that the distance between the dots results in at least 12 dots in the subsets used to calculate the displacements, and that the diameter of the dots is sufficient to ensure that less than 25% of the speckle area is uniform.

### 4.3. Effect of image contrast

Images that are used in the context of BOS are acquired with a digital camera, and thus the gray levels are affected by quantization

error. This error influences both the gray level of individual pixels of the images and the corresponding spatial derivatives, while both of these quantities are needed to estimate the displacements in Eq. (14). The relevant parameter to assess the ability of the algorithm to accurately process low-contrast images is the number of gray levels $n_{sg}$ between the darkest and the brightest point of the image. The darkest spots are found at the center of the dots and the brightest are in the background, or vice versa for speckles consisting of white dots and black background. This quantity is mainly determined by the contrast of the speckle itself and the bit depth of the camera sensor. For instance, a speckle which is printed with glossy black ink may reflect a significant portion of the incident light, and thus diminish the contrast. Also, for the same speckle and camera exposure, $n_{sg}$ will be 16 times larger for a camera recording 12-bits images rather than 8-bits ones.

A pair of reference and deformed images was generated to analyze the effect of $n_{sg}$ on the accuracy of the calculated displacement. A sample of the reference image is shown in Fig. 5(a). The value $n_{sg}$ is increased linearly from a value of 13 on the left of the image to 256 on the right side. The value $n_{sg} = 256$ is the maximum achievable in comBOS because it stores images as 8-bits .bmp files. The size of the full speckle is 1500 × 1500 pixel, mean $(\delta x)$ and rms $(\delta x'^2)^{1/2}$ displacement statistics can therefore be estimated for each value of $n_{sg}$ from the 60 vertically aligned 25 × 25 pixel large subsets. Fig. 5(b) shows that the expected displacement of 0.25 pixel is well captured by DaVis 10 for the full range of $n_{sg}$ values tested. For comBOS however the displacement is overestimated by up to $1 \times 10^{-3}$ pixel for $n_{sg} < 50$ in addition to a bias of $5 \times 10^{-4}$ to $1 \times 10^{-3}$ pixel constant over the range of $n_{sg}$ tested. This bias is due to the effect of the linearly varying background color across the subsets, which seems not to be handled correctly by comBOS. Both comBOS and DaVis 10 suffer from a significant increase in noise $(\Delta x'^2)^{1/2}$ for $n_{sg}$ smaller than 50 as shown in Fig. 5(c).
<figure><img src="image_6.png" /><figcaption>**Fig. 4.** Standard deviation of the displacement against the uniformity ratio (b) and dot density (c) for three different diameters of dots. (a): sample of a speckle with a density varying from 3 dots per subset to 66 dots per subset.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Fig. 5.** Average (b) and standard deviation (c) of the displacement computed for contrast varying from 13 to 256 gray level with $\Phi_{\text{dot}} = 3.3$ pixels and a density of 12 dots per subsets. (a): sample of the speckle.</figcaption></figure>

The noise level obtained with PIVlab follows the same trend as that obtained with comBOS or DaVis 10. However, the average error made by PIVlab is significantly more affected by the contrast, with a bias increasing by $5 \times 10^{-3}$ pixel between $n_{sg} = 200$ and $n_{sg} = 13$. Nonetheless, it is fair and important to remember that in the context of PIV, a bias of a few hundredths of a pixel would not be noticeable.

It is concluded from these observations that optimal results will be obtained with comBOS and DaVis 10 if the difference between the darkest and the lightest spots of the speckles in all subsets is at least 50 shades of gray, and ideally 100.

### 4.4. Effect of noise level

The synthetic images tested so far are all free of noise, which is not realistic. In fact, images captured by a camera are altered by various types of noise. The two main families are spatial noises and temporal noises. The former are caused by the specific response of each individual pixel in the camera, they are random in space but fairly constant in time, they can therefore be mitigated through a calibration. The latter include shot noise, dark current noise, and readout noise, which are random in time. These noise sources cannot be compensated

for, however the signal-to-noise ratio of the image can be improved by increasing the amount of light collected by the sensor, either by increasing the exposure time and aperture of the camera, or the illumination of the speckle. The sensitivity of the displacement estimation algorithms of comBOS, DaVis 10, and PIVlab, is tested with a reference speckle and a speckle shifted by $\delta_x = 0.25$ pixel. They are both altered by adding an increasing amount of random Gaussian noise from the left to the right of the images. A noise-to-signal ratio $N/S$ is defined as the local standard deviation of the noise divided by the $n_{sg}$ value of the speckle. Here, $n_{sg}$ is set constant and equal to $n_{sg} = 153$ before the noise is added. The noise level is increased from $N/S = 0\%$ on the left side of the image to $N/S = 7\%$ on the right side of the images. A sample of the image is displayed in Fig. 6(a). The statistics for $\langle \delta x \rangle$ and $\langle \delta x'^2 \rangle^{1/2}$ are obtained from 80 vertically aligned 25 × 25 pixel large subsets. Fig. 6(b) shows that DaVis 10 is more robust to noise than comBOS which suffers from a significant underestimation of $\langle \delta x \rangle$ for noise levels higher than 4% of $n_{sg}$. However, it is clear from Fig. 6(c) that the noise in the images should be kept as low as possible because it linearly propagates to the computed displacements. For example, $\langle \delta x'^2 \rangle^{1/2}$ is multiplied by 7 with DaVis 10 and by 9 with comBOS for a noise-to-signal ratio of 4% in comparison to a clean image. A similar noise propagation is observed
<figure><img src="image_8.png" /><figcaption>**Fig. 6.** Average (b) and standard deviation (c) of the displacement computed for Noise to Signal ratio from 0 to 7% with $\Phi_{\text{dot}} = 3.3$ pixels and a density of 12 dots per subsets. (a): sample of the speckle.</figcaption></figure>

for PIVlab, with a slope more than twice as steep as for comBOS and DaVis 10.

## 5. Performance comparison between comBOS and LaVision DaVis 10

The effect of the speckle parameters on the quality of the displacement estimation has been assessed in the previous section with a unidirectional and uniform displacement of 0.25 pixel. This constitutes an ideal case and is useful for establishing a set of general rules for the design of a speckle. In practical applications, deformations due to the gradients of the refractive index effect can be large and associated with normal and shear strains. The ability of comBOS and DaVis 10 to compute accurate displacements is assessed in the following for these realistic conditions.

### 5.1. Error in the estimation of large displacements

Three pairs of speckles are generated to determine the accuracy of comBOS and DaVis 10 with horizontal displacements up to 2 pixels. In all cases, the imposed horizontal displacement $\Delta x^*$ is increased linearly from $\Delta x^* = 0$ pixel on the left side of the 3000 × 1500 pixel large image to $\Delta x^* = 2$ pixel at the right border. The vertical displacement is kept equal to zero for two pairs of speckles while the dot diameter is 3.3 pixels in one case and 5 pixels in the second. The third pair consists of 3.3 pixel dots, but a vertical displacement equal to the local horizontal displacement has been imposed to ensure that the displacement vector always points 45° upward.

The results in terms of computed displacement error and noise are respectively shown in Figs. 7(a) and (b). These figures depict a periodic pattern in the error and noise with a period of one pixel. This result is expected with comBOS and PIVlab because integer pixel displacements are compensated before the sub-pixel displacements are estimated, and although the information is not available in the DaVis 10 documentation, a similar strategy is probably implemented as the noise presented in Fig. 7(b) follows the same behavior. For comBOS and DaVis 10, the maximum error and noise in $\Delta x$ is observed for displacements close to $\Delta x = 0.5 + n$ pixel where $n$ is an integer. With PIVlab, the error is nearly symmetrical with respect to a zero bias with an amplitude similar as that observed for comBOS with the smaller dots. It is minimum for $\Delta x = n/2$. The noise level obtained with PIVlab is less sensitive to $\Delta x$ than with the other two algorithms, and it is at an intermediate level between DaVis 10 and comBOS. Overall, DaVis 10 performs better than comBOS and PIVlab in terms of noise, and is not noticeably biased by the displacement in average.

Focusing on the comBOS results first shows that the vertical and horizontal displacements are well decoupled because the results obtained for unidirectional and bidirectional displacements are nearly identical. A maximum bias of $-6 \times 10^{-3}$ pixel is observed, it corresponds to -1.2% of $\Delta x^* = 0.5$ pixel. It is also clear that increasing the size of the dots reduces the error, for instance for $\Phi_{\text{dot}}$ increasing from 3.3 pixel to 5 pixel, the error drops by a factor of 3. This behavior is a consequence of the reduction in the magnitude of the curvature of the gray level map with larger dots, and therefore the improved accuracy in estimating the gray level derivatives. One can also observe that the maximum error and maximum noise do not occur at exactly $\Delta x^* = 0.5 + n$ pixel but for slightly higher $\Delta x^*$. This lack of symmetry, which is not detrimental, is a consequence of the software implementation of the bilinear interpolation of the deformed subsets between the two passes of the sub-pixel displacement calculation.

### 5.2. Effect of normal and shear strains of the speckle

In all previous examples, the displacement of the deformed image was uniform, or quasi-uniform at the scale of a subset for the benchmark proposed in the Section 5.1. However, unlike the original algorithm of Pan et al. [31] which was specifically designed to handle normal and shear strains of the speckle in the deformed image, the simplified version implemented in comBOS is not optimized for these conditions. A benchmark is therefore set up to test the performances of comBOS in comparison to DaVis 10 and PIVlab in terms of the resulting accuracy and noise affecting the displacement in the presence of normal and shear strains. A second objective is to determine whether the relevant parameter to quantify the robustness of the displacement calculation method under these conditions is the shear rate $\gamma = d\Delta x^*/dy$ and the normal strain $\sigma = d\Delta x^*/dx$, or the amplitude of the displacement within a subset $X_{\text{amp}} = \Delta x_{\text{max}}^* - \Delta x_{\text{min}}^*$. For that purpose, images were generated with five values $X_{\text{amp}}$ ranging from 0.01 pixel to 0.56 pixel for subsets of sizes 25, 35, and 51 pixels. For each condition two images are considered, one with $\Delta x^*$ varying vertically to test shear effects, and the second with $\Delta x^*$ varying horizontally for normal strains. The average displacement within a subset is always set to 0 pixel and the evolution from the minimum to the maximum displacements is linear. n addition, subsets are defined to contain at least 10 dots, which, is large enough according to the results presented in Section 4.2 to remove any significant variation in the noise level due to the number of dots per subset. Also, unless a specific diameter is mentioned, the diameter of the dots considered in the following is 3.3 pixels.

The Average and root mean square displacements ($\langle \delta_x \rangle$ and $\langle (\delta_x')^2 \rangle^{1/2}$) are calculated from a minimum of 50 displacement vectors for the 51 pixel large subsets, and 100 vectors for the 25 pixel large subsets.
<figure><img src="image_9.png" /><figcaption>Fig. 7. Average (a) and standard deviation (b) of the displacement error computed for imposed displacement varying from 0 to 2 pixels with a density of 12 dots per subsets.</figcaption></figure>

Results are presented in Fig. 8 as a function of $γ$ and $σ$ only. Indeed, the perfect match in average and noise between results obtained for the different subset sizes (see Fig. 8(b) and (d)) indicates that $γ$ and $σ$ are the relevant scaling parameters for comBOS and DaVis 10, and consequently, that $X_{amp}$ is not an appropriate scaling. As a second comment, the results in Fig. 8 show that the effect of normal and shear strain on the displacement calculation is similar. This observation was expected because all pixels of the subsets are considered individually by the algorithm and the derivatives are computed identically in both horizontal and vertical directions. Finally, Fig. 8 highlights a nearly linear relationship between the gradients in displacements $γ$ and $σ$, and the displacement noise. It should also be noted that despite an overall lower noise level with DaVis 10, all the three algorithms behave similarly. As a consequence, the regions of the BOS images affected by strong gradients of displacements, and thus regions that correspond to large curvatures of the refractive index, should be regarded with caution because they are affected by a locally large uncertainty. Unfortunately, regions of strong gradients are generally also regions of small spatial extent, such as mixing layers or shocks, so it is not possible to assess the local noise level as it would be in a larger region where the noise can be estimated by looking at the deviations of the displacements between adjacent subsets compared to the overall trend of the spatial evolution of the displacements. However, this uncertainty can be slightly reduced by locally decreasing the diameter of the dots according to the trend observed in Fig. 8(d), which also shows noise values for dots with diameters of 2.6 and 5 pixels. This suggests that strategies based on local optimization of the speckle properties from a priori knowledge of the flow features could help improve the accuracy of the method.

### 5.3. Realistic test case

This last assessment of the displacement computation algorithm is based on a realistic displacement map that contains multiscale features such as large scale modulations, turbulence, and shock waves. In fact, in this test case the reference image is deformed according to a conventional Schlieren image of a supersonic underexpanded jet (see Mercier et al. [34] for the details about the experimental configuration). The gray level values of the original Schlieren image have been scaled to have a mean of 0 and a standard deviation of 0.2. The reference image is a real camera image that is the same size as the Schlieren image, so it can be deformed directly by shifting each pixel by the scaled gray level value of the corresponding pixel. Unlike BOS, only one component of the density gradients, i.e. the displacements, is available with conventional Schlieren. Here, it is the horizontal component. To mimic a realistic test case where displacements are observed in both directions, the same displacement is applied in both $x$ and $y$ directions. This procedure deviates slightly from a real flow because the displacements in orthogonal directions are generally not so well correlated, but it should be still stress the algorithm in a realistic way.

The BOS displacement map is then computed from the reference speckle and the deformed one with $15 \times 15$ pixel large subsets that contains approximately 9 to 10 dots of diameter 3 pixel. A 50% overlap is also applied between adjacent subsets to increase the final resolution.

The conventional Schlieren image is presented in Fig. 9(a) and the corresponding BOS image computed by comBOS in Fig. 9(b). Overall, the main features of the conventional Schlieren image are well preserved in the BOS image, however fine scale structures are filtered out. No noticeable difference is found when observing the displacement map obtained with comBOS, DaVis 10 and PIVlab, so only the one obtained with comBOS is shown.

A more quantitative analysis of the BOS process is proposed in Fig. 9(c) which compares the profiles of displacements obtained with the three methods with the original Schlieren image along a probed line. The displacements are calculated at the center of each subset, so unlike a conventional Schlieren image, the BOS displacements are only known at discrete positions. A nearest-neighbor interpolation is chosen for the display of the calculated displacements to emphasize the lack of resolution of a BOS image compared to a conventional Schlieren image. Discrepancies in the order of 0.1 pixels are observed for the three methods in regions where the gray level of the Schlieren image depicts a large curvature such as at positions $x = 168$ mm or $x = 195$ mm. These discrepancies are large in comparison to the expected uncertainty observed in the ideal cases proposed in this study. They are the consequence of the spatial averaging imposed by the displacement calculation process, and should be a concern for the quantitative analyses of physical phenomena such as shock waves unless a modern single pixel BOS is employed as in Sugisaki et al. [16].

## 6. Conclusion

Background Oriented Schlieren (BOS) is one of the options in the fluid mechanist toolkit to visualize and analyze compressible or thermally driven flows. However, in the opinion of the authors, one of the great potentials of BOS is for teaching activities, as it offers the possibility to observe many of the concepts taught in fluid dynamics with a limited amount of equipment. This aspect of the BOS is certainly underutilized, perhaps because, to the authors' knowledge, there is no dedicated free and a low-noise solution for the processing of BOS images. This motivated the authors to develop an algorithm for the computation of the BOS displacements following the principles of the Lucas-Kanade method which rests upon linear least-squares optimization that can be handled with limited computational resources. This algorithm is included in the Matlab application comBOS that is found to perform with a similar level of accuracy as the module dedicated to BOS implemented in DaVis 10, and is more accurate and less prone to produce noise than the free PIV tool PIVlab in most applications. The discrepancies between the results of the two BOS algorithms are indeed limited to a maximum of 1.2% for displacements close to 0.5 pixel. PIVlab, however, is affected a deviation of 2% from the reference for
<figure><img src="image_10.png" /><figcaption>**Fig. 8.** Average (a, c) and standard deviation (b, d) of the displacement computed against the shear stress (a, b) and normal stress (c, d) for various size of subsets. $\phi_{\text{dot}} = 3.3$ except for the two specific cases in (d) in dashed line.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Fig. 9.** Imposed displacement map (a) from a Z-type Schlieren image of a supersonic jet [34] and the computed image from synthetic BOS images (b). (c): displacement profiles extracted from the lines in blue (a) and red (b). Subsets are 15 × 15 pixels large and contain 9 to 10 dots of diameter $\phi_{\text{dot}} = 3$.</figcaption></figure>

displacements of 0.25 pixel, together with a significantly higher level of noise than with comBOS and DaVis 10 when the speckles deviate from the optimal conditions. Determining these optimal conditions is also one of the objectives of this article, the analysis of which allowed to establish a set of general rules regarding the speckle and the quality of the images that are valid for both comBOS and DaVis 10. They can be used by students, lecturers, or researchers who want to perform BOS measurements using comBOS or least-squares based displacement calculations. The systematic comparison with PIVlab will also help those who want to perform BOS on images with larger distortions, for which a tool based on cross-correlation is better suited. First of all, the apparent diameter of the dots on the camera does not have any significant effect as long as it is larger than 3 pixels, and that it is not too large such that the subsets used for the computation of BOS displacements can contain at least 10 dots. This range of suitable diameters is found to be significantly wider than with PIV tools such as PIVlab. It is also demonstrated that the noise on the images directly propagates to the noise level of the results, so it is important to ensure that the images captured by the camera are bright enough to keep the noise at a low level. Ultimately, the accuracy of the displacements calculated by BOS or PIV methods has been shown to be questionable in regions where the displacements are subject to large gradients and even more where the curvature is substantial. In a future work it would be worthwhile to add a term in comBOS that handles the effect of the gradients in the formulation of the least-squares problem as done by Pan et al. [31]. This could be set as an option in the algorithm to the avoid the cost of increasing the computational complexity, which is detrimental to educational uses, while providing the best accuracy for research applications in complex flows.

CRediT authorship contribution statement

Bertrand Mercier: Writing, Methodology, Software. Tom Lacassagne: Conceptualization, Writing.
Declaration of competing interest

The authors declare that they have no known competing finan-
cial interests or personal relationships that could have appeared to
influence the work reported in this paper.

Data availability

Data will be made available on request.

Acknowledgments

This work was performed within the framework of the SPINEL
project partially funded by UNIT(Université Numérique Ingénierie et
Technologie) PUNCHY project, and the GFA (German-French Academy
for the industry of the future) call for education and knowledge.

References

[1] S.B. Dalziel, G.O. Hughes, B.R. Sutherland, Whole-field density measurements by 'synthetic schlieren', Exp. Fluids 28 (4) (2000) 322-335, http://dx.doi.org/10.1007/s003480050391.

[2] M. Raffel, H. Richard, G. Meier, On the applicability of background oriented optical tomography for large scale aerodynamic investigations, Exp. Fluids 28 (5) (2000) 477-481, http://dx.doi.org/10.1007/s003480050408.

[3] M.J. Hargather, G.S. Settles, Natural-background-oriented schlieren imaging, Exp. Fluids 48 (1) (2010) 59-68, http://dx.doi.org/10.1007/s00348-009-0709-3.

[4] A. Bauknecht, C.B. Merz, M. Raffel, Airborne application of the background oriented schlieren technique to a helicopter in forward flight, in: 17th International Symposium on the Application of Laser and Imaging Techniques to Fluid Mechanics, 2014, pp. 07-10.

[5] E.A. Haering Jr., M.A. Hill, Background oriented schlieren using celestial objects, 2017, US Patent 9, 599, 497.

[6] L. Venkatakrishnan, G. Meier, Density measurements using the background oriented schlieren technique, Exp. Fluids 37 (2) (2004) 237-247, http://dx.doi.org/10.1007/s00348-004-0807-1.

[7] E. Goldhahn, J. Seume, The background oriented schlieren technique: sensitivity, accuracy, resolution and application to a three-dimensional density field, Exp. Fluids 43 (2) (2007) 241-249, http://dx.doi.org/10.1007/s00348-007-0331-1.

[8] A. Schröder, B. Over, R. Geisler, A. Bulit, R. Schwane, J. Kompenhans, Measurements of density fields in micro nozzle plumes in vacuum by using an enhanced tomographic background oriented schlieren (BOS) technique, in: 9th International Symposium on Measurement Technology and Intelligent Instruments, Saint-Petersburg, 2009.

[9] U. Hartmann, R. Adamczuk, J. Seume, Tomographic background oriented schlieren applications for turbomachinery, in: 53rd AIAA Aerospace Sciences Meeting, (2015-1690) 2015, http://dx.doi.org/10.2514/6.2015-1690.

[10] F. Nicolas, V. Todoroff, A. Plyer, G. Le Besnerais, D. Donjat, F. Micheli, F. Cham-
pagnat, P. Cornic, Y. Le Sant, A direct approach for instantaneous 3D density
field reconstruction from background-oriented schlieren (BOS) measurements,
Exp. Fluids 57 (1) (2016) 1-21, http://dx.doi.org/10.1007/s00348-015-2100-x.

[11] S. Cai, Z. Wang, F. Fuest, Y.J. Jeon, C. Gray, G.E. Karniadakis, Flow over an espresso cup: inferring 3-D velocity and pressure fields from tomographic background oriented schlieren via physics-informed neural networks, J. Fluid Mech. 915 (2021) A102, http://dx.doi.org/10.1017/jfm.2021.135.

[12] B. Atcheson, W. Heidrich, I. Ihrke, An evaluation of optical flow algorithms for background oriented schlieren imaging, Exp. Fluids 46 (3) (2009) 467-476, http://dx.doi.org/10.1007/s00348-008-0572-7.

[13] L.P. Mendes, A.M. Ricardo, A.J. Bernardino, R.M. Ferreira, A comparative study of optical flow methods for fluid mechanics, Exp. Fluids 63 (1) (2022) 7, http://dx.doi.org/10.1007/s00348-021-03357-7.

[14] B.O. Cakir, S. Lavagnoli, B.H. Saracoglu, C. Fureby, Assessment and application of optical flow in background-oriented schlieren for compressible flows, Exp. Fluids 64 (1) (2023) 11, http://dx.doi.org/10.1007/s00348-022-03553-z.

[15] T. Shimazaki, S. Ichihara, Y. Tagawa, Background oriented schlieren technique with fast Fourier demodulation for measuring large density-gradient fields of fluids, Exp. Therm Fluid Sci. 134 (2022) 110598, http://dx.doi.org/10.1016/j.expthermflusci.2022.110598.

[16] H. Sugisaki, C. Lee, Y. Ozawa, K. Nakai, Y. Saito, T. Nonomura, K. Asai, Y. Matsuda, Single-pixel correlation applied to background-oriented schlieren measurement, Exp. Fluids 63 (1) (2022) 1-7, http://dx.doi.org/10.1007/s00348-021-03373-7.

[17] G. Elsinga, B. Van Oudheusden, F. Scarano, D. Watt, Assessment and application of quantitative schlieren methods: Calibrated color schlieren and background oriented schlieren, Exp. Fluids 36 (2) (2004) 309-325, http://dx.doi.org/10.1007/s00348-003-0724-8.

[18] M.J. Hargather, G.S. Settles, A comparison of three quantitative schlieren techniques, Opt. Lasers Eng. 50 (1) (2012) 8-17, http://dx.doi.org/10.1016/j.optlaseng.2011.05.012.

[19] M. Raffel, Background-oriented schlieren (BOS) techniques, Exp. Fluids 56 (3) (2015) 1-17, http://dx.doi.org/10.1007/s00348-015-1927-5.

[20] K. Hayasaka, Y. Tagawa, Mobile visualization of density fields using smartphone background-oriented schlieren, Exp. Fluids 60 (11) (2019) 1-15, http://dx.doi.org/10.1007/s00348-019-2817-z.

[21] F. Cozzi, E. Göttlich, Enhanced background oriented schlieren (EBOS), J. Phys.: Conf. Ser. 1249 (1) (2019) 012017, http://dx.doi.org/10.1088/1742-6596/1249/1/012017.

[22] L. Verso, A. Liberzon, Background oriented schlieren in a density stratified fluid, Rev. Sci. Instrum. 86 (10) (2015) 103705, http://dx.doi.org/10.1063/1.4934576.

[23] H. Ben-Gida, R. Gurka, A. Liberzon, OpenPIV-MATLAB—An open-source software for particle image velocimetry; test case: Birds’ aerodynamics, SoftwareX 12 (2020) 100585, http://dx.doi.org/10.1016/j.softx.2020.100585.

[24] W. Thielicke, R. Sonntag, Particle image velocimetry for MATLAB: Accuracy and enhanced algorithms in PIVlab, J. Open Res. Softw. 9 (1) (2021) http://dx.doi.org/10.5334/jors.334.

[25] T. Liu, D.M. Salazar, OpenOpticalFlow-PIV: An open source program integrating optical flow method with cross-correlation method for particle image velocimetry, J. Open Res. Softw. 9 (1) (2021) http://dx.doi.org/10.5334/jors.326.

[26] B. Mercier, S. Hamidouche, R. Gautier, T. Lacassagne, Educational background oriented schlieren based on a matlab app and a smartphone camera, in: 20th International Symposium on Applications of Laser and Imaging Techniques to Fluid Mechanics, Lisbon, 2022.

[27] O. Yevtikhiyeva, N. Skornyakova, A. Udalov, An investigation of the error of the background schlieren method, Meas. Tech. 52 (12) (2009) 1300, http://dx.doi.org/10.1007/s11018-010-9437-6.

[28] N. Vinnichenko, A. Uvarov, Y.Y. Plaksina, Accuracy of background oriented schlieren for different background patterns and means of refraction index reconstruction, in: 15th International Symposium on Flow Visualization, 2012, p. 34.

[29] J.E. Hartberger, Background-Oriented Schlieren Pattern Optimization, Tech. rep., Air Force Institute of Technology, 2011.

[30] W. Merzkirch, 2. Density sensitive flow visualization, in: Fluid Dynamics, in: Methods in Experimental Physics, vol. 18, Academic Press, 1981, pp. 345-403, http://dx.doi.org/10.1016/S0076-695X(08)60547-5.

[31] B. Pan, A. Asundi, H. Xie, J. Gao, Digital image correlation using iterative least squares and pointwise least squares for displacement field and strain field measurements, Opt. Lasers Eng. 47 (7-8) (2009) 865-874, http://dx.doi.org/10.1016/j.optlaseng.2008.10.014.

[32] B.D. Lucas, T. Kanade, An iterative image registration technique with an application to stereo vision, in: Proceedings of Imaging Understanding Workshop, 1981, pp. 121-130.

[33] D.A. Belsley, E. Kuh, R.E. Welsch, Regression Diagnostics: Identifying Influential Data and Sources of Collinearity, John Wiley & Sons, 2005.

[34] B. Mercier, T. Castelain, C. Bailly, Experimental characterisation of the screech feedback loop in underexpanded round jets, J. Fluid Mech. 824 (2017) 202-229, http://dx.doi.org/10.1017/jfm.2017.336.