#  Background-oriented Schlieren technique with vector tomography for
measurement of axisymmetric pressure fields of laser-induced underwater shock
waves

Sayaka Ichihara Department of Mechanical Systems Engineering, Tokyo University
of Agriculture and Technology, Koganei Campus 6-204, 2-24-16 Nakacho, Koganei,
Tokyo, Japan  Takaaki Shimazaki  Department of Mechanical Systems Engineering,
Tokyo University of Agriculture and Technology, Koganei Campus 6-204, 2-24-16
Nakacho, Koganei, Tokyo, Japan  Yoshiyuki Tagawa  Department of Mechanical
Systems Engineering, Tokyo University of Agriculture and Technology, Koganei
Campus 6-204, 2-24-16 Nakacho, Koganei, Tokyo, Japan  Institute of Global
Innovation Research, Tokyo University of Agriculture and Technology, Koganei
Campus 6-204, 2-24-16 Nakacho, Koganei, Tokyo, Japan

##### Abstract

This study aims to overcome the problems that existing background-oriented
schlieren (BOS) techniques based on computed tomography (CT-BOS) face when
measuring pressure fields of laser-induced underwater shock waves. To do this,
it proposes a novel BOS technique based on vector tomography (VT-BOS) of an
axisymmetric target. The remarkable feature of the proposed technique is the
reconstruction of an axisymmetric vector field with nonzero divergence, such
as the field of a laser-induced underwater shock wave. This approach is based
on an approximate relation between the projection of the axisymmetric vector
field and the reconstructed vector field. For comparison, the pressure fields
of underwater shock waves are measured with VT-BOS, CT-BOS, and a needle
hydrophone. It is found that VT-BOS is significantly better than CT-BOS in
terms of better convergence, less dependence on the spatial resolution of the
acquired images, and lower computational cost. The proposed technique can be
applied not only to fluid dynamical fields, but also to other axisymmetric
targets in other areas, such as electromagnetics and thermodynamics.
Keywords: Background-oriented schlieren (BOS) technique, vector tomography,
noncontact measurement, pressure field, shock wave

##  1 Introduction

Laser-induced underwater shock waves[1, 2] are widely used in clinical
medicine for various techniques, such as needle-free injection [3, 4, 5, 6,
7], lithotripsy [8, 9], skull base tumor removal [10], and treatment of
diseases such as femoral head osteonecrosis [11, 12, 13], osteochondritis
dissecans [14], osteochondral lesions ([15]), and bone marrow edema [16].
These techniques have the advantage that incisions are not necessary, and thus
there is a lower risk of complications developing [17, 12]. It is therefore of
great value to investigate underwater shock waves in the context of medical
applications. Previous studies have used point measurement systems such as
hydrophones, but these disturb the flow field or impose restrictions on the
experimental conditions [18]. Therefore, noncontact pressure measurements are
required to enable further developments in this field.

Existing noncontact pressure measurement methods include pressure-sensitive
paint (PSP) [19, 20, 21], particle image velocimetry (PIV)[22, 23, 24], and
background-oriented schlieren (BOS)[25, 26] techniques. The BOS technique has
the advantage of allowing measurement with simple equipment such as a
background image and a camera [27, 28, 29]. It has been used to measure a
variety of density fields, for example, those in flames [30], supersonic flows
[31, 32], and surface dielectric barrier discharges [33]. The prominent
feature of the BOS technique is that it can be used for noncontact
measurements of the pressure field of underwater shock waves [26], and, as
described below, it has great scope for further development in this context.

The procedure of the BOS technique for calculating the density and pressure of
a target is shown in Fig.1. The apparent displacement is calculated by
comparing background images with and without the measurement target present.
The calculated apparent displacement is proportional to the integral of the
density gradient with respect to the optical axis of the camera. Therefore, it
is necessary to calculate the three-dimensional (3D) density gradient field
from the apparent displacement by 3D reconstruction.

In conventional 3D reconstruction, instead of the vector field of apparent
displacement $\bm{w}$, reconstruction is performed on a scalar field, namely,
the divergence $\nabla\cdot\bm{w}$ of $\bm{w}$ [34, 35]. Existing BOS
techniques use 3D reconstruction methods, such as filtered back-projection
(FBP)[34, 36, 37] and the algebraic reconstruction technique (ART)[35, 38,
31], which are types of computed tomography (CT). The BOS technique with
computed tomography is hereinafter referred to as CT-BOS. After CT calculation
in CT-BOS, the Laplacian $\nabla^{2}\rho$ ($\propto\nabla\cdot\bm{w}$) of the
density field $\rho$ is obtained and used in the Poisson equation for $\rho$,
which is then solved iteratively. From the linear relationship between density
and pressure, the pressure field is calculated [26]. When the Poisson equation
is solved from the 3D reconstruction of a scalar field by iteration in CT-BOS,
problems arise with regard, for example, to the accuracy of pressure
calculation and convergence of the solution [34, 26], the dependence on
resolution [26, 25], and the high computational cost [39, 40].

To solve these problems, we apply vector tomography (VT), which can directly
reconstruct a 3D vector field (hereinafter referred to as the reconstructed
distribution) from apparent displacements (hereinafter referred to as the
projected distribution), in the BOS technique. VT is a reconstruction method
defined by Norton[41] and Prince[42, 43]. This method extends the central-
slice theorem [44] used in 3D reconstruction of scalar field and realizes the
reconstruction of a vector field when the divergence of the reconstructed
distribution is zero [41, 45]. VT is used, for example, to reconstruct 3D
magnetization vector fields [46], coronal magnetic fields [47], and
electrostatic potential fields [48]. However, the reconstructed distribution
of vector fields such as laser-induced shock waves measured in this paper does
not have a divergence of zero, making it difficult to apply conventional VT to
such fields.

Here we derive an approximate matrix equation relating the reconstructed and
projected distributions by assuming that the reconstructed distribution is an
axisymmetric vector field with only radial components. Based on this equation,
we propose a novel VT that can perform reconstruction without the need to
satisfy the above conditions [49] . It is possible to reconstruct the density
gradient field of an axisymmetric shock wave [50] with respect to the optical
axis of a laser pulse.

By applying this VT method to the BOS technique, the density gradient field
can be obtained from the displacement, which is a vector field. Integration is
performed on the density gradient field to obtain the density field. The
pressure field, which has a linear relationship with the density field, can
then be calculated. This new BOS technique using VT is hereinafter referred to
VT-BOS. We believe that VT-BOS can avoid the calculation procedures that cause
problems in CT-BOS and can overcome these problems.

The remainder of this paper is organized as follows. The theoretical basis of
the proposed VT-BOS method is described in Sec. 2. Experimental measurements
of underwater shock waves are presented in Sec. 3). Using these experimental
data and data from previous studies [26], we investigate the validity of VT-
BOS by comparing the accuracy of the calculated pressure, the convergence of
the calculations, the dependence on spatial resolution, and the computational
cost of VT-BOS with those of CT-BOS in Sec. 4. Our conclusions are presented
in Sec. 5

Figure 1: Calculation procedures for CT-BOS (the previous technique) and VT-
BOS (the proposed technique).

##  2 Theory of BOS technique

###  2.1 Background-oriented schlieren technique

Figure 2: Principle of BOS technique.

Fig. 2 shows a conceptual diagram of the BOS technique [51, 28, 25]. In the
BOS method, the camera, the measurement target, and the background are placed
on a straight line. The camera takes the background image with and without the
measurement target. The apparent local displacement vector $(u,v)$ on the
image is obtained from the two images using an image analysis technique such
as digital image correlation (DIC) [52, 53], optical flow [54, 55], or fast
checkerboard demodulation (FCD) [56, 57, 58]. An $(x,y,z)$ coordinate system
is adopted, with the background image in the $xy$ plane and the camera’s line
of sight in the $z$ direction (Fig. 2). The relationship between the apparent
displacement $u$ in the $x$ direction and the refractive index gradient is
given by

$u=\frac{Z_{D}}{n_{0}}\int^{Z_{D}+\Delta Z_{D}}_{Z_{D}-\Delta
Z_{D}}\frac{\partial n}{\partial x}\,dz,$ (2.1)

[28], where $n$ is the refractive index of the object being measured, $n_{0}$
is the refractive index of the surrounding fluid, $Z_{D}$ is the distance from
the background to the object being measured, and $\Delta Z_{D}$ is the half-
width of the object. The Gladstone–Dale equation [59, 60] holds between the
refractive index $n$ and the density $\rho$:

$n=\rho K+1,$ (2.2)

where $K$ ($=3.14\times 10^{-4}$ m3/kg) is the Gladstone–Dale constant.
Substituting Eq. 2.2 into Eq. 2.1, the equation relating the $x$-direction
displacement to the density gradient is obtained as

$u=\frac{Z_{D}K}{n_{0}}\int^{Z_{D}+\Delta Z_{D}}_{Z_{D}-\Delta
Z_{D}}\frac{\partial\rho}{\partial x}\,dz.$ (2.3)

Similarly, the equation relating the $y$-direction displacement $v$ to the
density gradient is

$v=\frac{Z_{D}K}{n_{0}}\int^{Z_{D}+\Delta Z_{D}}_{Z_{D}-\Delta
Z_{D}}\frac{\partial\rho}{\partial y}\,dz.$ (2.4)

