Supercomputer Simulations in Design of Ultrasound Tomography Devices

The paper considers the use of supercomputers in design of medical ultrasound tomography devices. The mathematical models describing the wave propagation in ultrasound tomography should take into account such physical phenomena as diffraction, multiple scattering, and so on. The inverse problem of wave tomography is posed as a coefficient inverse problem with respect to the wave propagation velocity and the absorption factor. Numerous simulations made it possible to determine the optimal parameters of an ultrasound tomograph in order to obtain a spatial resolution of 1.5 mm suitable for early-stage breast cancer diagnosis. The developed methods were tested both on model problems and on real data obtained at the experimental test bench for tomographic studies. The computations were performed on GPU devices of Lomonosov-2 supercomputer at Lomonosov Moscow State University.


Introduction
Modern medical tomographs are complex and expensive devices that can not be designed without extensive mathematical modeling. In ultrasound tomography, the mathematical models used to describe the wave propagation process should take into account such physical phenomena as diffraction, refraction, multiple scattering, and so on. Tomographic image reconstruction involves solving nonlinear large-dimensional inverse problems. The methods developed in the 1970s-1990s for solving inverse problems [8,9] are the most striking mathematical results of the last century. The inverse problem considered in this paper is a coefficient inverse problem [3]. The solution algorithms rely on the processing power of modern GPU clusters.
The concept of resolving power is widely used in the practice of tomographic studies [1]. In quantitative ultrasound tomography [10], both the spatial resolution and the reconstruction accuracy are important. The resolving power of an ultrasound tomograph depends on a large number of parameters, such as the wavelength, the number of sound sources, the number of detectors and distance between them, the frequency spectrum of sounding pulses. For a given set of these parameters, the resolving power depends on the wavefield measurement precision. In this paper we assess the reconstructed image resolution using mathematical modeling, as well as real data obtained in physical experiments.
The images reconstructed from experimental data showed that a spatial resolution of ≈ 1.5 mm is attainable in a low-frequency setup with a central wavelength of ≈ 3 mm, which can be implemented in practice. The algorithms developed by the authors in [3,5] for layer-by-layer ultrasound tomography schemes were used to solve the inverse problems.

Formulation of the Inverse Problem of Wave Tomography
A scalar wave model based on a second-order hyperbolic equation is a simple mathematical model that takes into account the effects of ultrasound diffraction and absorption. According to this model, acoustic pressure u(r, t) satisfies the equation: (1) Here, c(r) = 1/v 2 (r), where v(r) is the speed of sound; r ∈ R 2 is the point in the imaging plane; a(r) is the absorption factor; f (t) describes the sounding pulse emitted from point q; ∆ is the Laplace operator with respect to r. The initial conditions are zero. ∂ n u(r, t)| ST is the normal derivative to the surface S of the domain Ω (Fig. 1), where (r, t) ∈ S × (0, T ); function p(r, t) is known. It is assumed that v(r) = v 0 = const, a(r) = 0 outside of the studied object. This wave propagation model can be used to describe ultrasound waves in soft tissues.
Emitters G S Ω Figure 1. Tomographic examination scheme In this study, we use a layer-by-layer tomography scheme shown in Fig. 1. Object G is insonified using the ultrasound emitters that successively produce sounding pulses. Acoustic pressure U (s, t) is measured at points s of the boundary S for the time interval (0; T ).
The inverse problem of reconstructing the unknown coefficients c(r) and a(r) in equation (1). This inverse problem is ill-posed, and thus we formulate it as a problem of minimizing the residual functional between the measured and numerically simulated wavefields with respect to its argument {c, a}: Here, u(s, t) is the solution of the direct problem (1). We use the iterative gradient method [4] to minimize the functional. Representations of the gradient Φ c (u(c, a)), Φ a (u(c, a)) were obtained in [3,5]. Finite difference time-domain method [6] was used to compute the wavefields.

Numerical Simulations and Experimental Results
By means of numerical simulation, we studied the dependence of the resolving power on the wavelength, the number of emitters and detectors, and other factors. Optimal parameters of an ultrasound tomograph producing a resolving power of ≈ 1.5 mm were determined. Such resolution is sufficient for early-stage breast cancer diagnosis [10]. Figures 2(a) and 2(b), respectively, show the reconstructed speed of sound (SoS) and attenuation images obtained in a numerical simulation. The parameters of the simulated phantom were chosen to match the silicone phantom used in a physical experiment. The minimum distance between the objects is ≈ 1 mm, and the sound speed difference between objects is ≈ 5%. The speed of sound is reconstructed better than the absorption factor. This is a natural consequence of the fact that at sufficiently high frequencies the coefficient of the second derivative in equation (1) is reconstructed better than the coefficient of the first derivative. Figure 2(c) shows the SoS image obtained in a numerical experiment with simulated measurement errors. Figure 2(d) shows the SoS image reconstructed from experimental data. The computations were performed on NVidia Tesla K40s devices of Lomonosov-2 supercomputer [7]. GPU computing has proven to be the most promising technology for solving inverse problems of wave tomography [2]. It takes ≈ 15 minutes to obtain images for a single plane on a single GPU device. The finite difference grid contained 1000 × 1000 points at a pitch of 0.33 mm. The simulation time step was 0.15 µs for a total of 2400 steps. The frequency range of sounding waves was 100-600 kHz with a central wavelength of 3 mm.

Conclusion
It is impossible to design ultrasound tomographs without extensive mathematical modeling. At the design stage it is necessary to perform reconstructions for various phantoms and for a large number of physical parameters of the device. Numerous simulations and processing of experimental data on a supercomputer made it possible to determine the optimal parameters that provide high spatial resolution of ≈ 1.5 mm and high reconstruction accuracy. The physical experiments showed that a total of 20-30 ultrasound emitters are sufficient to achieve such resolution, while the distance between detectors should be approximately half the wavelength.