Medical Ultrasound Tomography Problem: Experimental Data Processing with High-Performance Computing ? ?? Victoria Filatova[0000−0003−3029−1550] and Leonid Pestov[0000−0001−6831−3669] Immanuel Kant Baltic Federal University, Kaliningrad, Russia ViFilatova@kantiana.ru, LPestov@kantiana.ru Abstract. The article is devoted to imaging in Ultrasound Tomography. We describe processing of experimental data obtained from an ultrasonic tomographic prototype in order to reconstruct an image of an acoustical medium. For visualization, we use the well-known geophysical method Reverse time migration (RTM). This work is the continuation of the re- search, earlier the method was tested on synthetic data, the stability of the method was studied and numerical experiments were carried out for 2-dimensional and 3-dimensional models of acoustic medium. The article presents the results of reconstruction based on the data measured for a phantom object using an experimental sample of the ultrasound tomo- graph developed at the Faculty of Physics, Lomonosov Moscow State University. The numerical solution of the problem is a resource consum- ing computational task. We use the standard parallel computing library (OpenMP) and a cluster system to reduce the processing time. Keywords: Ultrasound tomography · Visualization · Reverse time mi- gration · Real data processing · Breast imaging. 1 Introduction The main purpose of the medical ultrasound tomography (UST) is the breast cancer diagnostics. Despite the existing known diagnostic methods (mammog- raphy, computed tomography, MRI), scientists are developing a safe, cheap and simple methods for detecting breast neoplasms. Currently, both the ultrasonic tomographic devices and the methods of processing the data obtained from them are under active research. It is worth noticing several scientific groups provid- ing modern investigation in ultrasound tomography: from Germany under the leadership of N. Ruiter [11], from the USA under the leadership of N. Duric [5] and from Russia under the direction of O. D. Rumiantseva [3]. Existing two- dimensional and three-dimensional ultrasound tomographs use 256 transducers ? This work has been supported by the Russian Science Foundation (grant No. 16-11- 10027). ?? Copyright c 2020 for this paper by its authors. Use permitted under Creative Com- mons License Attribution 4.0 International (CC BY 4.0). 58 V. Filatova, L. Pestov or more. In the two-dimensional case the transducers are located on a circle of radius from 0.1 m. The dominant impulse frequency is from 1 to 8 MHz. From a mathematical point of view, the UST problem is an inverse problem for the wave equation associated with determining the parameters of acoustic medium (speed of sound, attenuation, density) using boundary measurements. The object is probed by an acoustic impulse emitted from the transducer. The inverse data is the measured acoustic field for different positions of the transduc- ers. For medical diagnostics, it is important to reconstruct both the image of the acoustical medium and the acoustical parameters values. The main methods for solving inverse problem for the wave equation are the ray method, the lineariza- tion method (Born approximation) and the optimization method. Typical UST systems use ray-based tomographic algorithms to reconstruct this information from the first-arrival times [10]. A two-step algorithm based on the linearized inverse kinematic problem and the Born approximation was applied in [2]. An- other algorithm, which is considered by O. D. Rumiantseva and her group, is the Grinevich-Novikov algorithm [2], where scattering data of the plane waves in a monochromatic mode are used. N. Duric’s group are developing their own ultrasound tomograph, they uses the Born approximation [20] too. The group of methods from geophysics adapted for ultrasound tomography is used. The ap- proach based on the Kirchhoff migration method is applied in [19]. An iterative method full-waveform inversion (FWI) has been applied to 2D and 3D ultra- sound datasets [1, 17]. The research presented in [13, 14] use the optimization method. We continue our research conducted before. In the papers [6, 7], we proposed a method based on visualizing the acoustical medium using the RTM procedure and determining the values in the inclusions using the kinematic approach. Pre- viously, numerical research (stability study, 2D and 3D experiments) was carried out for simulated data [6–9]. In the current work, the first results of experimental data processing are presented. We consider experimental dataset measured in a real tomographic experiment with phantom object (Department of Acoustics of Moscow State University, [16]). We thank the group O. D. Rumyantseva for the opportunity to work with the real data. The paper [3] presents the results of the same data processing. 2 Medical Ultrasound Imaging Problem 2.1 The Problem Statement Consider the acoustic dynamical system with zero Cauchy data in R2 × (0, T ):  1  c2 pt = div v + δ(x − xs ) r(t),   vt = ∇p, (1)    p|t=0 = 0, v|t=0 = 0. Here p(x, t; xs ) is the acoustic pressure; v(x, t; xs ) is the particle velocity vector field; δ(x − xs ) is Dirac delta function, which models a point source located at Medical ultrasound tomography problem 59 xs ∈ Γ = {x : |x| = R}; r(t) is an impulse (see below Fig. 2); c = c(x) is the speed of sound. Outside the disc Ω = {x : |x| 6 R} the function c is equal to a constant. The waves generated by the boundary sources are reflected/scattered from some objects and recorded on the circle Γ . Let T ∗ is the “acoustical” radius of disc Ω, i. e. minimal time that is required to fill disc Ω with waves initiated by all boundary sources. The registration time T must be greater than T ∗ . The imaging problem consists in the reconstruction of sound-speed image based on the pressure measurements on the boundary for different positions of the sources (inverse data): p0 (x, t; xs ) = p(x, t; xs ), x, xs ∈ Γ, t ∈ [0, T ]. Absorption reconstruction is important for diagnosis as well, but we do not consider this problem here. Remark 1. We consider forward problem (1) without boundary conditions. In practice, an absorbing tire is installed in the ultrasound tomograph, which ab- sorbs the waves. For modeling, we use perfectly matched layers to constrain the computational domain and attenuate waves reflected from the boundary of the inner square (details see in Section 4). 2.2 Visualization The RTM-image is a result of the following procedures: i. We calculate the pair pf (x, t; xs ), v f (x, t; xs ) for known approximation of the sound-speed c0 (x) (forward problem (1)). ii. For a given p0 (inverse data), we solve the reversal time problem:  1  2 pt = div v,   c0 (2)   vt = ∇p + p0 δ(|x| − R)νΓ ,  p|t=T = 0, v|t=T = 0, where νΓ is the outward unit normal vector to Γ , t ∈ (0, T ). Let the pair pb , v b is the solution of the problem (2). iii. For the imaging process, we calculate the so-called “Energy Imaging condi- tion” IE (x) [18]: T XZ IE (x) = [pf pb + c20 (v f , v b )](x, t; xs ) dt, (3) xs 0 where (·, ·) is the inner product in Rn . 60 V. Filatova, L. Pestov We omit the mathematical base of the formula (3). Notice only, that reverse time migration is related to linearized inversion [15]. To explain this formula (naively) let us assume that we deal with only one diffractor. Using kinematic argument it easy to see that waves in the direct and reversal time meet each other right over this diffractor (in the space-time domain, see Fig. 1). Notice also, that the conventional formula for RTM image does not contain integrand c20 (v f , v b ). Fig. 1. The waves in the direct and reversal time in the space-time domain, see for- mula (3). 3 Experiment Description The data were measured in a real tomographic experiment with a phantom object in the form of a boiled egg with two wires threaded (perpendicular to the tomography plane) through it. The phantom object is placed in water (cwater = 1.494 m/ms). Here is a brief description of the tomograph prototype (see details in [3]). There are 256 quasi-point transducers (each works as both a source and a re- ceiver) with a central operating frequency of about 1.25 MHz and located on a circle with a radius of 0.1536 m. When one transducer emits, all the others receive. The emitting impulse r(t) is presented in Fig. 2. The frequency of dig- itizing the signal is 5 MHz, i.e. it is strictly equal to four times the dominant frequency (∆t = 0.0002 ms). Sampling of each signal starts from the 341st sam- ple from the moment of signal emission (the moment of emission is considered the first sample), in order not to accept an uninformative initial part. An ex- ample of the raw signal received by transducers induced by a single source is illustrated in Fig. 3. Figures 3 and 4 have a qualitative meaning, showing the direct and reflected waves with noise. Medical ultrasound tomography problem 61 Fig. 2. Impulse r(t). Fig. 3. Raw data (pressure, single source). 62 V. Filatova, L. Pestov 4 Numerical Experiment The spatial domain is a square 0.32 m ×0.32 m with ∆x = 0.0001 m grid step. The number of nodes is 10240000. The transducers are uniformly located on a circle of radius R = 0.1536 m centered at the origin. The number of transducers is 256. We processed the raw signal using muting the direct waves, band-pass filtering and interpolating in time from 0.0002 ms to 0.00001 ms (see Fig. 4). Time step (∆t) is equal to 0.00001 ms. Fig. 4. Processed data (pressure, single source). The numerical solutions of both the forward problem (1) and the reversal time problem (2) were based on the explicit conditionally convergent Finite- Difference Time-Domain method (FDTD) with a staggered grid [12]. The 12th- order accuracy approximation with respect to space was used and the second- order to time. Spatial and time steps were determined according to the Courant condition, ∆x ∆t < √ , kcmax 2 where the coefficient k depends on the order of the approximation (for the 12th order, k ≈ 1.34). The function c is a constant outside the disc Ω and, therefore, no waves reflect from R2 \ Ω. Thus, we restrict the computational domain by using the well- known technique of perfectly matched layers (PML; see [4]). The computational domain (a square) contains an absorption layer with parameters specially chosen to attenuate waves reflected from the boundary Γ0 (see Fig. 5). We use the Dirichlet condition p|Γ0 =0. Medical ultrasound tomography problem 63 Fig. 5. Scheme of the computational domain. The background speed of sound is known (c0 = cwater ). We applied the RTM procedure and obtained the image (Fig. 6). The RTM image has a lot of oscillation and this is due to the emitting impulse is long-periodic. The image does not contain the wire puncture site and the egg. Fig. 6. RTM image. 64 V. Filatova, L. Pestov Then we applied empirical formula for the imaging process: τ (x,x Z s )+δ X I(x) = [(pb )2 + c20 |v b |2 ](x, t) dt, (4) xs τ (x,xs )−δ where τ (x, xs ) is the travel time between the points x and xs . This formula is explained by the following reasoning. If the impulse brought by waves in the direct and reversal time is long then the integral (3) is nonzero in big enough neighborhood of the diffractor. The shorter impulse, the smaller neighborhood is. But we deal with long-periodic impulse. So we need to use cut-off function in (3) which is non-zero only near the travel time τ (x, xs ). We calculated (4) for δ = 50∆t (see Fig. 7a), δ = 500∆t (see Fig. 7b), and δ = 700∆t (see Fig. 7c). Wire punctures are clearly visible in the image, with the left puncture displayed clearly and the right puncture weaker. This is explained by the fact that one puncture was filled with water, and in the other there was a wire. As δ increases, an egg is detected in the image (see Fig. 7b, c). 5 Conclusions The authors conducted an experimental data processing. A numerical experi- ment was conducted on the cluster system using OpenMP technology. The RTM image contains a lot of noise and it is impossible to identify the egg and wire. The oscillations on the image are explained by the principle of RTM imaging and the long-periodic impulse emitted by the source. We proposed an empirical formula for the imaging process. As a result, the egg and punctures of the wire were detected on the image. In [3], first large-scale inhomogeneities with a low spatial resolution of 0.5 – 1 cm (egg) are reconstructed, then the fine structure of the scatterers with a spatial resolution of 0.25 mm (wire). In the experiment described in this work, both the egg and the wire are reconstructed with a spa- tial resolution corresponding to a grid step of 0.1 mm. The sizes and locations of the scatterers (eggs and wire) are in accordance with values indicated in the work [3]. Acknowledgements The work was supported by the Russian Science Foundation (grant No. 16-11- 10027). We are deeply grateful to O. D. Rumyantseva and D. I. Zotov (Depart- ment of Acoustics of Moscow State University) for sharing the 2D data. Medical ultrasound tomography problem 65 Fig. 7. Resulting images: δ = 50∆t (a), δ = 500∆t (b), δ = 700∆t (c). 66 V. Filatova, L. Pestov References 1. Agudo, O.C., Guasch, L., Huthwaite, P., Warner, M.: 3D imaging of the breast using full-waveform inversion. In: SPIE Proceedings, pp. 99–110. KIT Scientific Publishing (2017) 2. Burov, V.A., Grishina, I.M., Lapshenkina, O.I., Morozov, S.A., Rumyant- seva, O.D., Suhov, E.G: Recovery of the fine structure of the acoustic scatterer in the background of the distorting influence of its large-scale components. Acous- tical Physics 49(6), 738–750 (2003) 3. Burov, V.A., Zotov, D.I., Rumyantseva, O.D.: Reconstruction of the sound ve- locity and absorption spatial distributions in soft biological tissue phantoms from experimental ultrasound tomography data. Acoustical Physics 61, 231–248 (2015) 4. Colino, F., Tsogka, C.: Application of PML absorbing layer model to the linear elastodynamic problem in anisotropic heterogeneous media. Geophysics 66(1), 294–307 (2001) 5. Duric, N., Littrup, P.: Breast ultrasound tomography. In: Kuzmiak, C. (ed.) Breast Imaging, vol. 6. IntechOpen (2017). https://www.intechopen.com/books/breast- imaging/breast-ultrasound-tomography 6. Filatova, V.M., Nosikova, V.V.: Determining sound speed in weak inclusions in the ultrasound tomography problem. Eurasian Journal of Mathematical and Computer Applications 6(1), 11–20 (2018) 7. Filatova, V.M., Nosikova, V.V., Pestov, L.N.: Application of Reverse Time Migra- tion (RTM) procedure in ultrasound tomography, numerical modeling. Eurasian Journal of Mathematical and Computer Applications 4(4), 5–13 (2016) 8. Filatova, V.M., Pestov, L.N., Nosikova, V.V., Rudnitskii, A.G.: Breast ultrasound tomography problem, simulation with noisy model. In: 2018 Days on Diffraction (DD), pp. 106–111. St. Petersburg, Russia (2018) 9. Filatova, V.M., Pestov, L.N., Sedaikina, V.A., Danilin, A.N.: Breast ultrasound tomography problem: 3D simulation. In: 2019 Days on Diffraction (DD), pp. 42– 45. St. Petersburg, Russia (2019) 10. Gemmeke, H., Hopp, T., Zapf, M., Kaiser, C., Ruiter, N.V.: 3D ultrasound com- puter tomography: Hardware setup, reconstruction methods and first clinical re- sults. Nuclear Instruments and Methods in Physics Research A: Accelerators, Spec- trometers, Detectors and Associated Equipment 873, 59–65 (2017) 11. Gemmeke, H., Ruiter, N.V.: 3D ultrasound computer tomography for medical imaging. Nuclear Instruments and Methods in Physics Research A 580(2), 1057–1065 (2007) 12. Ke, B., Zhao, B., Liu, C., Fang, Y.: Wave-equation datuming based on a single shot gather. In: SEG Technical Program Expanded Abstracts, pp. 2156–2159 (2004) 13. Natterer, F., Wubbeling, F.: A propagation-backpropagation method for ultra- sound tomography. Inverse Problems 11(6), 1225–1232 (1995) 14. Natterer, F.: Reflectors in wave equation imaging. Wave Motion 45(6), 776–784 (2008) 15. Op’t Root, T., Stolk, C., de Hoop, M.: Linearized inverse scattering based on seismic reverse time migration. Journal de Mathématiques Pures et Appliquées 98(2), 211–238 (2012) 16. Parkhomenko, P.P., Karavaj, M.F., Sukhov, E.G., Faleev, B.A., Dmitriev, O.V., Drozdov, S.A., Komarov, O.V., Babin, L.V., Popov, A.S., Burov, V.A., Rat- tehl, M.I., Bobov, K.N., Konjushkin, A.L., Rumjantseva, O.D.: Ultrasound tomo- graph and ring antenna array for ultrasound tomograph. Patent RU No. 2145797 C1 . Date of registration: 23.06.1999. Moscow Medical ultrasound tomography problem 67 17. Pratt, R.G., Huang, L., Duric, N., Littrup, P.: Sound-speed and attenuation imag- ing of breast tissue using waveform tomography of transmission and ultrasound data. In: SPIE Proceedings, pp. 1–12. KIT Scientific Publishing (2007) 18. Rocha, D., Tanushev, N., Sava, P.: Acoustic wavefield imaging using the energy norm. In: 2015 SEG Annual Meeting, pp. 49–68. New Orleans, Louisiana (2015) 19. Schmidt, S., Duric, N., Li, C., Roy, O., Huang, Z.F.: Modification of Kirchhoff migration with variable sound speed and attenuation for acoustic imaging of media and application to tomographic imaging of the breast. Medical Physics 38(2), 998–1007 (2011) 20. Simonetti, F., Huang, L., Duric, N., Rama, O.: Imaging beyond the Born ap- proximation: An experimental investigation with an ultrasonic ring array. Physical Review E 76(3) 036601 (2007)