There are two methods for obtaining the density from $u$ and $v$: vector field
3D reconstruction (vector tomography, VT) and scalar field 3D reconstruction
(computed tomography, CT), which are described in Secs. 2.2 and 2.3,
respectively.

Using the calculated density, the pressure is calculated by Tait’s formula
[61, 62]

$\frac{p+B}{p_{0}+B}=\left(\frac{\rho}{\rho_{0}}\right)^{\alpha},$ (2.5)

where $p_{0}$ is the atmospheric pressure, $\rho_{0}$ ($=998$ kg/m3) is the
density of the ambient fluid (in this case the density of standard-state
water), and $B$ ($=314$ MPa) and $\alpha$ ($=7$) are constants [61].

###  2.2 Density field calculation using vector tomography

Figure 3: (a) and (b) Projected values on the $X$ axis and reconstructed
distribution of a vector field in $(x,z)$ coordinates. The blue vectors and
blue zone in (b) show the reconstructed distribution of the vector field
$\partial w/\partial x$. The gray zone in (a) shows the distribution of the
scalar field $w$. The red vectors in (b) show the reconstructed distribution
of the vector field $\partial w/\partial r$ in $(r,\theta)$ coordinates. (c)
Enlarged cross section from the reconstructed distribution and projected
values in (b), where $x\geq 0$ and $X\geq 0$. Figure 4: Relationship between
reconstructed distribution and projected values in a $5\times 5$ matrix. The
enlargement shows specific formulas for linear interpolation on the $x=3$
axis.

In this section, we describe a new VT method for computing density fields.
Three coordinate systems are used: $(X,Y)$, $(x,y,z)$, and $(r,\theta)$.
First, we explain the relationship between these coordinate systems and the
properties of the target vector field.

As shown in Fig. 3(a), the 2D vector field projected from the 3D vector field
is represented in $(X,Y)$ coordinates, and the 3D reconstructed distribution
vector field is represented in $(x,y,z)$ coordinates. The $X$ and $x$ axes and
the $Y$ and $y$ axes are parallel, and the origin of the $(X,Y)$ coordinates
is on the $z$ axis. The vector field in $(X,Y)$ coordinates is the projection
of the reconstructed distribution in $(x,y,z)$ coordinates onto the $z$ axis.
The reconstructed distribution is assumed to be $y$-axis symmetric. Therefore,
the projected vector field is $Y$-axis symmetric. Furthermore, the
reconstructed distributions of all $xz$ sections are assumed to be vector
fields with $z$-axis symmetry and only a radial component. The radial
component is the $r$-direction component in $(r,\theta)$ coordinates, which
are the polar coordinates corresponding to the $(x,z)$ Cartesian coordinates
(see Fig. 3(b)). Based on the above assumptions, the reconstructed
distribution on the $xz$ cross section is not related to the $Y$ component of
the projected distribution, which is a vector field. Therefore, the VT
proposed in this paper focuses only on the $X$ component of the projected
vector field.

Let the reconstructed distribution (vector field) be represented by a
potential (scalar field) $w(x,z)$ (the gray area in Fig. 3(a)) distributed on
the $xz$ plane. That is, the $x$ component of the reconstructed distribution
is represented by the partial derivative of $w(x,z)$ with respect to $x$, that
is, $\partial w/\partial x$. Furthermore, let the $\partial w/\partial x$
distribution be the reconstructed distribution on the $xz$ plane and
$\bm{s}(X)$ be the projected field of $\partial w/\partial x$ along the $z$
axis. Since $\partial w/\partial x$ has no $Y$ component, neither does the
projected field $\bm{s}(X)$. Assuming that $X=x$, the relationship between the
reconstructed distribution and the projected values is

$\bm{s}(x)=\int\frac{\partial w}{\partial x}\,dz.$ (2.6)

From the characteristics of the reconstructed distribution, $\partial
w/\partial x$ is a vector field symmetric about the origin and about the $z$
axis on the $xz$ plane (see Figs. 3(a) and (b)). Also, $\partial w/\partial x$
can be expressed in terms of the radial vector $\partial w/\partial r$ in the
2D $(r,\theta)$ polar coordinate system and the angle $\theta$ with the $x$
axis (Fig. 3(b)): $\partial w/\partial x=(\partial w/\partial r)\cos\theta$.
Here, $w(r,\theta)$, which is a polar transformation of a scalar field
$w(x,z)$, has no $\theta$ component, and so $\partial w/\partial\theta$ is
zero. Therefore, hereinafter, $\partial w(r,\theta)/\partial r$ will be
written as $\partial w(r)/\partial r$. In terms of the Cartesian coordinate
system, $r$ can be expressed as $r=\sqrt{x^{2}+z^{2}}$ and $\theta$ as
$\cos\theta=x/\sqrt{x^{2}+z^{2}}$, and so Eq. 2.6 can be rewritten as

$\bm{s}(x)=\int\frac{\partial}{\partial
r}w\left(\sqrt{x^{2}+z^{2}}\right)\frac{x}{\sqrt{x^{2}+z^{2}}}\,dz.$ (2.7)

As shown in Fig. 3(b), $\partial w/\partial x$ is distributed symmetrically
about the origin and along the $z$ axis, and so we focus only on $xz$ sections
where $x$ is greater than or equal to 0, as shown in Fig. 3(c). We then
consider the final discretization of Eq. 2.7 and express it in terms of an
approximate matrix relationship

$\begin{bmatrix}\bm{s}(1)\\\ \bm{s}(2)\\\ \vdots\\\ \bm{s}(\bar{r})\\\
\vdots\\\ \bm{s}(N)\\\
\end{bmatrix}\propto\begin{bmatrix}\dfrac{\partial}{\partial r}w(1)\\\\[9.0pt]
\dfrac{\partial}{\partial r}w(2)\\\ \vdots\\\ \dfrac{\partial}{\partial
r}w(\bar{r})\\\ \vdots\\\ \dfrac{\partial}{\partial r}w(N)\\\ \end{bmatrix}.$
(2.8)

Note that $\bar{r}$ and $N$ are natural numbers greater than 1, with
$1\leq\bar{r}\leq N$. First, Eq. 2.7 is rewritten in discrete instead of
continuous form:

$\bm{s}(\bar{x})=\frac{\partial}{\partial
r}w\left(\bar{x}\right)+2\sum_{\bar{z}=1}^{N}\frac{\partial}{\partial
r}w\left(\sqrt{\bar{x}^{2}+\bar{z}^{2}}\right)\frac{\bar{x}}{\sqrt{\bar{x}^{2}+\bar{z}^{2}}}.$
(2.9)

Note that $\bar{x}$ and $\bar{z}$ are natural numbers, and the coordinates
$(\bar{x},\bar{z})$ are discrete coordinates as shown in Fig. 4. The relation
between the discretized coordinates $(\bar{x},\bar{z})$ and the original
coordinates $(x,z)$ is

$\begin{split}x=c\bar{x},\\\ z=c\bar{z},\end{split}$ (2.10)

where $c$ [pixel-1] is a constant.

At a certain position $\bar{x}=\bar{x}_{0}$, Eq. 2.9 can be written as

$\bm{s}(\bar{x}_{0})=\frac{\partial}{\partial
r}w(\bar{x}_{0})+2\sum_{\bar{z}=1}^{N}\frac{\partial}{\partial
r}w\left(\sqrt{\bar{x}_{0}^{2}+\bar{z}^{2}}\right)\frac{\bar{x}_{0}}{\sqrt{\bar{x}_{0}^{2}+\bar{z}^{2}}}.$
(2.11)

In Eq. 2.11, $\sqrt{\bar{x}_{0}^{2}+\bar{z}^{2}}$ in the term
$(\partial/\partial r)w(\sqrt{\bar{x}_{0}^{2}+\bar{z}^{2}})$ is not
necessarily natural numbers. On the other hand, in the term
$(\partial/\partial r)w(\bar{r})$ of the reconstructed distribution in Eq.
2.8, the variable $\bar{r}$ ($=1,2,\dots,N$) is a natural number. When Eq.
2.11 is written in a discrete matrix form such as Eq. 2.8, the real number
$(\partial/\partial r)w\left(\sqrt{\bar{x}_{0}^{2}+\bar{z}^{2}}\right)$ in Eq.
2.11 must be expressed using $(\partial/\partial r)w(\bar{r})$
($\bar{r}=1,2,\dots,N$). By linear interpolation, the real number
$(\partial/\partial r)w\left(\sqrt{\bar{x}_{0}^{2}+\bar{z}^{2}}\right)$ can be
written as

$\begin{split}\frac{\partial}{\partial
r}w\left(\sqrt{\bar{x}_{0}^{2}+\bar{z}^{2}}\right)=\frac{\partial}{\partial
r}w\left(b_{\bar{x}_{0},\bar{z}}\right)\left(b_{\bar{x}_{0},\bar{z}}+1-\sqrt{\bar{x}_{0}^{2}+\bar{z}^{2}}\right)\\\
+\frac{\partial}{\partial
r}w\left(b_{\bar{x}_{0},\bar{z}}+1\right)\left(\sqrt{\bar{x}_{0}^{2}+\bar{z}^{2}}-b_{\bar{x}_{0},\bar{z}}\right),\end{split}$
(2.12)

