Medical Ultrasound Tomography Problem: Experimental Data Processing with High-Performance Computing ? ??

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 research, 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 tomograph developed at the Faculty of Physics, Lomonosov Moscow State University. The numerical solution of the problem is a resource consuming computational task. We use the standard parallel computing library (OpenMP) and a cluster system to reduce the processing time.

Keywords

• Ultrasound tomography Visualization gration Real data processing Breast imaging

-

The main purpose of the medical ultrasound tomography (UST) is the breast cancer diagnostics. Despite the existing known diagnostic methods (mammography, 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 scienti c groups providing 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 twodimensional and three-dimensional ultrasound tomographs use 256 transducers 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 eld for di erent positions of the transducers. 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 linearization method (Born approximation) and the optimization method. Typical UST systems use ray-based tomographic algorithms to reconstruct this information from the rst-arrival times [ 10 ]. A two-step algorithm based on the linearized inverse kinematic problem and the Born approximation was applied in [ 2 ]. Another 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 approach based on the Kirchho migration method is applied in [ 19 ]. An iterative method full-waveform inversion (FWI) has been applied to 2D and 3D ultrasound 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. Previously, numerical research (stability study, 2D and 3D experiments) was carried out for simulated data [6{9]. In the current work, the rst 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 2.1

Medical Ultrasound Imaging Problem

The Problem Statement

Consider the acoustic dynamical system with zero Cauchy data in R2 (0; T ): 8 1 >>< c2 pt = div v + (x

vt = rp; > > :pjt=0 = 0; vjt=0 = 0: xs) r(t); (1) Here p(x; t; xs) is the acoustic pressure; v(x; t; xs) is the particle velocity vector eld; (x xs) is Dirac delta function, which models a point source located at xs 2 = fx : jxj = Rg; r(t) is an impulse (see below Fig. 2); c = c(x) is the speed of sound. Outside the disc = fx : jxj 6 Rg the function c is equal to a constant. The waves generated by the boundary sources are re ected/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 ll 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 di erent positions of the sources (inverse data):

p0(x; t; xs) = p(x; t; xs); x; xs 2 ; t 2 [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 absorbs the waves. For modeling, we use perfectly matched layers to constrain the computational domain and attenuate waves re ected 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); vf (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: 8 1 >>< c02 pt = div v; > > : pjt=T = 0; vjt=T = 0; vt = rp + p0 (jxj

R) ; (2) where is the outward unit normal vector to , t 2 (0; T ). Let the pair pb; vb is the solution of the problem (2). iii. For the imaging process, we calculate the so-called \Energy Imaging condition" IE (x) [ 18 ]:

IE (x) =

T X Z xs 0 where ( ; ) is the inner product in Rn.

[pf pb + c02(vf ; vb)](x; t; xs) dt; (3)

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 di ractor. Using kinematic argument it easy to see that waves in the direct and reversal time meet each other right over this di ractor (in the space-time domain, see Fig. 1). Notice also, that the conventional formula for RTM image does not contain integrand c20(vf ; vb). 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 receiver) 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 digitizing 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 sample from the moment of signal emission (the moment of emission is considered the rst sample), in order not to accept an uninformative initial part. An example 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 re ected waves with noise.

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 ltering and interpolating in time from 0:0002 ms to 0:00001 ms (see Fig. 4). Time step ( t) is equal to 0:00001 ms.

The numerical solutions of both the forward problem (1) and the reversal time problem (2) were based on the explicit conditionally convergent FiniteDi erence Time-Domain method (FDTD) with a staggered grid [ 12 ]. The 12thorder accuracy approximation with respect to space was used and the secondorder to time. Spatial and time steps were determined according to the Courant condition, t < x p ; kcmax 2 where the coe cient 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 re ect from R2 n . Thus, we restrict the computational domain by using the wellknown technique of perfectly matched layers (PML; see [ 4 ]). The computational domain (a square) contains an absorption layer with parameters specially chosen to attenuate waves re ected from the boundary 0 (see Fig. 5). We use the Dirichlet condition pj 0 =0.

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. Then we applied empirical formula for the imaging process:

I(x) = X (x;xs)+

Z xs (x;xs) [(pb)2 + c02jvbj2](x; t) dt; (4) 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 di ractor. The shorter impulse, the smaller neighborhood is. But we deal with long-periodic impulse. So we need to use cut-o 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 lled 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 experiment 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 ], rst large-scale inhomogeneities with a low spatial resolution of 0.5 { 1 cm (egg) are reconstructed, then the ne 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 spatial 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-1110027). We are deeply grateful to O. D. Rumyantseva and D. I. Zotov (Department of Acoustics of Moscow State University) for sharing the 2D data. = 50 t (a), = 500 t (b), = 700 t (c).

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 Scienti c Publishing ( 2017 )

2. Burov , V.A. , Grishina , I.M. , Lapshenkina , O.I. , Morozov , S.A. , Rumyantseva , O.D. , Suhov , E.G: Recovery of the ne structure of the acoustic scatterer in the background of the distorting in uence of its large-scale components . Acoustical Physics 49 ( 6 ), 738 { 750 ( 2003 )

3. Burov , V.A. , Zotov , D.I. , Rumyantseva , O.D. : Reconstruction of the sound velocity 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/breastimaging/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 Migration (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 Di raction (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 Di raction (DD) , pp. 42 { 45 . St . Petersburg, Russia ( 2019 )

10. Gemmeke , H. , Hopp , T. , Zapf , M. , Kaiser , C. , Ruiter , N.V.: 3D ultrasound computer tomography: Hardware setup, reconstruction methods and rst clinical results . Nuclear Instruments and Methods in Physics Research A: Accelerators , Spectrometers, 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 ultrasound tomography . Inverse Problems 11 ( 6 ), 1225 { 1232 ( 1995 )

14. Natterer , F.: Re ectors 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 Mathematiques Pures et Appliquees 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. , Rattehl , M.I. , Bobov , K.N. , Konjushkin , A.L. , Rumjantseva , O.D. : Ultrasound tomograph and ring antenna array for ultrasound tomograph . Patent RU No. 2145797 C1 . Date of registration: 23.06 . 1999 . Moscow

17. Pratt , R.G. , Huang , L. , Duric , N. , Littrup , P. : Sound-speed and attenuation imaging of breast tissue using waveform tomography of transmission and ultrasound data . In: SPIE Proceedings , pp. 1 { 12 . KIT Scienti c Publishing ( 2007 )

18. Rocha , D. , Tanushev , N. , Sava , P. : Acoustic wave eld 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. : Modi cation of Kirchho 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 approximation: An experimental investigation with an ultrasonic ring array . Physical Review E 76 ( 3 ) 036601 ( 2007 )