Ultrasonics 67 (2016) 136–150
Contents lists available at ScienceDirect
Ultrasonics journal homepage: www.elsevier.com/locate/ultras
A computer simulation study of soft tissue characterization using low-frequency ultrasonic tomography A.V. Goncharsky 1, S.Y. Romanov, S.Y. Seryozhnikov ⇑ Lomonosov Moscow State University, Vorobyevy Gory 1, Build. 4, Moscow 119991, Russia
a r t i c l e
i n f o
Article history: Received 26 February 2015 Received in revised form 30 July 2015 Accepted 12 January 2016 Available online 15 January 2016 Keywords: Ultrasonic tomography 3D coefficient inverse problems Wave equation Reflection and transmission tomography GPU clusters
a b s t r a c t We investigate the potential of using ultrasonic diffraction tomography technique for characterization of biological tissues. Unlike most of other studies where ultrasonic tomography operates at frequencies higher than 1 MHz, low-frequency tomography uses lower frequencies on the order of 0.3–0.5 MHz. Such a choice is due to low attenuation at these frequencies, resulting in higher precision of input data. In this paper we explore transmission and reflection schemes for both 2D (layer-by-layer) and 3D tomography. We treat inverse tomography problems as coefficient inverse problems for the wave equation. The time-domain algorithms employed for solving the inverse problem of low-frequency tomography focus on the use of GPU clusters. The results obtained show that a spatial resolution of about 2–3 mm can be achieved when operating at the wavelength of about 5 mm even using a stationary 3D scheme with a few fixed sources and no rotating elements. The study primarily focuses on determining the performance limits of ultrasonic tomography devices currently designed for breast cancer diagnosis. Ó 2016 Elsevier B.V. All rights reserved.
1. Introduction The differential diagnosis of breast cancer is a problem of prime importance in modern medicine. Common medical devices for ultrasonic examination usually employ a reflection-based scheme [1–4]. In the simplest case, where sounding is performed with the emitter at a fixed position, the doctor sees an image in the form of a time sweep of the ultrasonic signal reflected from internal organs within a narrow angle. A transducer array is usually employed to this end. It seems a natural solution to collect the reflected signals by moving the transducer array around the object studied. One can use the resulting data to try to reconstruct the internal structure of the object studied. However, high-quality tomographic images cannot be reconstructed using reflection data alone. The analysis of various tomography schemes and their optimisation is a problem of vital importance whose solution was addressed in many publications [5–7]. Most of the published studies in this field are dedicated to the development of ultrasonic tomography devices operating at the frequencies of 1–3 MHz or higher. The use of ray-based models appears to be a reasonable approach in this frequency domain. Theoretical frameworks have been developed and prototypes have ⇑ Corresponding author. E-mail addresses:
[email protected] (A.V. Goncharsky), romanov60@gmail. com (S.Y. Romanov),
[email protected] (S.Y. Seryozhnikov). 1 Tel.: +7 4959392759. http://dx.doi.org/10.1016/j.ultras.2016.01.008 0041-624X/Ó 2016 Elsevier B.V. All rights reserved.
been made that implement these approaches in 2D layer-by-layer schemes [8–10] and in fully 3D schemes [11,12]. Inverse problems of ultrasonic tomographic diagnostics are nonlinear in contrast to those arising in X-ray tomography. Refraction puts into question the applicability of layer-by-layer models [13]. Attempts to solve inverse problems in the 3D formulation appear to offer better prospects. Three-dimensional inverse problems of ultrasonic tomography, which are much more computationally expensive, can be tackled with computing clusters equipped with graphics processing units (GPU) [14–16]. Attempts have been made to improve the results obtained using ray models [17]. The disadvantage of ray-based models is that they are incapable of describing such wave phenomena as diffraction and multiple scattering. Pratt et al. [18] attempted to estimate the impact of wave effects on the tomographic image reconstruction. Kretzek and Ruiteret al. [15] also try to solve ultrasonic tomography problems in the ray optics approximation. These authors use the velocity structure reconstructed in the ray approximation for reflectivity image reconstruction. The use of models based on wave equations in ultrasonic tomography is a more promising approach. Inverse problems of ultrasonic tomography in wave-based models are nonlinear and many authors use various linearised approximations to solve them. The most common approach involves the use of the so-called Born approximation of the wave equation [19,20]. From the practical point of view, linearised approximations have a rather limited
A.V. Goncharsky et al. / Ultrasonics 67 (2016) 136–150
potential for solving nonlinear problems and can be used only in the neighborhood of the required solution, which is unknown in the case of inverse problems. The study [21] analyses the limiting capabilities of linearised Born-type models. However, from the theoretical point of view, linearised approximations can be used for qualitative analysis of tomography schemes [6,22]. At present, prototype ultrasonic tomography devices have been developed for the use in breast examinations [23–25] with the data interpreted in terms of nonlinear wave models. Wiskin et al. [23] use narrow-band ultrasonic wave sources with frequencies above 1 MHz. Approximate algorithms have been developed for solving the inverse problems of 3D ultrasonic tomography in the wave approximation of the Helmholtz equation. These algorithms have been tested on a prototype tomographic device for breast cancer diagnosis. Both transmission and reflection data are used to solve the inverse problem. The algorithms are implemented as a twostage procedure. In the first stage, the spatial distribution of the speed of sound is reconstructed approximately using a transmission tomography scheme within the framework of a parabolic model in the form of the Helmholtz equation approximation for small refraction angles. In the second stage, an attempt is made to use the reflected signal to improve the derived approximations. Unlike Wiskin et al. [23], in this paper we demonstrate the feasibility of a low-frequency ultrasonic tomography device operating in the 300–500 kHz frequency range. We develop algorithms for solving three-dimensional nonlinear inverse problem in terms of the model of hyperbolic wave equation with no simplifying assumptions. The authors of some studies [24] attempt to reconstruct both the velocity structure and the density distribution. It is shown that the velocity structure can be recovered better than the density distribution. In real objects attenuation is always present. We showed [25] that in the case of low attenuation in a model incorporating both diffraction and attenuation effects the velocity structure can be reconstructed better than the attenuating properties. In the case of low attenuation the reconstructed velocity structure depends only slightly on the attenuation model employed. The velocity structure can be reconstructed fairly well even if the input data errors are of about 5%, whereas the distribution of attenuating properties in the medium cannot be recovered in this case. Similar results follow from the studies of Wiskin et al. [26]. The aim of this study is to develop efficient algorithms for solving 2D and 3D problems of acoustic tomography in terms of wave models. Inverse problems are considered as coefficient inverse problems of the reconstruction of the velocity structure in the diagnosed region. From the mathematical point of view there are two approaches to solving inverse problems of ultrasonic tomography. On the one hand, one can try to develop algorithms using finite-difference schemes for differential equations. An alternative approach is represented by the well-known Green function method, which allows reformulating the inverse problem as a nonlinear operator equation [7,27]. The inverse problem of ultrasonic tomography in the Green function representation has a number of advantages including the possibility of a simple formulation of the problem as a set of nonlinear operator equations. The undoubted advantage of such a formulation is that it requires no boundary conditions to be imposed except for the natural Sommerfeld radiation condition at infinity. However, the integral approach has a very important disadvantage due to the extremely computationally expensive nature of the algorithms involved. Lavarello and Oelze [7] showed that the number of operations in the iterative algorithm proposed for solving the nonlinear problem scales as O(N6) in the 3D case and O(N5) in the 2D case, where N is the number of grid points along one dimension. We obtained a similar result in our study [27], where we used an iterative process based on the Newton method.
137
The strong dependence of the number of operations on the number of grid points forced us to solve inverse problems on a 50 50 50 grid. In real-life applications a 400 400 400 or denser grid is needed to address the 3D problem of ultrasonic tomography, resulting in the increase of the computational time by a factor of several tens of thousands. The situation cannot be remedied even by using a supercomputer. Even with a factor of 1000 speedup compared to a PC, a supercomputer would allow the number of grid points to be increased by a factor of ffiffiffiffiffiffiffiffiffiffiffiffi p 6 1000 3 along each dimension. This is a typical problem for the integral approach. It follows from the above that the development of algorithms that can be run within practically feasible time on modern computers is of great importance in three-dimensional inverse problems of ultrasonic tomography. The breakthrough results in the solution of inverse problems of wave tomography are associated with the studies that demonstrate the possibility of exact computation of the gradient of the residual functional by solving a ‘conjugate’ problem [19,28,29]. In this paper we use gradient-based minimisation methods to solve the nonlinear problem as a coefficient inverse problem for differential equation. The algorithms employed are based on finite-difference time-domain (FDTD) numerical methods and require only O(N4) operations to solve the three-dimensional problem.
2. Formulation of the direct and inverse problems of ultrasonic tomography and numerical algorithms employed in the two- and three-dimensional case Fig. 1 illustrates the arrangement of the sources and detectors for the three-dimensional inverse problem. Number 1 denotes the sources and number 2 denotes the detectors of ultrasonic radiation, which are located on the faces of the cube X. We assume that the object G is located inside the cube X. The remaining space L is filled with water with known sound speed v0. Fig. 2 illustrates the arrangement of the sources and detectors for the layer-by-layer tomography scheme, where the three-dimensional problem is replaced by a set of two-dimensional problems. Numbers 1 and 2 in Fig. 2 denote the sources and detectors, respectively; G is the domain under study, and L is the domain with known sound speed v0. In this study, we address the inverse problem using the wave approximation in the time-domain formulation with point sources. Acoustic pressure field u(r, t) in the domain X RN (N = 2, 3) produced by a point source located at point r0 and generating a pulse described by function f(t), obeys the wave equation:
cðrÞutt ðr; tÞ Duðr; tÞ ¼ dðr r 0 Þ f ðtÞ;
ð1Þ
uðr; t ¼ 0Þ ¼ ut ðr; t ¼ 0Þ ¼ 0;
ð2Þ
@ n ujCT ¼ pðr; tÞ;
ð3Þ
where t is the time, 0 < t < T; u is the acoustic pressure; c(r) = v2(r), v(r) is the sound speed in the medium; D is the Laplace operator with respect to r 2 RN (N = 2, 3), C is the boundary of the domain N; @ nu|CT is the derivative along the normal to the surface C in the region C (0, T), and p(r, t) is a known function. The two-dimensional wave equation (for N = 2) describes a three-dimensional problem that is independent of one of the coordinates. For example, if the object does not vary along one of the dimensions and is sounded with cylindrical waves. The twodimensional version of the model can be used as an approximation for 3D computations if the object under study varies only slightly along one of the dimensions. In medical tomography, where the structure and parameters of an irregularity have to be
138
A.V. Goncharsky et al. / Ultrasonics 67 (2016) 136–150
Let us now consider the problem that is inverse to that defined by Eqs. (1)–(3) and consists of determining some unknown sound speed v(r) based on the experimental data provided by measurements of the acoustic pressure U(c, t) at the boundary C of the domain X during the time (0, T). The time T is usually chosen long enough so that almost all of the reflected and refracted waves would reach the detectors and extend beyond the boundary of the computational domain. Theoretically, the greater is the time T — the higher is the accuracy of the reconstruction. In our numerical simulations, T is chosen to make the received signal U(c, t) negligible at t > T. In the case of exact data U(c, t), the solution of the inverse problem v(r) when substituted into the direct problem (1)–(3) should yield a wave field u(v(r)), so that u(v(r))|CT = U(c, t). Because of the ill-posed nature of the inverse problem and inevitable errors in the experimental measurements, we formulate the inverse problem as the minimisation of the residual functional with respect to an unknown coefficient c(r):
1 2
UðuðcÞÞ ¼ kujCT Uk2 ¼
Fig. 1. The layout of the 3D experiment.
1 2
Z
T
0
Z
C
ðuðc; tÞ Uðc; tÞÞ2 dcdt
ð4Þ
Here U(c, t) is the acoustic pressure measured at the boundary C of the domain X during the time (0, T) and u(r, t) is the solution of the direct problem (1)–(3) for a given c(r) = 1/v2(r). According to [19,27–29], the gradient of the residual functional (4) can be computed using the formula
U0C ðuðcÞÞ ¼
Z
T
wt ðr; tÞut ðr; tÞdt;
ð5Þ
0
where u(r, t) is the solution of the problem (1)–(3) and w(r, t) is the solution of the following ‘conjugate’ problem for a given c(r):
cðrÞwtt ðr; tÞ Dwðr; tÞ ¼ 0;
ð6Þ
wðr; t ¼ TÞ ¼ wt ðr; t ¼ TÞ ¼ 0; @ n wjCT ¼ ujCT U:
ð7Þ ð8Þ
The boundary condition @ nw|CT = 0 is imposed at the points of the boundary C where no experimental data U(c, t) are available. Thus, to compute the gradient U0 C(u), one has to solve both the main (1)–(3) and the ‘conjugate’ (6)–(8) problems. Given U0 C computed using the formula (5), gradient-based methods such as [28] can be used to minimise the residual functional (4). To accommodate the case in which multiple ultrasound sources are used, the following scheme is employed. The experimental data U(l)(c, t) for each source (l) are obtained, and the gradients of residuals U(l) are computed according to (5) separately for each source. The residual U for the whole data set is equal to the sum of the residuals U(l):
U¼
X
UðlÞ ¼
l
Fig. 2. The layer-by-layer scheme.
reconstructed with high precision, the 3D model should be preferred, which treats the wave phenomena in three dimensions [23]. The scalar wave model can describe the propagation of waves in non-attenuating media rather well. It describes the simplest case of pressure waves. In the case of soft tissue studies the scalar model is an approximation, because it assumes that the wave propagation speed depends only on the spatial coordinates, but not on the direction. For fibrous (e.g., muscle) tissues, the wave propagation speed along and across the fibre direction may differ. Such phenomena can be described more accurately by a tensor model. Solving inverse problems in terms of tensor models is too challenging a task even for modern supercomputers.
2 X1 ðlÞ ðlÞ u jCT U : 2 l
Therefore, the gradient U0 C can be computed as the sum of the P gradients obtained for each source: grad U ¼ l grad UðlÞ . For the numeric implementation of the iterative algorithms we used the finite-difference time-domain (FDTD) method. In this formulation, solving differential wave equations reduces to solving difference equations. Let us now present the numerical approximation of the problem (1) for the three-dimensional case. We introduce the following uniform discrete mesh in the domain of spatial arguments (x, y, z) and time t:
vijlk
8 9 ðxi ; yj ; zl ; t k Þ : > > > > > > > > > > ¼ ih; 0 6 i < n; x > > i < = ¼ yj ¼ jh; 0 6 j < n; ; > > > > > > z ¼ lh; 0 6 l < n; > > l > > > > : ; t k ¼ ks; 0 6 k < m
ð9Þ
139
A.V. Goncharsky et al. / Ultrasonics 67 (2016) 136–150
where h is the mesh spacing in the spatial coordinates and s, the time step. We performed our simulations using h = 0.3 mm and s = 78 ns. In the 2D case, we used n = 640 and m = 3200, and in the 3D case, we used n = 320 and m = 1600. We approximate the partial derivatives in the formula (1) using the following secondorder approximation finite differences:
utt ðr; tÞ ¼
k k1 ukþ1 ijl 2uijl þ uijl
s2
; uxx ðr; tÞ ¼
ukiþ1jl 2ukijl þ uki1jl 2
h
:
Similar schemes are used for uyy(r, t) and uzz(r, t). We obtain the following finite-difference scheme for the differential equation (1) in the domain where no sources are present:
cijl
kþ1 k1 uijl 2ukijl þ uijl
s
ukijlþ1
2
2ukijl 2
þ
ukiþ1jl 2ukijl þ uki1jl
ukijl1
h
h
2
ukijþ1l 2ukijl þ ukij1l h
2
¼ 0;
a pre-defined sound speed c(r). The values of u|CT at the detectors were then used as input data U(c, t) for the inverse problem of reconstructing c(r). In all the simulations we used short initial pulses, as shown in Fig. 3(a). Fig. 3(b) shows the frequency spectrum of the initial pulse with a centre wavelength of k = 5 mm. The spectral energy density E is plotted along the vertical axis and the frequency in MHz, along the horizontal axis. As is evident from the figure, the intensity of spectral components at frequencies above 0.6 MHz does not exceed 5% of the peak intensity, which is located at the centre frequency of 0.3 MHz. We used k = 5 mm for 3D simulations and k = 10 mm for 2D simulations. We use the following iterative process to solve the inverse problem. We start with the initial approximation c(0) = c0 = const, which corresponds to the sound speed in pure water v0 = 1.5 km s1. At each iteration (m), the following actions are performed. For each source (l):
or, by isolating the term ukþ1 of the (k + 1)-th time step, we derive ijl
(1) The initial pulse is rendered. (2) The direct problem (1)–(3) is solved using the current iterative approximation c(m). The acoustic pressure u(l)(r, t) is computed using the formula (10). (3) The residual U(l) = U(u(l)) is computed using the formula (4). (4) The ‘conjugate’ problem (6)–(8) for w(l)(r, t) is solved. (5) The gradient grad U(l)(r) is computed using the formula (11).
the following formula for computing the acoustic pressure u(r, t) sequentially in time (the ‘forward-time’ computation):
ukþ1 ¼ ijl
1 2 k s Duijl þ 2ukijl uk1 ijl ; cijl
where Dukijl ¼
ukiþ1jl 2ukijl þuki1jl h2
þ
ð10Þ
ukijþ1l 2ukijl þukij1l h2
þ
ukijlþ1 2ukijl þukijl1 h2
is the dis-
crete Laplacian; ukijl is the value of u(r, t) at point (i, j, l) at time k, and cijl is the value of c(r) at point (i, j, l). The parameters h and s are related by the Courant–Friedrichs–Lewy (CFL) stability condipffiffiffi tion c0:5 s < h= 3 for the 3D problem. We use similar difference formulas to solve the ‘conjugate’ problem (6)–(8). The gradient of the residual functional is computed using the following difference formula:
ðgrad UÞijl ¼
kþ1 m2 X uijl k¼0
k ukijl wkþ1 ijl wijl
s
s
s:
ð11Þ
The difference formulas for the 2D problem can be derived in a similar way, and they differ from the above ones by dropping the subscript l and the variable z. Our numerical simulations use the non-reflecting boundary condition [30] when solving the direct problem (1)–(3). In this case, the boundary conditions (3) acquire the following form:
@ n ujCT ¼ c0:5 @ t ujCT : The numerical simulations consisted of solving the direct problem of computing the acoustic pressure u|CT at the boundary C of the computational domain, where the detectors are located, using
P
The residual UðmÞ ¼
P l
UðlÞ and the gradient grad UðmÞ ðrÞ ¼
ðlÞ
l grad U ðrÞ are computed and the current iterative approximation is updated: c(m+1)(r) = c(m)(r) + k(m) grad U(m)(r). The process returns to the stage 2. The step k(m) of the gradient descent is chosen based on a priori considerations. A more accurate determination of the step of the steepest descent requires extra iterations, which would increase the computing time by a factor of two or more. If the residual U(m+1) at the next iteration is greater than U(m), the step k(m) is decreased by a factor of 1.5. Therefore, the condition U(m+1) 6 U(m) is always fulfilled. The inverse problems of wave tomography considered here are ill-posed. For the approximate solution of the inverse problem with the input data error d, we can take an element c(m(d)) of the iterative sequence c(m), where m(d) is the first number such that inequality U(m(d)) 6 bd2 is fulfilled. Here, b is a fixed number greater than 1. In our simulations, we set b equal to 1.1. According to the theory of iterative regularisation [31,32], such iterative algorithms have regularising properties in the sense that the constructed approximate solution tends to the exact solution of the problem as the error d tends to zero.
Fig. 3. The initial pulse: (a) waveform and (b) spectrum of the pulse for k = 5 mm.
140
A.V. Goncharsky et al. / Ultrasonics 67 (2016) 136–150
The root mean squared (RMS) error d(l) is computed for (l)-th source using the following formula:
d2ðlÞ ¼
2 1 X t ~ tijk ; u u qT i;j;k2C ijk t¼0T
where utijk are the exact values of u(r, t) at the detectors, which are ~ tijk is the perturbed data determined by solving the direct problem; u used as experimental data for solving the inverse problem, q is the number of detectors, and T is the number of time steps. The RMS error d for the entire data set is computed using the formula P d2 ¼ 1L l d2ðlÞ , where L is the number of sources. The proposed algorithms can be efficiently parallelised on computer clusters. More than 90% of the time is spent for computing the gradient. If the number of computing nodes in the cluster is equal to that of the ultrasound sources, then, given that P grad UðmÞ ðrÞ ¼ l grad UðlÞ ðrÞ, the computation of the gradient can be performed in parallel for all the sources. The time-domain method is well-suited for GPU computing. Finite-difference methods operate on large 2D and 3D data sets. According to our analysis, for 2D problems of ultrasonic tomography a single graphics processing unit (GPU) reduces the computing time approximately by a factor of 30 compared to a personal computer. Using 30 GPU units to compute the gradient in parallel for 30 sources provides a speedup factor of 30 30 = 900. We performed the computations on the GPU cluster that is part of the Lomonosov Moscow State University supercomputer complex. The computing time was about 5 min for the 2D case and about 2 h for the 3D case. We used one GPU NVidia Tesla X2070 for each ultrasound source. 3. Numerical simulations of two-dimensional problems of ultrasonic tomography 3.1. Two-dimensional full-range tomography scheme For computations with full-range data we used the acoustic pressure data U measured over the entire boundary of the computational domain to solve the inverse problem. The size of the computational domain was 200 200 mm and the size of the finite difference grid was 640 640 pixels. The 2D simulations were performed for the centre wavelength of the initial pulse k = 10 mm and k = 5 mm, so the total number of detectors equals 316 (80 on each side). The range of sound speed v(r) inside the phantom was 1.43–1.6 km s1, and the ambient sound speed, v0 = 1.5 km s1. The signal recording time and the sampling frequency were equal to T = 250 ls and 12.8 MHz, respectively. The simplest finitedifference scheme used in our implementation requires the time step to be sufficiently small to ensure stability (the Courant condition). This can possibly be avoided by using implicit schemes or the finite elements method. An advantage of simple finite-difference schemes is that they allow easy parallelization of the numerical algorithms. Fig. 4 shows the scheme of our 2D numerical simulations. Simulations shown in Figs. 6–13 were performed for k = 10 mm and the inter-detector distance of k/4 = 2.5 mm. We obtained the input data for the inverse problem by solving the direct problem using the pre-defined sound speed phantom v(r) shown in Fig. 5. Fig. 6 shows the reconstructed images in the case of exact input data, using 8 ultrasound sources. These images have been obtained at the 400th iteration. Fig. 7 shows enlarged parts of the reconstructed images obtained using exact input data. It is evident from these figures that inclusions as small as 2 mm in size and just 4 mm apart can be easily discerned. Increasing the number of sources improves reconstruction quality.
Fig. 4. Scheme of the 2D simulations with full-range data.
The method proposed in this paper is a gradient technique, for which U(m+1) 6 U(m). Fig. 8 shows the history of the residual U(m) as a function of the iteration number m for simulations with exact input data. The convergence of the iterative process at initial iterations can possibly be improved by using randomly combined groups of sources instead of single sources [33]. We solved the inverse problems both with exact input data U(c, t) and with randomly perturbed data U(c, t) + dE(c, t), where E (c, t) has a standard normal distribution. The error level d was set to 1.25% and 5% of the peak magnitude on the detectors maxc2C ðmaxt Uðc; tÞ mint Uðc; tÞÞ 0.4. Fig. 9 shows plots of the perturbed data U(c, t)+dE(c, t) with 1.25% (a) and 5% (b) error level at one of the detectors c. The solid line shows the exact value U(c, t) obtained by solving the direct problem. Fig. 10(a) shows the image reconstructed from noisy data with 1.25% error. The mean squared error d2 = 2.5 105. The residual functional value after 48 iterations becomes U(48) = 2.55 105.
Fig. 5. The sound speed phantom used to solve the direct problem.
A.V. Goncharsky et al. / Ultrasonics 67 (2016) 136–150
Fig. 6. The 2D image reconstructed using exact input data and eight sources.
Fig. 7. Enlarged parts of the reconstructed images: (a) the exact sound-speed image and (b) the image reconstructed using eight sources and k = 10 mm.
141
5% error. Simulations showed that the accuracy of the sound speed reconstruction can be improved by increasing the number of sources. Fig. 12 plots the velocity profile along the A–A line (shown in Fig. 6) for the image reconstructed using 8 sources and exact input data (a) and noisy data with 1.25% error (b). The dashed lines represent the exact sound speed. As is evident from these figures, our technique allows not only the boundaries of the inclusions but also the velocity structure to be reconstructed fairly well, as long as the input data error is sufficiently small. Fig. 13 shows the images reconstructed with different spacing between the detectors. In Fig. 13(b) and (c) the spacing between the detectors is equal to 2.5 and 5 mm, respectively. Fig. 13(a) shows the sound speed phantom containing a small-scale structure. As is evident from the figure, the image reconstructed with detectors spaced 2.5 mm apart is better than the image reconstructed with the detectors spaced 5 mm apart. Both images are reconstructed from exact input data obtained by solving the direct problem. As is evident from a comparison of Fig. 13(b) and (c), the optimum spacing between the detectors should be between k/4 and k/2. Basically, the smaller is the spacing between the detectors, the better is the quality of the reconstructed image, however, there is little point in placing detectors closer than k/4 from each other because of the increased complexity and cost of the hardware involved. The number of sources should be increased with reducing the wavelength. We performed simulations for k = 5 mm with 16 sources of ultrasonic emission and exact data at the detectors. Fig. 14(a) shows the solution of our simulated problem for k = 5 mm with detectors spaced k/2 = 2.5 mm apart. Fig. 14(b) shows the velocity structure reconstructed with detectors spaced 5 mm apart. As is evident from a comparison of Fig. 14 (a) and (b), the quality of the image reconstructed with interdetector spacing of about k is substantially inferior to that obtained with detectors spaced about k/2 apart. We found similar results for simulated problems with input data with nonzero measurement errors. We used NVidia GeForce GTX 460 and NVidia Tesla X2070 GPUs for our 2D numerical simulations. According to our results, GTX 460 and NVidia Tesla X2070 GPUs reduce the computing time by a factor of 20 and 40, respectively, compared to computations on a PC without using a GPU [29]. A GPU-equipped supercomputer is needed for solving 3D problems of ultrasonic tomography.
3.2. Two-dimensional transmission tomography scheme
Fig. 8. Residual U(m) as a function of the iteration number m.
Hence, one can adopt the image obtained at the 48th iteration as an approximate solution of the inverse problem. It takes approximately 5 min to solve a 2D inverse problem with eight ultrasound sources using eight GPU units in parallel. Fig. 10(b) shows the image reconstructed from noisy data with 5% error. The mean squared error and the residual functional values in this case are equal to d2 = 4 104 and U(12) = 4.02 104, respectively. Hence c(12) can be adopted as an approximate solution. For comparison, we show enlarged parts of the reconstructed images in Fig. 11. The smallest inclusions are 2 mm in size and are located 4 mm apart. These inclusions are visible even in the case of
Let us now consider the inverse problem of transmission tomography. Fig. 15 shows the arrangement of the sources and detectors in the transmission tomography scheme. Eight sources are located at the boundary of the computational domain. The detectors corresponding to each source position S are located on the opposite side of the object. We draw the line S O S0 passing through a source S and through the centre O of the computational domain. For each source position S, only the detectors D for which the angle D O S0 is no greater than a are used in the transmission tomography scheme. We have performed the numerical simulations for a = 60°. In our transmission tomography simulations, we used the sound speed phantom shown in Fig. 5 to solve the direct problem and to obtain the input data for the inverse problem. Fig. 16 shows the images reconstructed for a = 60° using exact input data (a) and noisy data with 1.25% error (b). It is evident from a comparison of Fig. 16 and Figs. 10–12 that the smallest
142
A.V. Goncharsky et al. / Ultrasonics 67 (2016) 136–150
Fig. 9. Exact and noisy data at a fixed detector position c: (a) for 1.25% error and (b) for 5% error.
Fig. 10. The 2D images reconstructed using 8 sources and noisy data: (a) with 1.25% error and (b) with 5% error.
Fig. 11. Enlarged parts of the reconstructed images: (a) using exact data, (b) 1.25% error and (c) 5% error.
inclusions, which are 2 mm in size, are practically indiscernible in the transmission tomography image, and hence, the full-range tomography scheme yields more accurate results than the transmission tomography scheme. Fig. 17 plots the velocity profiles along the A–A line reconstructed using the transmission tomography scheme with a = 60°, for exact data (a) and for noisy data with 1.25% error (b). The dashed line corresponds to the exact sound speed image used to solve the direct problem. It is evident that in the case of transmission tomography, the sound speed is reconstructed with substantial errors and much less
accurately than in the case of full-range tomography. Inaccuracy of the input data decreases the reconstruction quality even further.
3.3. Two-dimensional reconstruction using backscatter data Reflection-based diagnostic systems are widely used both in industrial diagnostics and in medicine. The reflection problem arises, for example, when a transmitted wave cannot be recorded because of the strong attenuation of high-frequency ultrasound in the medium or when it is impossible to adequately arrange
A.V. Goncharsky et al. / Ultrasonics 67 (2016) 136–150
143
Fig. 12. The A–A velocity profile reconstructed using full-range tomography scheme with 8 sources: (a) for exact input data obtained from the direct problem and (b) for noisy data with 1.25% error.
Fig. 13. The images containing a fine structure: (a) exact sound speed image, (b) reconstructed image, detectors spaced 2.5 mm apart and (c) reconstructed image, detectors spaced 5 mm apart.
Fig. 14. The 2D images reconstructed for k = 5 mm: (a) detectors spaced at 2.5 mm and (b) detectors spaced at 5 mm.
the detectors to record the transmitted waves and only the backscatter data can be obtained (e.g., in seismic studies [34,35]). One can try to reconstruct the velocity structure in a tomographic scheme using only the backscatter data, which can be
obtained with common ultrasonic facilities. We show the scheme of our backscatter tomography simulation in Fig. 18. The detectors corresponding to each source position S are located on the same side of the object as the source. For each source S, only the
144
A.V. Goncharsky et al. / Ultrasonics 67 (2016) 136–150
the backscatter data, one can attempt to reconstruct only the boundaries of irregularities, and even such reconstructions are fraught with substantial errors. Experimental studies [35,36] comply with these results. 4. Numerical simulations of 3D problems of ultrasonic tomography
Fig. 15. Scheme of the 2D transmission tomography simulations.
detectors R for which the angle R O S does not exceed a pre-defined value a are used. We used 32 source positions S and detectors spaced 0.3 mm apart, i.e., at a = 45° a total of 640 detectors are used for each source. Fig. 19 shows the image reconstructed using only the backscatter data with a = 45°. We performed our numerical simulations with exact data obtained by solving the direct problem. The range of sound speed v(r) inside the phantom was 1.43–1.6 km s1. Fig. 19 shows an image that could be obtained using the backscatter data under ideal conditions with no measurement error. The quality of the reconstructed image depends on the number of detectors and sources. As is evident from Fig. 19, even a configuration with 32 sources and many detectors (about 2560) does not allow high quality of the reconstructed image to be achieved. Fig. 20 plots the reconstructed velocity profile along the A–A line. As is evident from Figs. 19 and 20, the velocity structure inside the medium cannot be reconstructed using the backscatter data alone. The accuracy of the sound speed reconstruction in this case even with exact input data (the plot in Fig. 20) is comparably worse than the result obtained in the case of full-range tomography (Fig. 12) and transmission tomography (the plots in Fig. 17). Using
In the previous section, we studied the possibility of solving inverse problems of ultrasonic tomography using the layer-bylayer scheme. The actual problem of ultrasonic tomography is three-dimensional with the sources and detectors located on a surface surrounding a three-dimensional object. Inverse problems in three-dimensional representation are physically more similar to reality. Singling out a two-dimensional layer in X-ray tomography is justified because the X-ray refraction index of any material is close to unity. From the physical point of view, singling out a 2D layer in ultrasonic tomography problems is a considerable simplification of the model because of the refraction, reflection and scattering that are actually present in three dimensions. In this section, we attempt to assess the potentialities of different tomography schemes in the explicitly three-dimensional version. Three-dimensional inverse problems of ultrasonic tomography are very computationally expensive. We have chosen a smaller 3D object to reduce the computation time. The cube-shaped computational domain was 100 100 100 mm in size and the finite difference grid was 320 320 320 pixels. Fig. 21 illustrates the layout of the 3D experiment. In our numerical simulations we used 24 ultrasound sources S located in two horizontal planes Z = h0 = 0.33h and Z = h1 = 0.66h, where h = 100 mm is the height of the computational domain. The wavelength k was equal to 5 mm, and the sound speed inside the simulated phantom varied from 1.43 to 1.6 km s1, whereas that in the ambient medium was 1.5 km s1. The signal recording time and sampling frequency were equal to T = 125 ls and 12.8 MHz, respectively. 4.1. Three-dimensional full-range tomography simulations In the full-range 3D tomography simulation, the detectors are located over the entire surface of the cube-shaped computational domain and are spaced 2.5 mm apart, so that 40 40 = 1600 detectors are located on each face of the cube for a total of 8976
Fig. 16. The 2D images reconstructed via the transmission tomography scheme with a = 60°: (a) using exact input data and (b) using noisy data with 1.25% error.
A.V. Goncharsky et al. / Ultrasonics 67 (2016) 136–150
145
Fig. 17. The A–A velocity profiles reconstructed using the transmission tomography scheme with a = 60°: (a) for exact input data and (b) for noisy data with 1.25% error.
Fig. 20. Plot of the A–A velocity profile reconstructed using the backscatter data.
Fig. 18. Scheme of 2D backscatter tomography simulation.
Fig. 21. The layout of the 3D simulation.
Fig. 19. The 2D image reconstructed using only the backscatter data.
detectors. Fig. 22 shows the cross sections of the reconstructed three-dimensional velocity structure. Fig. 23 shows the images reconstructed from noisy data with 1.25% error. Fig. 24 shows enlarged parts of the images reconstructed from exact data, noisy data with 1.25% error and noisy data with 5% error. It follows from Fig. 24 that an ultrasonic examination with a wavelength of k = 5 mm allows inclusions as small as 1 mm that are spaced 2 mm apart to be resolved if the input data error is
146
A.V. Goncharsky et al. / Ultrasonics 67 (2016) 136–150
Fig. 22. Cross sections of the 3D velocity structure reconstructed via full-range tomography scheme using k = 5 mm and exact input data.
Fig. 23. Cross sections of the 3D velocity structure reconstructed via full-range tomography scheme using k = 5 mm and noisy data with 1.25% error.
A.V. Goncharsky et al. / Ultrasonics 67 (2016) 136–150
147
Fig. 24. Cross sections of the 3D velocity structure reconstructed via full-range tomography scheme using k = 5 mm: (a) exact input data, (b) noisy data with 1.25% error and (c) noisy data with 5% error.
Fig. 25. Cross sections of the 3D velocity structure reconstructed via transmission tomography scheme using k = 5 mm and exact input data.
sufficiently small. These results are based on numerical simulations and determine some of the ultimate capabilities of physical experiments. The choice of the 2 mm resolution is determined by requirements of medical professionals, who need to diagnose breast cancer at early disease stages when the tumor has the size of 2–4 mm (no greater than 4 mm). Currently, prototypes exist with declared resolution of about 2 mm (for transmissiontomography scheme) when using acoustic sources with central frequencies of 1.5–2.5 MHz. The results obtained in this study show that ultrasonic tomography devices can be developed to provide a resolution of about 2–3 mm when operating at much lower frequencies of 0.3–0.5 MHz. The results of numerical simulations suggest that the use of lower frequencies has a number of advantages due primarily to much weaker attenuation. Another advantage of the use of low frequencies is the possibility of high-precision registration of the phase of acoustic waves.
4.2. Three-dimensional transmission tomography simulations In our transmission tomography scheme, the detectors for each source S are located only on the face of the cube opposite to the source, for a total of 40 40 = 1600 detectors for each source. The wavelength is k = 5 mm, and the sources are arranged as shown in Fig. 21. Fig. 25 shows the cross sections of the 3D image reconstructed via a transmission tomography scheme. As is evident from Figs. 25 and 22, the quality of the images reconstructed via transmission tomography is inferior to that of
the images reconstructed via full-range tomography. The small details become indiscernible. As for the reflection (backscatter) scheme, in Section 3.3 we showed that the velocity structure c(r) cannot be reconstructed using the backscatter data alone. 4.3. Three-dimensional limited data tomography simulations One of the typical inverse problems of wave tomography is the limited data tomography problem, where the object cannot be sounded from all sides. In the case of a tomographic examination aimed at breast cancer diagnosis, neither the sources nor the detectors can be placed on the thorax side. However, a reasonably chosen arrangement of the sources and detectors allows the velocity structure to be reconstructed rather accurately. Consider now a simulated problem of limited data tomography, where the sources are located, like in the above case, in the Z = h0 = 33 mm and Z = h1 = 66 mm planes, as shown in Fig. 21, and the detectors are located on all the faces of the cube-shaped computational domain except the top face Z = h. Fig. 26 shows the cross sections of the 3D velocity structure reconstructed using the limited data tomography scheme and the exact input data obtained by solving the direct problem. These cross sections correspond to the Z = 95 mm and Z = 85 mm planes. As is evident from Fig. 26, although the exact input data are used for the inverse problem, the quality of the reconstructed image significantly decreases near the Z = h plane, where no detectors are present. However, at some distance from the Z = h plane (e.g., at Z < 0.85h), the image reconstruction quality is fairly good.
148
A.V. Goncharsky et al. / Ultrasonics 67 (2016) 136–150
Fig. 26. Cross sections of the 3D velocity structure reconstructed using the limited data tomography scheme and exact input data: (a) Z = 95 mm and (b) Z = 85 mm. The sources are arranged as shown in Fig. 21.
Fig. 27. Scheme of the 3D limited data tomography simulations.
The reconstruction quality in this problem can be improved by placing more sources near the Z = h plane, e.g., by increasing the height h1 of one of the planes where the sources are located to 0.95h, as shown in Fig. 27. Fig. 28 shows the reconstructed images corresponding to the arrangement where the sources are positioned in the h0 = 0.33h and h1 = 0.95h planes. Exact input data are used for the reconstruction. Fig. 29 shows the same cross sections reconstructed using noisy input data with 1.25% error. As is evident from Figs. 28 and 29, the reconstructed cross sections near the top boundary have practically the same quality as the images reconstructed using full-range tomography scheme, shown in Figs. 22 and 23. 5. Discussion and conclusion The problem of the analysis of various ultrasonic tomography schemes was addressed in a number of studies [6,7,23]. Most of
them discuss the reconstruction of 3D objects in terms of the layer-by-layer approach. In this paper we analyse the inverse problem of ultrasonic tomography not only in the layer-by-layer case, but also in the three-dimensional formulation where the threedimensional function is reconstructed directly from the input data. The inverse problem is treated as a three-dimensional coefficient inverse problem for a second-order differential equation. The algorithms for solving the inverse problem are based on the possibility of directly computing the gradient of the residual functional. We show that the full-range tomography scheme should be preferred both in the case of the two- and three-dimensional approach. Three-dimensional nonlinear problems of ultrasonic diagnosis involving the reconstruction of a 3D velocity structure are too challenging to be solved using commodity PC the formulation considered in this paper. Iterative gradient methods for solving inverse problems of wave tomography are easily parallelizable, allowing these problems to be solved on modern supercomputers. In three-dimensional inverse problems of ultrasonic breast cancer diagnosis the inverse problem with incomplete data must be solved, because there are neither sources nor detectors present on the patient’s thorax side. We show in this paper that a special setup of the sources allows bona fide reconstruction of the threedimensional velocity structure. Detailed measurements of both the amplitude and phase of the wave allow the velocity structure to be reliably reconstructed with a resolution of about 2–3 mm using low-frequency sources operating in the 300–500 kHz range. The accuracy with which the detectors measure the wave field is an important factor and it directly depends on the detector spacing. We show in this paper that a configuration with few sources operating in the frequency interval 300–500 kHz allows the velocity structure to be reliably reconstructed with detectors spaced k/4 to k/2 apart. At higher frequencies performing detailed phase measurements with transducers with sizes no greater than 1 mm appears to be a very challenging task primarily because of the large number of the detectors involved. The inverse problem analysed in this paper is ill-posed. The methods of solving ill-posed problems developed in late last century allow one to address the fundamental problems that arise when interpreting the data of physical experiments. Modern methods of the solution of ill-posed problems allow an approximate solution of the inverse problem to be constructed that converges to the exact solution of the inverse problem as the experimental
A.V. Goncharsky et al. / Ultrasonics 67 (2016) 136–150
149
Fig. 28. Cross sections of the 3D velocity structure reconstructed using the limited data tomography scheme and exact input data: (a) Z = 95 mm and (b) Z = 90 mm. The sources are arranged as shown in Fig. 26.
Fig. 29. Cross sections of the 3D velocity structure reconstructed using the limited data tomography scheme and noisy input data with 1.25% error: (a) Z = 95 mm and (b) Z = 90 mm. The sources are arranged as shown in Fig. 26.
error goes to zero. From the point of view of this theory, arbitrarily high resolution can be achieved if the data errors are sufficiently small [31,37,38]. In reality, the resolution of a tomography device depends on both the errors of the input data and the frequency of sounding pulses. In this study we choose low frequencies for sounding pulses in order to reduce attenuation and minimise the input data error. The quality of reconstructed images directly depends on the degree of detail with which the phase of acoustic waves can be recorded. The optimum setup is to place the detectors with a spacing of about k/2, which for frequencies 300–500 kHz is equal to about 2 mm. The currently existing prototype ultrasonic tomography devices operate at frequencies of about 1.5 MHz making recording ultrasonic waves with a discretization of 0.5 mm too challenging a task. According to the results of our numerical simulations, in the chosen frequency range of 300–500 kHz a resolution of about 2–3 mm is quite achievable for the velocity contrast of the object no greater than 10%. Tomographic reconstruction of ultrasonic
images with such a resolution can already be used for differential diagnosis of breast cancer [39]. It goes without saying that the adequacy of the model employed to the real physical experiment is of great importance. The answer to this question can be obtained only experimentally. Three-dimensional problems of ultrasonic diagnostics cannot be solved without using supercomputers. The parameters of general-purpose supercomputers, such as their size, cost and power consumption, which amounts to several megawatts, make them impossible to be incorporated into tomography facilities. In this paper we propose an actual solution to this problem. The approach to the low-frequency ultrasonic tomography considered in the paper involves the use of a relatively small number of ultrasound sources and many detectors. Optimum parallelization can be achieved if the number of sources is equal to that of the computing nodes of the GPU cluster. With 30 ultrasound sources processed in parallel using currently available GPU units, the speedup is a factor of about one thousand compared to a PC. NVidia announced the
150
A.V. Goncharsky et al. / Ultrasonics 67 (2016) 136–150
forthcoming release of NVidia Volta graphics processors, which are to become available on the market in 2015–2016 and which will reduce the computing time down to about 30 min for the threedimensional problem. Such a computing time is quite acceptable for practical use in ultrasonic tomography devices. Acknowledgements We are pleased to thank Prof. Vl. V. Voevodin, a corresponding member of the Russian Academy of Sciences, for numerous discussions of the work. We are also grateful to Academician V.A. Kubyshkin, Director of The Vishnevsky Institute of Surgery of the Russian Academy of Medical Sciences, for advice concerning the medical aspects of the diagnosis of oncological diseases. This work was supported by the Russian Foundation for Basic Research (Project No. 14-07-00078-a). References [1] J.A. Jensen, Ultrasound imaging and its modeling, Top. Appl. Phys. 84 (2002) 135–166. [2] M.F. Schiffner, G. Schmitz, Plane wave pulse-echo ultrasound diffraction tomography with a fixed linear transducer array, Acoust. Imag. 31 (2012) 19– 30, http://dx.doi.org/10.1007/978-94-007-2619-2_3. [3] A. Sanpanich, P. Greesuradej, S. Aootaphao, C. Pintavirooj, M. Sangworasil, P. Tosranon, 3D ultrasonic reflection tomography with matrix linear array transducer, ISBME 44 (2008) 351–355. [4] C.R. Hill, J.C. Bamber, Gail R. ter Haar, Physical Principles of Medical Ultrasonics, second ed., Wiley, 2005, http://dx.doi.org/10.1002/0470093978. [5] H. Gemmeke, A. Menshikov, D. Tchernikovski, L. Berger, G. Göbel, M. Birk, M. Zapf, N. Ruiter, Hardware setup for the next generation of 3D ultrasound computer tomography, NSS/MIC IEEE (2010) 2449–2454, http://dx.doi.org/ 10.1109/NSSMIC.2010.5874228. [6] F. Natterer, Possibilities and limitations of time domain wave equation imaging, Contemp. Math. 559 (2011) 151–162. [7] R.J. Lavarello, M.L. Oelze, Tomographic reconstruction of three-dimensional volumes using the distorted Born iterative method, IEEE Trans. Med. Imag. 28 (10) (2009) 1643–1653, http://dx.doi.org/10.1109/TMI.2009.2026274. [8] K.J. Opielinski, T. Gudra, Three-dimensional reconstruction of biological objects’ internal structure heterogeneity from the set of ultrasonic tomograms, Ultrasonics 42 (2004) 705–711. [9] L. Huang, Y. Quan, Sound-speed tomography using first-arrival transmission ultrasound for a ring array, Proc. SPIE, Med. Imag.: Ultrason. Imag. Signal Process. 6513 (2007) 651306, http://dx.doi.org/10.1117/12.709647. [10] N. Duric, P. Littrup, C. Li, O. Roy, S. Schmidt, R. Janer, X. Cheng, J. Goll, O. Rama, L. Bey-Knight, W. Greenway, Breast ultrasound tomography: bridging the gap to clinical practice, Proc. SPIE, Med. Imag.: Ultrason. Imag., Tomogr., Ther. 8320 (2012) 83200O, http://dx.doi.org/10.1117/12.910988. [11] S. Li, K. Mueller, M. Jackowski, D. Dione, L. Staib, Physical-Space RefractionCorrected Transmission Ultrasound Computed Tomography Made Computationally Practical, MICCAI, Part II, LNCS, vol. 5242, 2008, pp. 280–288. [12] R. Jirik, I. Peterlík, N. Ruiter, J. Fousek, R. Dapp, M. Zapf, J. Jan, Sound-speed image reconstruction in sparse-aperture 3-D ultrasound transmission tomography, IEEE Trans. Ultrason. Ferroelectr. Freq. Control 59 (2012) 254– 264, http://dx.doi.org/10.1109/TUFFC.2012.2185. [13] K.J. Opielinski, T. Gudra, Ultrasound transmission tomography image distortions caused by the refraction effect, Ultrasonics 38 (2000) 424–429, http://dx.doi.org/10.1016/S0041-624X(99)00137-7. [14] O. Roy, I. Jovanovic, A. Hormati, R. Parhizkar, M. Vetterli, Sound speed estimation using wave-based ultrasound tomography: theory and GPU implementation, Proc. SPIE Med. Imag.: Ultrason. Imag., Tomogr., Ther. 7629 (2010) 762909, http://dx.doi.org/10.1117/12.844691. [15] E. Kretzek, N.V. Ruiter, GPU based 3D SAFT reconstruction including phase aberration, Proc. SPIE 9040 (2014) 90400W, http://dx.doi.org/10.1117/ 12.2042669.
[16] A.V. Goncharsky, S.Y. Romanov, S.Y. Seryozhnikov, Inverse problems of 3D ultrasonic tomography with complete and incomplete range data, Wave Motion 51 (2014) 389–404, http://dx.doi.org/10.1016/j.wavemoti.2013. 10.001. [17] S. Mensah, R. Ferriere, Diffraction tomography: a geometrical distortion free procedure, Ultrasonics 42 (2004) 677–682, http://dx.doi.org/10.1016/j. ultras.2003.11.012. [18] R.G. Pratt, L. Huang, L.N. Duric, P. Littrup, Sound-speed and attenuation imaging of breast tissue using waveform tomography of transmission ultrasound data, Proc. SPIE, Med. Imag.: Phys. Med. Imag. 6510 (2007) 65104S, http://dx.doi.org/10.1117/12.708789. [19] F. Natterer, F. Wubbeling, A propagation–backpropagation method for ultrasound tomography, Inverse Probl. 11 (1995) 1225–1232, http://dx.doi. org/10.1088/0266-5611/11/6/007. [20] A.H. Rohde, M. Veidt, L.R.F. Rose, J. Homer, A computer simulation study of imaging flexural inhomogeneities using plate-wave diffraction tomography, Ultrasonics 48 (2008) 6–15, http://dx.doi.org/10.1016/j.ultras.2007.09.002. [21] F. Simonetti, L. Huang, N. Duric, O. Rama, Imaging beyond the Born approximation: an experimental investigation with an ultrasonic ring array, Phys. Rev. E 76 (2007) 036601, http://dx.doi.org/10.1103/PhysRevE.76. 036601. [22] F. Natterer, Ultrasonic Image Reconstruction via Plane Wave Stacking: Possibilities and Limitations of Time Domain Wave Equation Imaging, Technical Report, Fachbereich Mathematik und Informatik, Universitat Munster, 2005. [23] J. Wiskin, D. Borup, M. Andre, S. Johnson, J. Greenleaf, Y. Parisky, J. Klock, Three-dimensional nonlinear inverse scattering: quantitative transmission algorithms, refraction corrected reflection, scanner design, and clinical results, J. Acoust. Soc. Am. 133 (2013) 3229, http://dx.doi.org/10.1121/1.4805138. [24] R.J. Lavarello, M.L. Oelze, Density imaging using inverse scattering, J. Acoust. Soc. Am. 125 (2) (2009) 793–802, http://dx.doi.org/10.1121/1.3050249. [25] A.V. Goncharsky, S.Y. Romanov, Inverse problems of ultrasound tomography in models with attenuation, Phys. Med. Biol. 59 (2014) 1979–2004, http://dx.doi. org/10.1088/0031-9155/59/8/1979. [26] J. Wiskin, D. Borup, S. Johnson, M. Berggren, T. Abbott, R. Hanover, Full wave, non-linear, inverse scattering, Acoust. Imag. 28 (2007) 183–194. [27] A.V. Goncharskii, S.Y. Romanov, Two approaches to the solution of coefficient inverse problems for wave equations, Comput. Math. Math. Phys. 52 (2012) 245–251, http://dx.doi.org/10.1134/S0965542512020078. [28] L. Beilina, M.V. Klibanov, Approximate Global Convergence and Adaptivity for Coefficient Inverse Problems, Springer, New York, 2012, http://dx.doi.org/ 10.1007/978-1-4419-7805-9. [29] A.V. Goncharsky, S.Y. Romanov, Supercomputer technologies in inverse problems of ultrasound tomography, Inverse Probl. 29 (2013) 075004, http://dx.doi.org/10.1088/0266-5611/29/7/075004. [30] B. Engquist, A. Majda, Absorbing boundary conditions for the numerical simulation of waves, Math. Comput. 31 (1977) 629. [31] A.B. Bakushinsky, A.V. Goncharsky, Ill-posed Problems: Theory and Applications, Kluwer Acad. Publ., Dordrect, 1994. [32] A.B. Bakushinsky, M.Y. Kokurin, Algorithmic Analysis of Irregular Operator Equations, LENAND, Moscow, 2012 (in Russian). [33] K. Wang, T. Matthews, F. Anis, C. Li, N. Duric, M. Anastasio, Waveform inversion with source encoding for breast sound speed reconstruction in ultrasound computed tomography, IEEE Trans. Ultrason., Ferroelectr. Freq. Contr. 62 (3) (2014), http://dx.doi.org/10.1109/TUFFC.2014.006788. 12/2014. [34] G.E. Murphy, B.P. Amoco, S.H. Gray, Manual seismic reflection tomography, Geophysics 64 (1999) 1546–1552, http://dx.doi.org/10.1190/1.1444658. [35] P. Lasaygues, J.P. Lefebvre, Cancellous and cortical bone imaging by reflected tomography, Ultrason. Imag. 23 (2001) 55–70, http://dx.doi.org/10.1177/ 016173460102300104. [36] A. Kak, M. Slaney, Principles of Computerized Tomographic Imaging, SIAM, 2001, http://dx.doi.org/10.1137/1.9780898719277. [37] A.N. Tikhonov, A.V. Goncharsky, V.V. Stepanov, A.G. Yagola, Numerical Methods for the Solution of Ill-posed Problems, Kluwer Acad. Publ., Dordrecht, Boston, London, 1995, ISBN 0-7923-3583-X. [38] A.N. Tikhonov, The solution of ill-posed problems and the regularization method, DAN SSSR 151 (3) (1963) 501–504 (in Russian). [39] J. Wiskin, D.T. Borup, S.A. Johnson, M. Berggren, Non-linear inverse scattering: high resolution quantitative breast tissue tomography, J. Acoust. Soc. Am. 131 (5) (2012) 3802–3813, http://dx.doi.org/10.1121/1.3699240.