where $b_{\bar{x}_{0},\bar{z}}$ is a natural number satisfying
$b_{\bar{x}_{0},\bar{z}}\leq\sqrt{\bar{x}_{0}^{2}+\bar{z}^{2}}<b_{\bar{x}_{0},\bar{z}}+1$.
When the linear interpolation in Eq. 2.12 is applied to Eq. 2.11, the latter
becomes

$\bm{s}(\bar{x}_{0})=\frac{\partial}{\partial
r}w(\bar{x}_{0})+2\sum_{\bar{z}=1}^{N}\left[\frac{\partial}{\partial
r}w\left(b_{\bar{x}_{0},\bar{z}}\right)\left(b_{\bar{x}_{0},\bar{z}}+1-\sqrt{\bar{x}_{0}^{2}+\bar{z}^{2}}\right)+\frac{\partial}{\partial
r}w\left(b_{\bar{x}_{0},\bar{z}}+1\right)\left(\sqrt{\bar{x}_{0}^{2}+\bar{z}^{2}}-b_{\bar{x}_{0},\bar{z}}\right)\right]$
(2.13)

If the coefficients of $(\partial/\partial r)w(b_{\bar{x}_{0},\bar{z}})$ are
taken to be constants $\alpha_{\bar{x}_{0},\bar{z}}$, Eq. 2.13 becomes

$\bm{s}(\bar{x}_{0})=\sum_{\bar{z}=1}^{N}\alpha_{\bar{x}_{0},\bar{z}}\frac{\partial}{\partial
r}w(b_{\bar{x}_{0},\bar{z}}).$ (2.14)

From Eq. 2.13, at some position $\bar{x}=\bar{x}_{0}$, the term
$b_{\bar{x}_{0},\bar{z}}$ in the reconstructed distribution satisfies
$\bar{x}_{0}\leq b_{\bar{x}_{0},\bar{z}}\leq N$, and so the terms
$(\partial/\partial r)w(1)$ to $(\partial/\partial r)w(\bar{x}_{0}-1)$ are not
included in the summation. Therefore, the coefficients
$\alpha_{\bar{x}_{0},\bar{z}}$ of the terms $(\partial/\partial
r)w(b_{\bar{x}_{0},\bar{z}})$ in the reconstructed distribution are 0 for
$1\leq\bar{z}\leq\bar{x}_{0}-1$. Also, the coefficients
$\alpha_{\bar{x}_{0},\bar{z}}$ are determined only by the position in the
discrete coordinates $(\bar{x},\bar{z})$.

In the projection $\bm{s}(\bar{x})$ in Eq. 2.14, the variable $\bar{x}$ takes
values among all the natural numbers $1,2,\dots$, and we can then show that

$\bm{s}(\bar{x})=\sum_{\bar{z}=1}^{N}\alpha_{\bar{x},\bar{z}}\frac{\partial}{\partial
r}w(b_{\bar{x},\bar{z}}),$ (2.15)

which can be written in matrix form as

$\displaystyle\bm{S}$ $\displaystyle=\bm{A}\bm{W},$ (2.16a) or, explicitly,
$\displaystyle\begin{bmatrix}\bm{s}(1)\\\ \bm{s}(2)\\\ \bm{s}(3)\\\ \vdots\\\
\bm{s}(n)\end{bmatrix}$
$\displaystyle=\begin{bmatrix}\alpha_{1,1}&\alpha_{1,2}&\alpha_{1,3}&\cdots&\alpha_{1,N}\\\
0&\alpha_{2,2}&\alpha_{2,3}&\cdots&\alpha_{2,N}\\\
0&0&\alpha_{3,3}&\cdots&\alpha_{3,N}\\\ \vdots&\vdots&\vdots&\ddots&\vdots\\\
0&0&0&\cdots&\alpha_{N,N}\end{bmatrix}\begin{bmatrix}\dfrac{\partial}{\partial
r}w(1)\\\\[9.0pt] \dfrac{\partial}{\partial r}w(2)\\\\[9.0pt]
\dfrac{\partial}{\partial r}w(3)\\\\[9.0pt] \vdots\\\
\dfrac{\partial}{\partial r}w(N)\end{bmatrix}.$ (2.16b)

As mentioned earlier, the coefficients $\alpha_{\bar{x}_{0},\bar{z}}$ in Eq.
2.14 for the reconstructed distribution at some location $\bar{x}=\bar{x}_{0}$
are zero for $1\leq\bar{z}\leq\bar{x}_{0}-1$. Similarly, the coefficient
$\alpha_{2,\bar{z}}$ in the reconstructed distribution at $\bar{x}=2$ is 0 for
$\bar{z}=1$, and the coefficient $\alpha_{3,\bar{z}}$ of the reconstructed
distribution at $\bar{x}=3$ is 0 for $\bar{z}=1,2$. It can be seen that the
coefficient matrix $\bm{A}$ in Eq. 2.16 for the reconstructed distribution for
all natural numbers $\bar{x}$ ($\bar{x}\geq 1$) is always an upper triangular
matrix with nonzero and nonsingular diagonal components. It therefore has an
inverse matrix $\bm{A}^{-1}=\bm{\widetilde{A}}/|\bm{A}|$, where
$\bm{\widetilde{A}}$ and $|\bm{A}|$ are the respectively adjugate matrix and
determinant of $\bm{A}$ [63]. Multiplication of both sides of Eq. 2.16a by the
inverse matrix $\bm{A}^{-1}$ then gives

$\bm{A}^{-1}\bm{S}=\bm{A}^{-1}\bm{AW}=\bm{W}.$ (2.17)

Thus, the reconstructed distribution in matrix form $\bm{W}$ can be calculated
by simply multiplying the matrix of projected values $\bm{S}$ by the inverse
of the coefficient matrix $\bm{A}$, and the integral relation in Eq. 2.6
between the projected values of the vector field and the reconstructed
distribution of the vector field has been transformed into the matrix form in
Eq. 2.17 by VT.

The integral equation 2.3 shows the relationship between the projected value
of the vector field and the reconstructed distribution of the vector field in
the BOS technique. When applying VT to Eq. 2.3, the displacement $u$ in Eq.
2.3 corresponds to the projected vector field $\bm{s}(x)$ in Eq. 2.6, and the
density gradient field $\partial\rho(x,y)/\partial x$ in Eq. 2.3 corresponds
to the reconstructed distribution $\partial w/\partial x$ in Eq. 2.6. The
displacement $u$ and density gradient $\partial\rho(x,y)/\partial x$ are
represented in a discretized coordinate system. Applying VT using Eq. 2.6 to
the displacement $\bm{u}(\bar{x},\bar{y_{0}})$ and density gradient
$\partial\rho(\bar{r})/\partial\bar{r}$
($=\partial\rho(\bar{x},\bar{y_{0}})/\partial\bar{x}$) at the position
$\bar{y}=\bar{y_{0}}$ in Eq. 2.3 gives

$\frac{\partial\rho(\bar{r})}{\partial
r}=\frac{n_{0}}{Z_{D}K}\bm{A}^{-1}\bm{u}(\bar{x},\bar{y_{0}}).$ (2.18)

Note that $\bar{r}$ here denotes the radial component of the discretized polar
coordinates. By applying VT to all the $y$ displacements
$\bm{u}(\bar{x},\bar{y})$ on the $y$ axis, the density gradient distribution
$\partial\rho(\bar{x},\bar{y})/\partial x$ in the $xy$ section at $z=0$ can be
calculated. The density field $\rho$ is calculated by integrating in the $x$
direction over the calculated density gradient field
$\partial\rho(x,y)/\partial x$:

$\rho=\int^{\infty}_{0}\frac{\partial\rho(x,y)}{\partial x}\,dx$ (2.19)

Here, the density of the surrounding fluid $\rho_{0}$ is taken to be that of
the fluid outside the measurement target (i.e., the shock wave).

###  2.3 Density field calculation using computed tomography

In CT, reconstruction is performed from the divergence of the projected field
[34, 35, 37, 38, 26]. In this paper, we use the analytical method of filter
back-projection (FBP) [64, 65] for CT, following Hayasaka [26]. The underwater
shock wave is $y$-axis symmetric [9, 50] because it was generated by pulsed
laser irradiation in the $y$-axis direction in the coordinate system shown in
Fig. 2. Therefore, the projection onto the plane perpendicular to the $z$ axis
is assumed to be the same from all directions, and reconstruction is performed
using only images taken from one direction, as shown in Fig. 2. Taking the
divergences of the displacements $u$ (Eq. 2.3) and $v$ (Eq. 2.4) and applying
FBP, we obtain the Poisson equation for the density on the plane containing
the axis of symmetry (the $y$ axis) [26] (Fig. 1):

$\frac{\partial^{2}\rho}{\partial x^{2}}+\frac{\partial^{2}\rho}{\partial
y^{2}}=-\frac{2(1+K\rho_{0})}{\Delta Z_{D}K(\Delta
Z_{D}+2Z_{D})}\left(\frac{\partial U}{\partial x}+\frac{\partial V}{\partial
y}\right),$ (2.20)

Let $U$ and $V$ be the displacements of the reconstructed distribution in the
$x$ and $y$ directions, respectively. The FBP calculation uses the iradon
function implemented in MATLAB. The density field is calculated by solving the
Poisson equation 2.20 iteratively (through successive over-relaxation, SOR)
(Fig. 1). The convergence condition for this iterative procedure is defined as

