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

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.


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 w, reconstruction is performed on a scalar field, namely, the divergence ∇ · w of w [34,35]. Existing BOS techniques use 3D     [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 (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 ∆Z D is the half-width of the object. The Gladstone-Dale equation [59,60] holds between the refractive index n and the density ρ: where K (= 3.14 × 10 −4 m 3 /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 Similarly, the equation relating the y-direction displacement v to the density gradient is 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] where p 0 is the atmospheric pressure, ρ 0 (= 998 kg/m 3 ) is the density of the ambient fluid (in this case the density of standard-state water), and B (= 314 MPa) and α (= 7) are constants [61].  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, θ). 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, θ) 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.

Density field calculation using vector tomography
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, ∂w/∂x. Furthermore, let the ∂w/∂x distribution be the reconstructed distribution on the xz plane and s(X) be the projected field of ∂w/∂x along the z axis. Since ∂w/∂x has no Y component, neither does the projected field s(X). Assuming that X = x, the relationship between the reconstructed distribution and the projected values is From the characteristics of the reconstructed distribution, ∂w/∂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, ∂w/∂x can be expressed in terms of the radial vector ∂w/∂r in the 2D (r, θ) polar coordinate system and the angle θ with the x axis ( Fig. 3(b)): ∂w/∂x = (∂w/∂r) cos θ. Here, w(r, θ), which is a polar transformation of a scalar field w(x, z), has no θ component, and so ∂w/∂θ is zero. Therefore, hereinafter, ∂w(r, θ)/∂r will be written as ∂w(r)/∂r. In terms of the Cartesian coordinate system, r can be expressed as r = √ x 2 + z 2 and θ as cos θ = x/ √ x 2 + z 2 , and so Eq. 2.6 can be rewritten as As shown in Fig. 3(b), ∂w/∂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 . . . (2.8) Note thatr and N are natural numbers greater than 1, with 1 ≤r ≤ N . First, Eq. 2.7 is rewritten in discrete instead of continuous form: Note thatx andz are natural numbers, and the coordinates (x,z) are discrete coordinates as shown in Fig. 4. The relation between the discretized coordinates (x,z) and the original coordinates (x, z) is At a certain positionx =x 0 , Eq. 2.9 can be written as (2.11) In Eq. 2.11, x 2 0 +z 2 in the term (∂/∂r)w( x 2 0 +z 2 ) is not necessarily natural numbers. On the other hand, in the term (∂/∂r)w(r) of the reconstructed distribution in Eq. 2.8, the variabler (= 1, 2, . . . , N ) is a natural number. When Eq. 2.11 is written in a discrete matrix form such as Eq. 2.8, the real number (∂/∂r)w x 2 0 +z 2 in Eq. 2.11 must be expressed using (∂/∂r)w(r) (r = 1, 2, . . . , N ). By linear interpolation, the real number (∂/∂r)w x 2 0 +z 2 can be written as where bx 0,z is a natural number satisfying bx 0 ,z ≤ x 2 0 +z 2 < bx 0 ,z + 1. When the linear interpolation in Eq. 2.12 is applied to Eq. 2.11, the latter becomes (2.13) If the coefficients of (∂/∂r)w(bx 0 ,z ) are taken to be constants αx 0,z , Eq. 2.13 becomes (2.14) From Eq. 2.13, at some positionx =x 0 , the term bx 0,z in the reconstructed distribution satisfiesx 0 ≤ bx 0,z ≤ N , and so the terms (∂/∂r)w(1) to (∂/∂r)w(x 0 − 1) are not included in the summation. Therefore, the coefficients αx 0,z of the terms (∂/∂r)w(bx 0,z ) in the reconstructed distribution are 0 for 1 ≤z ≤x 0 − 1. Also, the coefficients αx 0,z are determined only by the position in the discrete coordinates (x,z). In the projection s(x) in Eq. 2.14, the variablex takes values among all the natural numbers 1, 2, . . . , and we can then show that which can be written in matrix form as or, explicitly, . . .
As mentioned earlier, the coefficients αx 0 ,z in Eq. 2.14 for the reconstructed distribution at some location x =x 0 are zero for 1 ≤z ≤x 0 − 1. Similarly, the coefficient α 2,z in the reconstructed distribution atx = 2 is 0 forz = 1, and the coefficient α 3,z of the reconstructed distribution atx = 3 is 0 forz = 1, 2. It can be seen that the coefficient matrix A in Eq. 2.16 for the reconstructed distribution for all natural numbersx (x ≥ 1) is always an upper triangular matrix with nonzero and nonsingular diagonal components. It therefore has an inverse matrix A −1 = A/|A|, where A and |A| are the respectively adjugate matrix and determinant of A [63]. Multiplication of both sides of Eq. 2.16a by the inverse matrix A −1 then gives Thus, the reconstructed distribution in matrix form W can be calculated by simply multiplying the matrix of projected values S by the inverse of the coefficient matrix 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 s(x) in Eq. 2.6, and the density gradient field ∂ρ(x, y)/∂x in Eq. 2.3 corresponds to the reconstructed distribution ∂w/∂x in Eq. 2.6. The displacement u and density gradient ∂ρ(x, y)/∂x are represented in a discretized coordinate system. Applying VT using Eq. 2.6 to the displacement u(x,ȳ 0 ) and density gradient ∂ρ(r)/∂r (= ∂ρ(x,ȳ 0 )/∂x) at the positionȳ =ȳ 0 in Eq. 2.3 gives ∂ρ(r) ∂r Note thatr here denotes the radial component of the discretized polar coordinates. By applying VT to all the y displacements u(x,ȳ) on the y axis, the density gradient distribution ∂ρ(x,ȳ)/∂x in the xy section at z = 0 can be calculated. The density field ρ is calculated by integrating in the x direction over the calculated density gradient field ∂ρ(x, y)/∂x: Here, the density of the surrounding fluid ρ 0 is taken to be that of the fluid outside the measurement target (i.e., the shock wave).

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): 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 where n is the iteration step and (i, j) is the position in (x, z) coordinates. ρ n max is the maximum value of the density at the nth 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.

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 µ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 × 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 µ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 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.

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  u is taken to be zero and the density is taken to be that of water at 25 • C under atmospheric pressure, ρ 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 ∂u/∂x + ∂v/∂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.   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.

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.

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 µmm/pixel) multiplied by a factor η. 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.

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.  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 ≤ r/R ≤ 0.98 and 1.03 ≤ r/R ≤ 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 ≤ r/R ≤ 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 ≤ r/R ≤ 0.98 and 1.03 ≤ r/R ≤ 1.05.

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 × 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.  Table 1, and dataset iv comprises the data from [26].

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. 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 η = 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 η 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 η 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 η decreases from 1 to 1/4. In Fig. 13(iv), the maximum pressure from CT-BOS increases about 1.5 times as η 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 η 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 η is reduced to 1/4, the ratio increases by a factor of approximately 5 compared with the case η = 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.