$E(n)=\frac{\rho^{n}_{i,j}-\rho^{n-1}_{i,j}}{\rho^{n}_{\mathrm{max}}}\leq
E_{0},$ (2.21)

where $n$ is the iteration step and $(i,j)$ is the position in $(x,z)$
coordinates. $\rho_{\mathrm{max}}^{n}$ is the maximum value of the density at
the $n$th iteration step. When $E(n)$ is less than a certain value (the
convergence value $E_{0}$) at all positions, the iterative computation is said
to have converged. The determination of the convergence value is discussed in
Sec. 4.1.

##  3 Experiments

###  3.1 Experimental setup

The experimental methods follow those of [26] and [58]. The camera, target,
and background are placed on the same line, as shown in Fig. 5. The grid or
checker pattern background (grid width 20–40 $\mathrm{\SIUnitSymbolMicro m}$)
is placed perpendicular to the camera’s line of sight. The object to be
measured (i.e., an underwater shock wave) is generated between the camera and
the background. The distance from the center of the shock wave to the
background is $Z_{D}$. The underwater shock wave to be measured is a field
with large local displacement, and if the displacement field is too large, its
detection may become difficult. Since the displacement field increases in
proportion to $Z_{D}$ under given measurement conditions (see Eq. 2.3),
$Z_{D}$ should be about 1.5–2.5 mm to ensure that the displacement field does
not become too large. A pulsed laser for illumination (SI-LUX640, Specialised
Imaging, wavelength 640 nm, pulse width 20 ns) is placed behind the background
image. An Nd:YAG pulsed laser (Nano S PIV, Litron Lasers Ltd., wavelength 532
nm, pulse width 6 ns) is focused into the water by an objective lens
(SLMPLN20X, Olympus), generating shock waves from the focal point. A high-
spatial-resolution camera (EOS 80D, Canon, spatial resolution $4000\times
6000$ pixels) is used to take a background image (Fig. 6(a)) in the absence of
a shock wave and a background image (Fig. 6(b)) when a shock wave is generated
(shutter speed 0.5 s). The pulsed lasers for generating the shock wave and for
illuminating the background image are both connected to a delay generator
(Model 575, BNC), and the timing of shooting is adjusted so that the range of
the shock wave radius $R$ is 2.0–3.5 mm. [26] have pointed out that high-
resolution images are necessary for sufficient accuracy of pressure
calculation, and therefore the resolution of the captured images should be
0.9–1.4 $\mathrm{\SIUnitSymbolMicro m}$/pixel. A hydrophone (Muller-Platte
Needle Probe, Mueller Instruments, effective radius $<$0.25 mm) is placed on
the shock wave center axis perpendicular to the camera’s line of sight, about
0.2–0.3 mm from the shock wave radius, to measure the shock wave pressure. In
this experiment, the shock wave pressure peak $P_{\mathrm{max}}$ is measured
under three different conditions. We also use part of the data from [26] as
reference data (see Table 1). To simplify the analysis, the image is rotated
90∘ along the $y$ axis (Figs. 6(a) and 6(b)), as shown in Fig. 5. Two
90∘-rotated images (Figs. 6(a) and 6(b)) are analyzed.

Figure 5: Experimental setup. Table 1: Experimental conditions. Dataset | i | ii | iii | iv
---|---|---|---|---
|  |  |  | [26]
$p_{\mathrm{max}}$ [MPa] | 0.4 | 1.0 | 2.4 | 1.2
$Z_{D}$ [mm] | 1.93 | 2.01 | 1.50 | 8.00
$R$ [mm] | 2.1 | 2.3 | 3.5 | 5.0
Spatial resolution [$\mathrm{\SIUnitSymbolMicro m}$/pixel] | 1.16 | 0.91 | 1.31 | 5.08
Image resolution [pixels] | $2535\times 1499$ | $2947\times 1885$ | $2941\times 1796$ | $1979\times 1046$

###  3.2 Pressure calculation method

In this subsection, we describe the procedure for analyzing the BOS. Fast
checkerboard demodulation (FCD) ([57, 58]) is applied to the two background
images taken in the experiment (Figs. 6(a) and 6(b)) to obtain the
$x$-direction displacement vector field of the shock wave, as shown in Fig.
6(c). In VT-BOS, VT is applied to the displacement vector field in the $x$
direction (Fig. 7(a), left) to calculate the density gradient field. Outside
the displacement and shock wave (Fig. 7(a), center, in red), as boundary
conditions, the displacement field $\bm{u}$ is taken to be zero and the
density is taken to be that of water at 25 ∘C under atmospheric pressure,
$\rho_{0}$. This region outside the shock wave is defined as the region
outside a circle of radius $R+50$ pixels, where $R$ is the radius of the shock
wave. The density and pressure fields (Fig. 7(a), right) are calculated by
integrating the density gradient field. In CT-BOS, CT (Sec.2.3) is applied to
the divergence of the displacement vector field (a scalar field) (Fig. 7(b),
left). Outside the displacement and shock wave (Fig. 7(b), center, in red), as
the boundary condition, the divergence of the displacement $\partial
u/\partial x+\partial v/\partial y$ is taken to be zero. In CT-BOS, after
reconstruction, the density field is calculated by solving the Poisson
equation 2.20 from Sec.2.3 using the SOR method. In the iterative calculation,
the portion of the laser-induced bubble at the center of the shock wave is
taken to have a pressure of zero. The error that defines the convergence
condition 2.21 is also calculated at the same time. After the density field
has been calculated, the pressure field is calculated using Tait’s equation
2.5. The pressures of the shock wave calculated by both techniques are
compared with the pressure measured by the hydrophone.

Figure 6: (a) and (b) Captured background images respectively with and without
the target (shock wave). (c) Vector field of $x$ component measured by
displacement detection (FCD). Figure 7: (a) Left: Vector field of $x$
component (blue vectors). Center: Boundary conditions for VT-BOS. Right:
Pressure field measured by VT-BOS. (b) Left: Divergence of displacement.
Center: Boundary conditions for CT-BOS. Right: Pressure field measured by CT-
BOS.

###  3.3 Evaluation methods

Figure 8: (a) Effective measurement range of the hydrophone compared with the
pressure field measured by BOS techniques. (b) Comparison of pressures
measured by BOS and by a hydrophone on the $x$ axis. The horizontal axis shows
the distance from the center of the shock wave, nondimensionalized by the
radius $R$ of the shock wave. (c) The blue area $S_{B}$ represents the
difference between the pressure measured by the BOS technique and that
measured by the hydrophone, and the gray area $S_{h}$ represents the the
pressure measured by the hydrophone.

####  3.3.1 Method of comparison of pressure calculation accuracy and
convergence of solutions

The accuracies of CT-BOS and VT-BOS pressure calculations are compared based
on the pressure measured by the hydrophone. Unlike the BOS technique, the
hydrophone measures shock wave pressure at a small measurement range.
Therefore, in the pressure field of the shock wave calculated by the BOS
technique, the pressure field in the effective measurement range of the
hydrophone is averaged in the $y$-axis direction (the shaded area in Fig.
8(a)). Then, as shown in Fig. 8(b), the pressure calculated by the BOS
technique is compared with the hydrophone result. The closer the pressure
value from the BOS technique is to the hydrophone measurement, the greater the
measurement accuracy of the BOS technique.

For a quantitative comparison of the accuracy of pressure calculations, we
proceed as follows. We calculate the area $S_{B}$ between the pressure curve
calculated by the BOS method and the measured pressure curve (Fig. 8(c), blue
area) and the area $S_{h}$ between the measured pressure curve and the zero-
pressure line (Fig. 7(d), gray area). We then define the ratio
$S=(S_{B}+S_{h})/S_{h}$. The closer this ratio is to 1, the closer is the
hydrophone pressure measurement to the pressure calculated by the BOS method
and thus the more accurate is the latter calculation.

Investigation of convergence for CT-BOS requires an iterative approach. A
solution is said to be convergent if the convergence condition 2.21 is
satisfied at some iteration. VT-BOS does not require iteration, and the
density gradient can be uniquely calculated from the displacements using a
matrix equation.

####  3.3.2 Reduction of spatial resolution

To investigate the relationship between spatial resolution and the accuracy of
pressure calculation by CT-BOS and VT-BOS, the spatial resolution of the
original image is reduced and the results of pressure calculations by CT-BOS
and VT-BOS are compared. CT-BOS changes the spatial resolution of the
divergence of the displacement field after reconstruction, whereas VT-BOS
changes the spatial resolution of the displacement field itself. The spatial
resolution of the image with reduced resolution is given by that of the
original image (0.91–5.08 $\mathrm{\SIUnitSymbolMicro m}$m/pixel) multiplied
by a factor $\eta$. Specifically, the spatial resolution is reduced using the
interp2 function of MATLAB. Table 1 shows the image resolution and spatial
resolution for the different experimental conditions.

####  3.3.3  Calculation of computational cost

The computational costs are compared using MATLAB’s CPUtime function. This
function returns the total CPU time, allowing comparison of computational cost
independently of PC performance.

##  4 Results and discussion

###  4.1 Pressure calculation accuracy and convergence

####  4.1.1 Pressure calculation accuracy

Figure 9: Pressures of underwater shock wave obtained by BOS and hydrophone:
(a) datasets i–iii; (b) dataset iv.

The pressure waves from VT-BOS, CT-BOS, and the hydrophone are shown in Fig.
9. Note that the calculation criteria for the CT-BOS pressure calculation
results will be discussed later in Sec. 4.1.2.

Fig. 9(a) shows that the VT-BOS pressure wave is closer to the hydrophone
pressure wave than the CT-BOS pressure wave. Especially in the ranges
$0.95\leq r/R\leq 0.98$ and $1.03\leq r/R\leq 1.05$, the pressure wave from
CT-BOS has relatively high and low values compared with that from the
hydrophone. This trend is more pronounced in dataset iii. On the other hand,
in Fig. 9(b), the pressure wave from VT-BOS is similar to that from the
hydrophone in the range $0.95\leq r/R\leq 1.05$. This is because, unlike CT-
BOS, VT-BOS does not involve the application of a differential operator,
making its measurement accuracy less sensitive to noise. Further verification
is described in Sec. 4.2.2.

To quantitatively compare the accuracy of pressure calculation between VT-BOS
and CT-BOS, the results for the ratio $S$ introduced in Sec. 3.3.1 are shown
in Table 2. The ratio for VT-BOS is lower than that for CT-BOS. Furthermore,
the error between the maximum measured pressure and the maximum CT-BOS
pressure is less than 4.15%, as shown in Table 3. On the other hand, the error
for VT-BOS is less than 2.78%. These results show that VT-BOS has higher
accuracy in pressure calculation than CT-BOS. In particular, VT-BOS tends to
calculate pressure more accurately than CT-BOS in ranges of $r/R$ where the
pressure gradient drops rapidly, such as $0.95\leq r/R\leq 0.98$ and $1.03\leq
r/R\leq 1.05$.

Table 2: Ratio $S$ comparing pressure measurements from CT-BOS and VT-BOS with hydrophone measurements. Dataset | i | ii | iii | iv
---|---|---|---|---
|  |  |  | [26]
VT-BOS | 1.30 | 1.14 | 1.18 | 1.27
CT-BOS | 1.47 | 1.17 | 1.34 | 1.70
Table 3: Error [%] in the pressure measurements using CT-BOS and VT-BOS. Dataset | i | ii | iii | iv
---|---|---|---|---
|  |  |  | [26]
VT-BOS | 1.47 | 0.66 | 2.78 | 1.51
CT-BOS | 3.11 | 1.13 | 3.66 | 4.15

####  4.1.2 Convergence

The convergence of the CT-BOS solution is discussed with reference to Figs. 10
and 11 and Table 4. Fig. 10 shows the CT-BOS iteration results for a
convergence value $E_{0}=1.0\times 10^{-8}$[26]. However, for datasets i–iii,
when the convergence condition 2.21 is satisfied, the pressure waveforms
exhibit larger values compared with the hydrophone pressure waveforms. On the
other hand, for dataset iv, the convergence condition 2.21 is satisfied at
iteration step $n=19\,915$, and the pressure waveforms from the hydrophone and
CT-BOS are similar in shape. It is not feasible to plot the waveforms for all
iteration steps in Fig. 10. Therefore, for datasets i–iii, we have plotted the
results for the minimum value of $n$ such that the maximum pressure from CT-
BOS is greater than or equal to the maximum pressure from the hydrophone,
together with the results for two other values of $n$ respectively larger and
smaller than this value. For dataset iv, we have arbitrarily chosen a number
of iteration steps $n=19\,915$ such that the convergence condition 2.21 is
satisfied, together with two other values of $n$ less than 19 915. In all data
in Fig. 10, the pressure from CT-BOS increases with increasing number of
iteration steps. If the maximum pressure from CT-BOS is close to the maximum
pressure measured by the hydrophone, then number of iteration steps required
for convergence depends on the experimental conditions.

Table 4 shows the convergence values $E_{0}$ when the maximum pressure from
CT-BOS is closest to the maximum measured pressure (the red curves in Figs.
10(i)–10(iii)). It can be seen that $E_{0}$ is not the same for the four
datasets (i)–(iv).

Fig. 11 shows $E(n)$ for each number of iteration steps until the convergence
condition 2.21 is satisfied. For datasets i, ii, and iv, $E(n)$ changes only
by an amount of the order of $10^{-8}$ when the number of iteration steps
exceeds 4000, while for dataset iii, $E(n)$ changes only by an amount of the
order of $10^{-7}$ when the number of steps exceeds 1800. Therefore, the CT-
BOS iterations converge, but if the peak of the measured pressure and that of
the pressure from CT-BOS are to coincide, it may be necessary to take into
account not only the convergence condition 2.21, but other conditions as well.

Figure 10: Iteration results for a convergence value $E_{0}=1.0\times 10^{-8}$ for datasets i–iv. The experimental conditions of datasets i–iii are shown in Table 1, and dataset iv comprises the data from [26]. Figure 11: $E(n)$ of datasets i–iv at each iteration step until the convergence condition 2.21 is satisfied.  Table 4: Convergence value $E_{0}$ of each dataset. Dataset | i | ii | iii | iv
---|---|---|---|---
|  |  |  | [26]
$E_{0}$ | $3.59\times 10^{-8}$ | $7.25\times 10^{-8}$ | $5.19\times 10^{-7}$ | $1.0\times 10^{-8}$

###  4.2 Relationship to spatial resolution

The relationship between the accuracy of pressure calculation and the spatial
resolution of VT-BOS and CT-BOS is compared using the ratio $S$. The method of
varying the spatial resolution was described in Sec. 3.3.2. In Sec. 4.2.1, we
investigate the relationship between the spatial resolution and the
measurement accuracy of the BOS technique using the four sets of experimental
data given in Table 1. In Sec. 4.2.2, we use synthetic data to investigate the
reason for the strong dependence of CT-BOS measurement accuracy on spatial
resolution, in contrast to that of VT-BOS.

####  4.2.1 Comparison using experimental data

Fig. 12 shows the pressure field calculated by the BOS technique for the data
of [26]. The upper and lower diagrams show the pressure fields of the shock
wave calculated by VT-BOS and CT-BOS, respectively, from the displacement
fields corresponding to spatial resolution factors $\eta=1$, 1/2, 1/3, and 1/4
from right to left. For VT-BOS, the pressure field remains almost unchanged as
the spatial resolution decreases. On the other hand, for CT-BOS, when the
spatial resolution of the pressure field decreases, an increase in pressure is
observed at the center of the shock wave near the $y$ axis, which is unrelated
to the shock wave and should not occur.

In Fig. 13, the pressures from VT-BOS and CT-BOS are compared for each
dataset. In Figs. 13(i)–13(iv), the VT-BOS pressure waves remain almost
unchanged with changes in spatial resolution. The pressure waves measured by
the hydrophone and those from VT-BOS are almost identical. In Fig. 13(i), for
CT-BOS, as the spatial resolution factor $\eta$ is decreased in the range from
1 to 1/4, the maximum pressure increases about four times, and the measured
pressure at each decrease in spatial resolution no longer coincides with the
CT-BOS pressure. In Fig. 13(ii), the maximum pressure from CT-BOS increases
about 1.4 times as $\eta$ is decreased from 1 to 1/4. As the spatial
resolution decreases, the CT-BOS pressure no longer coincides with the
measured pressure. In Fig. 13(iii), the maximum pressure from CT-BOS increases
by a factor of about three as $\eta$ decreases from 1 to 1/4. In Fig. 13(iv),
the maximum pressure from CT-BOS increases about 1.5 times as $\eta$ decreases
from 1 to 1/4. Thus, in all cases, at each decrease in spatial resolution, the
measured pressure no longer coincides with the CT-BOS pressure.

Fig. 14 shows the ratio $S$, calculated as in Sec. 3.3.1, for CT-BOS and VT-
BOS as the spatial resolution factor $\eta$ is reduced in the range from 1 to
1/4. It can be seen that $S$ increases with decreasing spatial resolution for
CT-BOS. For example, for CT-BOS i, as $\eta$ is reduced to 1/4, the ratio
increases by a factor of approximately 5 compared with the case $\eta=1$.
Similarly, $S$ increases by factors of 1.3, 4 and 3 for CT-BOS ii, CT-BOS iii,
and CT-BOS iv, respectively. On the other hand, for VT-BOS, $S$ remains almost
constant for all dataset, regardless of the change in spatial resolution. In
summary, similar to the results in Fig. 13, the accuracy of the VT-BOS
pressure calculation is hardly affected by the resolution.

To further investigate the effect of a reduction in spatial resolution on the
results from VT-BOS, the ratio $S$ is shown in Fig. 15 as the spatial
resolution factor $\eta$ is reduced from 1 to 1/15. It can be seen that $S$
changes significantly as $\eta$ decreases below 0.2. For example, for VT-BOS
i, $S=1.28$ for $\eta=1$, while $S=1.38$ for $\eta=1/15$, an increase of 0.1.
For VT-BOS ii, $S=1.15$ for both $\eta=1$ and $1/15$. For VT-BOS iii and VT-
BOS iv, $S$ is increased by 0.14 and 0.15, respectively, for $\eta=1/15$
compared with $\eta=1$. As can be seen from Fig. 14, for CT-BOS, $S$ increases
by factors of 4–5 as $\eta$ decreases from 1 to 1/4, whereas for VT-BOS, even
when $\eta$ decreases below 0.2, although $S$ increases slightly, this change
can be regarded as almost negligible.