Comparison using experimental data
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 η is reduced from 1 to 1/15. It can be seen that S changes significantly as η decreases below 0.2. For example, for VT-BOS i, S = 1.28 for η = 1, while S = 1.38 for η = 1/15, an increase of 0.1. For VT-BOS ii, S = 1.15 for both η = 1 and 1/15. For VT-BOS iii and VT-BOS iv, S is increased by 0.14 and 0.15, respectively, for η = 1/15 compared with η = 1. As can be seen from Fig. 14, for CT-BOS, S increases by factors of 4-5 as η decreases from 1 to 1/4, whereas for VT-BOS, even when η 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.

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: where ρ 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, θ), as in Figs. 3(b) and 3(c). Eq. 4.1 represents a density field that does not depend on θ and is distributed in the xz section centered at the origin, as shown in Fig. 3. Fig. 16(a) is a plot of ρ(r), and Fig. 16(b) shows the quantities calculated using this density. The density gradient field ∂ρ/∂r is distributed on the xz section (Fig. 3) in the same way as ρ(r). The displacement u is calculated from Eq. 2.3 by integrating ∂ρ/∂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 × 10 −3 m, ξ = 0.5, γ = 1 × 10 −9 m, ε = 5.0 × 10 −4 π/m , and ζ = 2.7 × 10 4 . The displacement 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 ∇ · w of the displacement (i.e., after one differentiation operation). However, as can be seen from Fig. 16(b), the synthetic ∇ · w cannot be calculated from the density ρ, and so the divergence ∇ · w from CT-BOS and the synthetic ∇ · w are not compared. Next, the Laplacian ∇ 2 ρ of the density gradient from the CT-BOS calculation is compared with the synthetic ∇ 2 ρ calculated from the density ρ(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 ρ(r) (Eq. 4.1). (b) Synthetic data that can be generated from the density ρ(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 ∇ · w for reduced spatial resolution factors η = 1/2, 1/3, and 1/4 is lower than that for η = 1. In particular, at 0.95 ≤ r/R ≤ 1.02, the larger the value of η (the greater the spatial resolution), the larger is the divergence of the displacement 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 ∇ 2 ρ from CT-BOS and that from the synthetic data increases as the spatial resolution decreases. For a spatial resolution factor η = 1, the synthetic ∇ 2 ρ and the CT-BOS ∇ 2 ρ are almost identical, but as η decreases from 1/2 to 1/4, the magnitude of the CT-BOS ∇ 2 ρ decreases compared with that of the synthetic ∇ 2 ρ. In CT-BOS, ∇ 2 ρ 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 ∇ · w, (b) Laplacian of density ∇ 2 ρ.

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 η = 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 η = 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 η = 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 η = 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 η =1 and the case η =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 A once. In VT-BOS, the lengths of the rows and columns of the constant matrix 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.

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. 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. 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. 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.