Therefore, compared with CT-BOS, the accuracy of pressure calculation by VT-
BOS is nearly independent of spatial resolution, and VT-BOS is able to give
good results for the pressure field even when the spatial resolution is
decreased. In this context, the difference between the methods is the
calculation procedure. Differentiation operations are performed twice in CT-
BOS (see Fig. 1), while VT-BOS does not involve any differentiation at all.
Therefore, the dependence of CT-BOS measurement accuracy on resolution might
be attributable to the amplification of noise by the differentiation
operation. This will be investigated further in the next subsection.

Figure 12: Pressure fields of underwater shock wave from dataset iv measured
by VT-BOS (top) and CT-BOS (bottom) at different image resolutions.  Figure
13: Pressure measured by VT-BOS (left) and CT-BOS (right) at different spatial
resolutions. (i)–(iv) show results from datasets i–iv, respectively.  Figure
14: Ratio $S$ at different spatial resolutions (represented by the factor
$\eta$) for datasets i–iv.  Figure 15: Ratio $S$ as the spatial resolution
factor $\eta$ varies from 1 to 1/15 for datasets i–iv.

####  4.2.2 Comparison of results using synthetic data

To investigate in greater detail the cause of the difference between CT-BOS,
whose measurement accuracy depends on spatial resolution, and VT-BOS, whose
accuracy does not, we examine whether the presence of differentiation
operations in the calculation procedure does indeed cause noise amplification.
The experimental data contain noise not only in the displacement field
measured by the BOS technique, but also in the pressure measured by the
hydrophone. Therefore, when comparing datasets, it is difficult to determine
how much of the increased noise is due to amplification by differentiation
operations in the BOS technique and how much is due to the noise in the
hydrophone pressure. The use of synthetic data allows the displacement field
to be subjected to controllable noise, eliminating factors other than the
displacement field’s own noise.

##### Synthetic data

The synthetic data are provided by the following expression simulating the
density field of underwater shock waves:

$\rho(r)=\rho_{0}+\xi\exp\\!\left[-\frac{(r-r_{c})^{2}}{\gamma}\right]+\varepsilon\sin(\zeta
r),$ (4.1)

where $\rho_{0}$ is the density of the surrounding fluid, the exponential term
represents the change in density of the measurement target, and the sinusoidal
term represents noise. Here, $r$ is the radial coordinate in a polar
coordinate system $(r,\theta)$, as in Figs. 3(b) and 3(c). Eq. 4.1 represents
a density field that does not depend on $\theta$ and is distributed in the
$xz$ section centered at the origin, as shown in Fig. 3. Fig. 16(a) is a plot
of $\rho(r)$, and Fig. 16(b) shows the quantities calculated using this
density. The density gradient field $\partial\rho/\partial r$ is distributed
on the $xz$ section (Fig. 3) in the same way as $\rho(r)$. The displacement
$u$ is calculated from Eq. 2.3 by integrating $\partial\rho/\partial r$ along
the $z$ direction. When using the synthetic data from Eq. 4.1), in order to
obtain a density field similar to that of an underwater shock wave, the
following parameter values are adopted: $r_{c}=2.5\times 10^{-3}$ m,
$\xi=0.5$, $\gamma=1\times 10^{-9}$ m, $\varepsilon=5.0\times 10^{-4}\pi$/m ,
and $\zeta=2.7\times 10^{4}$.

The displacement $\bm{w}$ in Fig. 1 can be obtained by the above procedure. In
this subsection, we investigate the effects of the calculations of first and
second derivatives in CT-BOS on the spatial resolution. The spatial resolution
of the displacements generated using the synthetic data is reduced as in Sec.
4.2.1. CT-BOS is applied to these synthetic displacements to investigate the
effect of spatial resolution on the divergence $\nabla\cdot\bm{w}$ of the
displacement (i.e., after one differentiation operation). However, as can be
seen from Fig. 16(b), the synthetic $\nabla\cdot\bm{w}$ cannot be calculated
from the density $\rho$, and so the divergence $\nabla\cdot\bm{w}$ from CT-BOS
and the synthetic $\nabla\cdot\bm{w}$ are not compared. Next, the Laplacian
$\nabla^{2}\rho$ of the density gradient from the CT-BOS calculation is
compared with the synthetic $\nabla^{2}\rho$ calculated from the density
$\rho(r)$ to examine the effect of reducing the spatial resolution when two
differentiation operations have been performed. Finally, we compare the
pressure calculated by CT-BOS with the synthetic $p$ and discuss the impact of
reduced spatial resolution.

Figure 16: (a) Synthetic $\rho(r)$ (Eq. 4.1). (b) Synthetic data that can be
generated from the density $\rho(r)$.

##### Comparison with synthetic data

The CT-BOS calculation is applied to the synthetic data displacement to
investigate the effect of spatial resolution on the derivative calculation.
Fig. 17 show the relationships between the spatial resolution and the
divergence of the displacement (i.e., a first derivative) and the Laplacian of
the density (i.e., a second derivative). From Fig. 17(a), it can be seen that
the divergence $\nabla\cdot\bm{w}$ for reduced spatial resolution factors
$\eta=1/2$, 1/3, and 1/4 is lower than that for $\eta=1$. In particular, at
$0.95\leq r/R\leq 1.02$, the larger the value of $\eta$ (the greater the
spatial resolution), the larger is the divergence of the displacement
$\bm{w}$. The accuracy of the first derivative decreases in proportion to the
decrease in spatial resolution.

From Fig. 17(b), it can be seen that the difference between $\nabla^{2}\rho$
from CT-BOS and that from the synthetic data increases as the spatial
resolution decreases. For a spatial resolution factor $\eta=1$, the synthetic
$\nabla^{2}\rho$ and the CT-BOS $\nabla^{2}\rho$ are almost identical, but as
$\eta$ decreases from 1/2 to 1/4, the magnitude of the CT-BOS $\nabla^{2}\rho$
decreases compared with that of the synthetic $\nabla^{2}\rho$. In CT-BOS,
$\nabla^{2}\rho$ is not calculated by differentiation operations, but by
performing a scalar field 3D reconstruction, and the accuracy of this
reconstruction is decreased when the spatial resolution is reduced.

The above results indicate that when CT-BOS is applied to a displacement field
containing noise, the calculation accuracy decreases as the spatial resolution
decreases in the calculation procedure for the first derivative (taking the
divergence of the displacement) and second derivative (obtaining the Laplacian
of the density field by 3D reconstruction of a scalar field). Unlike CT-BOS,
VT-BOS does not involve differential operations and therefore is not subject
to degradation in measurement accuracy when the spatial resolution is reduced.

Figure 17: (a) Divergence of displacement $\nabla\cdot\bm{w}$, (b) Laplacian
of density $\nabla^{2}\rho$.

###  4.3 Computational cost

The computational costs incurred in calculating the pressure field from the
displacement field using CT-BOS and VT-BOS, respectively, are now compared.
The lower the CPU time, the lower is the computational cost.

Fig. 18 shows the CPU times of both techniques at different spatial
resolutions. For dataset i, the difference in CPU times between CT-BOS and VT-
BOS is about 30 000 s at any spatial resolution. At a spatial resolution
factor $\eta=1$, the CPU time of CT-BOS is 1800 times that of VT-BOS. For
dataset ii, the difference in CPU times is about 25 000 s at all spatial
resolutions. At $\eta=1$, the CPU time of CT-BOS is 710 times that of VT-BOS.
For dataset iii, the difference in CPU times is approximately 58 000 s at all
spatial resolutions. At $\eta=1$, the CPU time of CT-BOS is 1050 times that of
VT-BOS. For dataset iv, the difference in CPU times is about 74 000 s at all
spatial resolutions. At $\eta=1$, the CPU time of CT-BOS is 1300 times that of
VT-BOS. Thus, under all experimental conditions, the computational cost of VT-
BOS is significantly lower than that of CT-BOS.

With VT-BOS, the computational cost also decreases as the spatial resolution
decreases under all experimental conditions. For example, for dataset i, there
is a 15 s difference in CPU times between the case $\eta=$1 and the case
$\eta=$1/4, and for datasets ii and iii, the difference in CPU times between
the two cases is about 22–39 s. With CT-BOS, the computational cost is nearly
constant under all experimental conditions, independent of spatial resolution.

Thus, VT-BOS has an overwhelmingly lower computational cost than CT-BOS,
because, rather than 3D scalar field reconstruction, VT-BOS uses a matrix
equation (Eq. 2.18) for VT, which does not require iterative calculations.
Specifically, in VT, Eq. 2.18 can be computed by computing a square matrix of
constant $\bm{A}$ once. In VT-BOS, the lengths of the rows and columns of the
constant matrix $\bm{A}$ decrease as the spatial resolution decreases, and so
the computational cost decreases as the spatial resolution decreases.

We now discuss the reasons why the computational cost of CT-BOS is independent
of the resolution, referring to Fig. 19, which shows CPU times at different
spatial resolutions for dataset i. These CPU times are for 3D reconstruction
of the scalar field from the displacement field and for calculation of the
pressure from the Laplacian of the density via the Poisson equation. It can be
seen that more than 90% of the total computational cost incurred by CT-BOS is
for 3D reconstruction of the scalar field and that this cost remains constant
regardless of the resolution. On the other hand, the comparatively small cost
associated with the Poisson equation does decrease with decreasing resolution.

From Figs. 18 and 19 reveal that the computational cost of CT-BOS is almost
independent of the spatial resolution, whereas that of VT-BOS does vary with
resolution, depending on the experimental conditions.

Figure 18: CPU times of CT-BOS and VT-BOS at different spatial resolutions for
datasets i–iv.  Figure 19: CPU times of CT-BOS at different spatial
resolutions for dataset i. The pink bars represent the CPU time for 3D
reconstruction of the scalar field from the displacement, and the blue bars
represent the CPU time for calculation of the pressure field from the
Laplacian of the density.

##  5 Conclusions

The purpose of this study was to measure the noncontact pressure field of a
laser-induced underwater shock wave using the background-oriented schlieren
(BOS) technique with a new vector tomography (VT) technique. In VT, the
distribution of a vector field (in this case the refractive index gradient of
an underwater shock wave) is reconstructed from the projected value of a
vector field (in this case the displacement vector), which is integrated along
one direction. Conventional VT can only perform reconstruction if the
divergence of the reconstructed distribution is zero. Therefore, it has been
difficult to apply VT to measurement targets (e.g., underwater shock waves)
that do not satisfy this condition. In this paper, we have assumed that the
reconstructed distribution is axisymmetric and has only a radial component,
and we have constructed an approximate matrix equation for VT that relates the
projected values and the reconstructed distribution.

To evaluate VT-BOS, experimental data and data from a previous study [26] have
been used to compare it with computed tomography (CT)-BOS from three aspects:
accuracy and convergence of pressure calculation, dependence on spatial
resolution, and computational cost. The results can be summarized as follows:

  1. 1.

The accuracy of pressure calculation is higher with VT-BOS than with CT-BOS.
This is especially true in the range where the gradient of the measurement
target changes rapidly. In the case of VT-BOS, the use of the Poisson equation
is avoided and the pressure can be calculated uniquely using a matrix
equation.

  2. 2.

In comparison with the experimental data, the accuracy of the CT-BOS pressure
calculation decreases as the spatial resolution decreases, whereas the
accuracy of the VT-BOS pressure calculation is almost independent of spatial
resolution within the conditions tested in this research. Results obtained
using synthetic data show that the dependence of CT-BOS measurement accuracy
on spatial resolution is due to the use of differentiation operations and 3D
scalar field reconstruction. VT-BOS does not involve either differentiation
operations or 3D scalar field reconstruction, and instead uses linear
interpolation (Eq. 2.12), as a consequence of which the accuracy of pressure
calculation is virtually unaffected by the spatial resolution within the
conditions tested in this research.

  3. 3.

Compared with CT-BOS, VT-BOS incurs a far lower computational cost. This is
because VT-BOS does not require iterative calculations. Also, as a consequence
of point 2, VT-BOS has the advantage of further reducing the computational
load by reducing the requirements on spatial resolution, since its accuracy
remains almost constant even when the spatial resolution is reduced.

For targets that can be assumed to be axisymmetric, such as laser-induced
underwater shock waves, the BOS method with the new VT technique is more
suitable for calculating the pressure field than the BOS method with 3D scalar
field reconstruction. Unlike conventional VT, the technique proposed in this
paper does not require that the divergence of the reconstructed distribution
be zero. Therefore, it can be used not only for measurements in fluids, but
also for measurements of axisymmetric objects with only radial components in
electromagnetics, material mechanics, and other fields.

## References

  * [1] C. Bell and J. Landt, “Laser-induced high-pressure shock waves in water,” Applied Physics Letters, vol. 10, no. 2, pp. 46–48, 1967.
  * [2] W. Lauterborn and A. Vogel, “Shock wave emission by laser generated bubbles,” in Bubble dynamics and shock waves, pp. 67–103, Springer, 2013.
  * [3] Y. Tagawa, N. Oudalov, A. El Ghalbzouri, C. Sun, and D. Lohse, “Needle-free injection into skin and soft matter with highly focused microjets,” Lab on a Chip, vol. 13, no. 7, pp. 1357–1363, 2013.
  * [4] A. Kiyama, N. Endo, S. Kawamoto, C. Katsuta, K. Oida, A. Tanaka, and Y. Tagawa, “Visualization of penetration of a high-speed focused microjet into gel and animal skin,” Journal of Visualization, vol. 22, no. 3, pp. 449–457, 2019\.
  * [5] J. Krizek, F. De Goumoëns, P. Delrot, and C. Moser, “Needle-free delivery of fluids from compact laser-based jet injector,” Lab on a Chip, vol. 20, no. 20, pp. 3784–3791, 2020.
  * [6] Y. Miyazaki, M. Usawa, S. Kawai, J. Yee, M. Muto, and Y. Tagawa, “Dynamic mechanical interaction between injection liquid and human tissue simulant induced by needle-free injection of a highly focused microjet,” Scientific Reports, vol. 11, no. 1, pp. 1–10, 2021.
  * [7] K. Hayasaka, A. Kiyama, and Y. Tagawa, “Effects of pressure impulse and peak pressure of a shockwave on microjet velocity in a microchannel,” Microfluidics and Nanofluidics, vol. 21, no. 11, pp. 1–8, 2017.
  * [8] H. A. Razvi, J. D. Denstedt, S. S. Chun, and J. L. Sales, “Intracorporeal lithotripsy with the holmium: Yag laser,” The Journal of urology, vol. 156, no. 3, pp. 912–914, 1996.
  * [9] G. N. Sankin, Y. Zhou, and P. Zhong, “Focusing of shock waves induced by optical breakdown in water,” The Journal of the Acoustical Society of America, vol. 123, no. 6, pp. 4071–4081, 2008.
  * [10] Y. Ogawa, A. Nakagawa, K. Takayama, and T. Tominaga, “Pulsed laser-induced liquid jet for skull base tumor removal with vascular preservation through the transsphenoidal approach: a clinical investigation,” Acta neurochirurgica, vol. 153, no. 4, pp. 823–830, 2011.
  * [11] J. Ludwig, S. Lauber, H.-J. Lauber, U. Dreisilker, R. Raedel, and H. Hotzinger, “High-energy shock wave treatment of femoral head necrosis in adults,” Clinical Orthopaedics and Related Research®, vol. 387, pp. 119–126, 2001.
  * [12] C.-J. Wang, J.-H. Cheng, C.-C. Huang, H.-K. Yip, and S. Russo, “Extracorporeal shockwave therapy for avascular necrosis of femoral head,” International Journal of Surgery, vol. 24, pp. 184–187, 2015.
  * [13] J. Hausdorf, A. Lutz, S. Mayer-Wagner, C. Birkenmaier, V. Jansson, and M. Maier, “Shock wave therapy for femoral head necrosis—pressure measurements inside the femoral head,” Journal of biomechanics, vol. 43, no. 11, pp. 2065–2069, 2010.
  * [14] S. Heidersdorf, S. Lauber, H. Lauber, H. Hotzinger, J. Ludwig, U. Dreisilker, and R. Rodel, “Osteochondritis dissecans in musculoskeletal shockwave therapy,” Greenwich Medical Ltd ISBN, vol. 1, pp. 255–261, 2000.
  * [15] C. Zhang, H. Huang, L. Yang, and X. Duan, “Extracorporeal shock wave therapy for pain relief after arthroscopic treatment of osteochondral lesions of talus,” The Journal of Foot and Ankle Surgery, vol. 59, no. 1, pp. 190–194, 2020.
  * [16] J. Cao, C. Zhang, H. Huang, L. Yang, and X. Duan, “Bone marrow edema syndrome of the foot treated with extracorporeal shock wave therapy: A retrospective case series,” The Journal of Foot and Ankle Surgery, vol. 60, no. 3, pp. 523–528, 2021.
  * [17] S. Shrivastava et al., “Shock wave treatment in medicine,” Journal of biosciences, vol. 30, no. 2, pp. 269–275, 2005.
  * [18] E. Herbert, S. Balibar, and F. Caupin, “Cavitation pressure in water,” Physical Review E, vol. 74, no. 4, p. 041603, 2006.
  * [19] T. Liu, B. T. Campbell, S. P. Burns, and J. P. Sullivan, “Temperature-and pressure-sensitive luminescent paints in aerodynamics,” Applied Mechanics Reviews, vol. 50, no. 4, pp. 227–246, 1997.
  * [20] J. H. Bell, E. T. Schairer, L. A. Hand, and R. D. Mehta, “Surface pressure measurements using luminescent coatings,” Annual Review of Fluid Mechanics, vol. 33, no. 1, pp. 155–206, 2001.
  * [21] K. Mitsuo, K. Asai, M. Hayasaka, and M. Kameda, “Temperature correction of psp measurement using dual-luminophor coating,” Journal of Visualization, vol. 6, no. 3, pp. 213–223, 2003.
  * [22] B. Van Oudheusden, “Piv-based pressure measurement,” Measurement Science and Technology, vol. 24, no. 3, p. 032001, 2013.
  * [23] N. Fujisawa, Y. Nakamura, F. Matsuura, and Y. Sato, “Pressure field evaluation in microchannel junction flows through $\mu$piv measurement,” Microfluidics and Nanofluidics, vol. 2, no. 5, pp. 447–453, 2006.
  * [24] G. E. Elsinga, F. Scarano, B. Wieneke, and B. W. van Oudheusden, “Tomographic particle image velocimetry,” Experiments in fluids, vol. 41, no. 6, pp. 933–947, 2006.
  * [25] S. Yamamoto, Y. Tagawa, and M. Kameda, “Application of background-oriented schlieren (bos) technique to a laser-induced underwater shock wave,” Experiments in Fluids, vol. 56, no. 5, pp. 1–7, 2015.
  * [26] K. Hayasaka, Y. Tagawa, T. Liu, and M. Kameda, “Optical-flow-based background-oriented schlieren technique for measuring a laser-induced underwater shock wave,” Experiments in Fluids, vol. 57, no. 12, pp. 1–11, 2016.
  * [27] G. Meier, “New optical tools for fluid mechanics,” Sadhana, vol. 23, no. 5, pp. 557–567, 1998.
  * [28] L. Venkatakrishnan and G. Meier, “Density measurements using the background oriented schlieren technique,” Experiments in Fluids, vol. 37, no. 2, pp. 237–247, 2004.
  * [29] K. Hayasaka and Y. Tagawa, “Mobile visualization of density fields using smartphone background-oriented schlieren,” Experiments in Fluids, vol. 60, no. 11, pp. 1–15, 2019.
  * [30] S. J. Grauer, A. Unterberger, A. Rittler, K. J. Daun, A. M. Kempf, and K. Mohri, “Instantaneous 3d flame imaging by background-oriented schlieren tomography,” Combustion and Flame, vol. 196, pp. 284–299, 2018.
  * [31] M. Ota, K. Hamada, H. Kato, and K. Maeno, “Computed-tomographic density measurement of supersonic flow field by colored-grid background oriented schlieren (cgbos) technique,” Measurement Science and Technology, vol. 22, no. 10, p. 104011, 2011.
  * [32] M. Ota, K. Kurihara, K. Aki, Y. Miwa, T. Inage, and K. Maeno, “Quantitative density measurement of the lateral jet/cross-flow interaction field by colored-grid background oriented schlieren (cgbos) technique,” Journal of Visualization, vol. 18, no. 3, pp. 543–552, 2015.
  * [33] Y. Kaneko, H. Nishida, and Y. Tagawa, “Background-oriented schlieren measurement of near-surface density field in surface dielectric-barrier-discharge,” Measurement Science and Technology, vol. 32, no. 12, p. 125402, 2021.
  * [34] L. Venkatakrishnan, “Density measurements in an axisymmetric underexpanded jet by background-oriented schlieren technique,” AIAA journal, vol. 43, no. 7, pp. 1574–1579, 2005.
  * [35] B. Atcheson, I. Ihrke, W. Heidrich, A. Tevs, D. Bradley, M. Magnor, and H.-P. Seidel, “Time-resolved 3d capture of non-stationary gas flows,” ACM transactions on graphics (TOG), vol. 27, no. 5, pp. 1–9, 2008.
  * [36] T. Tipnis, M. V. Finnis, K. Knowles, and D. Bray, “Density measurements for rectangular free jets using background-oriented schlieren,” The Aeronautical Journal, vol. 117, no. 1194, pp. 771–785, 2013.
  * [37] F. Sourgen, F. Leopold, and D. Klatt, “Reconstruction of the density field using the colored background oriented schlieren technique (cbos),” Optics and Lasers in Engineering, vol. 50, no. 1, pp. 29–38, 2012.
  * [38] F. Leopold, M. Ota, D. Klatt, and K. Maeno, “Reconstruction of the unsteady supersonic flow around a spike using the colored background oriented schlieren technique,” Flow Control, Measurement $\&$ Visualization, vol. 1, no. 2, 2013.
  * [39] F. Nicolas, V. Todoroff, A. Plyer, G. Le Besnerais, D. Donjat, F. Micheli, F. Champagnat, P. Cornic, and Y. Le Sant, “A direct approach for instantaneous 3d density field reconstruction from background-oriented schlieren (bos) measurements,” Experiments in fluids, vol. 57, no. 1, p. 13, 2016.
  * [40] I. Ihrke and M. Magnor, “Image-based tomographic reconstruction of flames,” in Proceedings of the 2004 ACM SIGGRAPH/Eurographics symposium on Computer animation, pp. 365–373, 2004.
  * [41] S. J. Norton, “Tomographic reconstruction of 2-d vector fields: application to flow imaging,” Geophysical Journal International, vol. 97, no. 1, pp. 161–168, 1989.
  * [42] J. L. Prince, “Convolution backprojection formulas for 3-d vector tomography with application to mri,” IEEE Transactions on Image Processing, vol. 5, no. 10, pp. 1462–1472, 1996.
  * [43] J. L. Prince, “Tomographic reconstruction of 3-d vector fields using inner product probes,” IEEE transactions on image processing, vol. 3, no. 2, pp. 216–219, 1994.
  * [44] F. Natterer and F. Wübbeling, Mathematical methods in image reconstruction.  SIAM, 2001.
  * [45] I. Svetov, E. Y. Derevtsov, Y. S. Volkov, and T. Schuster, “A numerical solver based on b-splines for 2d vector field tomography in a refracting medium,” Mathematics and Computers in Simulation, vol. 97, pp. 207–223, 2014.
  * [46] C. Donnelly, S. Gliga, V. Scagnoli, M. Holler, J. Raabe, L. J. Heyderman, and M. Guizar-Sicairos, “Tomographic reconstruction of a three-dimensional magnetization vector field,” New Journal of Physics, vol. 20, no. 8, p. 083009, 2018.
  * [47] M. Kramar, B. Inhester, H. Lin, and J. Davila, “Vector tomography for the coronal magnetic field. ii. hanle effect measurements,” The Astrophysical Journal, vol. 775, no. 1, p. 25, 2013.
  * [48] D. Wolf, N. Biziere, S. Sturm, D. Reyes, T. Wade, T. Niermann, J. Krehl, B. Warot-Fonrose, B. Büchner, E. Snoeck, et al., “Holographic vector field electron tomography of three-dimensional nanomagnets,” Communications Physics, vol. 2, no. 1, pp. 1–9, 2019.
  * [49] S. Ichihara, T. Shimazaki, and Y. Tagawa, “Background oriented schlieren technique with reconstruction of three-dimensional vector fields for pressure measurement of laser induced underwater shock waves,” in Proceedings of the 20th international symposium on application of laser and imaging techniques to fluid mechanics, Lisbon, 2022.
  * [50] Y. Tagawa, S. Yamamoto, K. Hayasaka, and M. Kameda, “On pressure impulse of a laser-induced underwater shock wave,” Journal of Fluid Mechanics, vol. 808, pp. 5–18, 2016.
  * [51] G. Meier, “Computerized background-oriented schlieren,” Experiments in fluids, vol. 33, no. 1, pp. 181–187, 2002.
  * [52] P. Meunier and T. Leweke, “Analysis and treatment of errors due to high velocity gradients in particle image velocimetry,” Experiments in fluids, vol. 35, no. 5, pp. 408–421, 2003.
  * [53] F. Moisy, M. Rabaud, and K. Salsac, “A synthetic schlieren method for the measurement of the topography of a liquid interface,” Experiments in Fluids, vol. 46, no. 6, pp. 1021–1036, 2009.
  * [54] B. Atcheson, W. Heidrich, and I. Ihrke, “An evaluation of optical flow algorithms for background oriented schlieren imaging,” Experiments in fluids, vol. 46, no. 3, pp. 467–476, 2009.
  * [55] T. Brox, A. Bruhn, N. Papenberg, and J. Weickert, “High accuracy optical flow estimation based on a theory for warping,” in European conference on computer vision, pp. 25–36, Springer, 2004.
  * [56] K. Hatanaka and T. Saito, “Background oriented schlieren method using multi-scale periodic pattern,” in International Symposium on Shock Waves, pp. 453–458, Springer, 2013.
  * [57] S. Wildeman, “Real-time quantitative schlieren imaging by fast fourier demodulation of a checkered backdrop,” Experiments in Fluids, vol. 59, no. 6, pp. 1–13, 2018.
  * [58] T. Shimazaki, S. Ichihara, and Y. Tagawa, “Background oriented schlieren technique with fast fourier demodulation for measuring large density-gradient fields of fluids,” Experimental Thermal and Fluid Science, p. 110598, 2022\.
  * [59] W. Merzkirch, Flow visualization.  Elsevier, 2012.
  * [60] M. Raffel, “Background-oriented schlieren (bos) techniques,” Experiments in Fluids, vol. 56, no. 3, pp. 1–17, 2015.
  * [61] E. Brujan, Cavitation in Non-Newtonian fluids: with biomedical and bioengineering applications.  Springer Science & Business Media, 2010.
  * [62] J. Richardson, A. Arons, and R. Halverson, “Hydrodynamic properties of sea water at the front of a shock wave,” The Journal of Chemical Physics, vol. 15, no. 11, pp. 785–794, 1947.
  * [63] D. A. Harville, “Matrix algebra from a statistician’s perspective,” 1998.
  * [64] A. Konkani, “Springer handbook of medical technology (2012 edition),” Biomedical Instrumentation & Technology, vol. 47, no. 6, pp. 510–510, 2013\.
  * [65] A. Ziegler, T. Köhler, and R. Proksa, “Noise and resolution in images reconstructed with fbp and osc algorithms for ct,” Medical physics, vol. 34, no. 2, pp. 585–598, 2007.