An algorithm for 3D ultrasound tomography

Frank Natterer

Institut für Numerische und instrumentelle Mathematik, Westfälische Wilhelms-Universität, Einsteinstrasse 62, D-48149 Münster, Germany e-mail: frank natterer@math.uni-muenster.de

1 Introduction

We consider the inverse problem for the 3D Helmholtz equation

Δu^{j} + k^{2}(1 + f)u^{j} = 0

(1.1)

u^{j} = u_{i}^{j} + v^{j} ,

where u_{i}^{j}, j = 1, …, p are the incoming waves, v^{j} satisfies the Sommerfeld radiation condition and the function f vanishes outside the ball of radius ρ. We want to recover f numerically from knowing u^{j} = g^{j} on the sphere of radius ρ for j = 1, …, p and a fixed frequency k.

This is a model for ultrasonic tomography [6]. However we point out that in a real ultrasonic scanner the irradiating waves are no longer plane waves but standing waves in a finite container.

We start with a short survey on the extensive literature on numerical methods. With the exception of methods which use the Born or Rytov approximation [4], it seems that the only method which actually has been tested numerically in 3D is the Newton method combined with a finite Fourier expansion of f [5]. The other methods have been used in 2D only, even though an extension to 3D is possible in principle. The dual space method [2] reduces the problem from the whole space IR^{3} with the far field (which we do not use) as data to an overposed boundary value problem in a finite volume which in turn is solved by optimization.

Of course we can always try to compute the Born series [3]. For this purpose we write (1.1) as an integral equation

u^{j}(x) = u_{i}^{j}(x) − k^{2} ∫ G(x, y)f(y)u^{j}(y)dy .  (1.2)

Here, G is the Green’s function of Δ + k^{2} with the radiation condition at ∞. The Born series (f_{ℓ}) is now obtained by solving

g^{j}(x) = u_{i}^{j}(x) − k^{2} ∫ G(x, y)f_{ℓ}(y)u_{ℓ}^{j}(y)dy , |x| = ρ , j = 1, …, p

for f_{ℓ}, where u_{0}^{j} = u_{i}^{j} and

Δu_{ℓ+1}^{j} + k^{2}(1 + f_{ℓ})u_{ℓ+1}^{j} = 0

with u_{ℓ}^{j} = u_{i}^{j} + a function satisfying the radiation condition.

The generalized SOR-method of [7] also starts out from the integral equation (1.2). Writing for the integral operator in (1.2) simply G, this method minimizes the functional

∑_{j=1}^{p} {∥g^{j} + u_{i}^{j} + k^{2}Gfu^{j}∥_{L_{2}(|x|=ρ)}^{2} + γ∥u^{j} − u_{i}^{j} + k^{2}Gfu^{j}∥_{L_{2}(|x|<ρ)}^{2}} (1.3)

with some weight factor γ. The minimization is done by

u_{ℓ+1}^{j} = u_{ℓ}^{j} + α_{ℓ}r_{ℓ}^{j} , f_{ℓ+1} = f_{ℓ} + β_{ℓ}d_{ℓ}

where r_{ℓ}^{j}, d_{ℓ} are update directions and α_{ℓ}, β_{ℓ} are chosen so as to minimize (1.3).

A non iterative method has been suggested in [15]. With G_{f} the Green’s function of Δ + k^{2}(1 + f), we can rewrite (1.1) as

v^{j}(x) = −k^{2} ∫ G_{f}(x, y)f(x)u_{i}^{j}(y)dy .  (1.4)

Now form a linear combination of the incoming waves such that the resulting field peaks at z, i.e.

∑_{j=1}^{p} α_{j}(z)u_{i}^{j}(x) ∼ δ(x − z) .

Then, from (1.4),

− k^{2}G_{f}(x, z)f(z) ∼ ∑_{j=1}^{p} α_{j}(z)v^{j}(x)  (1.5)

is approximately known on |x| = ρ. From the identity

G_{f}(x, y) − G(x, y) = −k^{2} ∫_{|z|<ρ} G_{f}(x, z)f(z)G(z, y)dy

we get an approximation to G_{f}(x, y) for |x| = ρ which together with (1.5) determines f. The method still has to be tested.

A very efficient code has been given in [6]. It is very similar to ours in that it uses initial value techniques and ignores in the Jacobian entries which correspond to different incoming waves. The initial value technique is based on factoring the Helmholtz equation into a product of two first order differential operators.

The purpose of this note is to extend the method of [14] to 3D and to conduct 3D numerical experiments. We illuminate the object by plane waves u_{i}^{j} = e^{ikx·θ_{j}} where the unit vectors θ_{j} are lying in a plane. We use an initial value technique as in [6], but we do not rely on parabolic approximations [8] of the Helmholtz equation.

2 The initial value method

Let Ω_{j} be the cube circumscribed to the reconstruction region |x| ≤ ρ with edges parallel to the coordinate axes and aligned with the direction θ_{j} = (cos φ_{j}, sin φ_{j}, 0) of the j-th incoming wave. Let Γ_{j}^{±} be the face lying in the plane x · θ_{j} = ±ρ, and let Γ_{j} be the other faces of Ω_{j}. Rather than working with the scattered fields v^{j} we use the scaled scattered fields w^{j} = e^{−ikx·θ_{j}}v^{j} which satisfy

Δw^{j} + 2ikθ_{j} · ∇w^{j} + k^{2}(1 + w^{j})f = 0 in Ω_{j} .  (2.1)

Note that we do not make the parabolic approximation [8], i.e. we do not assume that the second derivative of w^{j} in direction θ_{j} is small compared to kθ_{j} · ∇w^{j}.

In [12] we haved analysed the stability of the initial value problem

w^{j} = g^{j} on Γ_{j} ∪ Γ_{j}^{−} ,  \frac{∂w^{j}}{∂ν} = \frac{∂g^{j}}{∂ν} on Γ_{j}^{−}  (2.2)

for this elliptic differential equation. Let

ŵ^{j}(ξ, s) = \frac{1}{2π} ∫_{x·θ_{j}=s} e^{−ix·ξ}w^{j}(x)dx , ξ ⊥ θ_{j}

be the 2D Fourier transform of w in the plane x · θ_{j} = s. We found that ŵ^{j}(ξ, s) depends in a perfectly stable way on the initial values for w^{j} on Γ_{j}^{−} for all frequencies |ξ| < κ where κ is some number depending essentially on k and, to a minor extent, on f. In ultrasound tomography we can choose κ slightly smaller than k, typically κ = 0.99k.

Thus we may define a nonlinear map R_{j} : L_{2}(|x| < ρ) → L_{2}(Γ_{j}^{+}) by putting

R_{j}(f) = w_{κ}^{j}|_{Γ_{j}^{+}} .  (2.3)

The inverse scattering problem now calls for the solution of the nonlinear system

R_{j}(f) = g_{j} , g_{j} = g^{j}|_{Γ_{j}^{+}} ,  j = 1, …, p .  (2.4)

As in [12] this is done by an ART-type procedure. Starting out from an initial approximation f^{0}, we put f_{0} = f^{0} and for j = 1, …, p

f_{j} = f_{j−1} − ωR′_{j}(f_{j−1})^{∗}C_{j}^{−1}(R_{j}(f_{j−1}) − g_{j}) .

The first approximation f^{1} is then defined to be f_{p}. For C_{j} we simply take the operator γI where γ is chosen such that, in the limit k → ∞, C_{j} ∼ R′_{j}(0)R′_{j}(0)^{∗}, i.e. γ = k^{2}/ρ [10]. The evaluation of R′_{j}(f)^{∗}g for some g ∈ L_{2}(Γ_{j}^{+}) can be done as follows: Solve the initial value problem

Δz + 2ikθ_{j} · ∇z − k^{2}\overline{f}z = 0 in Ω_{j}

(2.5)

z = 0 on Γ_{j} ∪ Γ_{j}^{+} , \frac{∂z}{∂ν} = g on Γ_{j}^{+} .

Then,

R_{j}^{∗}(f)g = k^{2}(1 + \overline{w}^{j})z

where w^{j} is the solution of (2.1) - (2.2). Of course the stability properties of (2.5) are exactly as discussed above.

3 The finite difference method

The numerical solution of the initial value problems (2.1) - (2.2) and (2.5) can efficiently and conveniently be done by a finite difference method. In view of our stability result (see previous section) this suggestions itsself, but it has been used already in [9]. We simply use the usual five point discretization on the grid Ω_{j}^{h} = {hℓθ_{j} + hmθ_{j}^{⊥} + hne_{3} : ℓ, m, n = −q, …, q} where θ_{j}^{⊥} = (− sin φ_{j}, cos φ_{j}, 0) and e_{3} = (0, 0, 1) and h = ρ/q. Then, the finite difference approximation to (2.1) - (2.2) reads (the superscript j is omitted)

w_{ℓ+1,m.n} + w_{ℓ−1,m,n} + w_{ℓ,m+1,n} + w_{ℓ,m−1,n} + w_{ℓ,m,n+1} + w_{ℓ,m,n−1}
−6w_{ℓ,m,n} + iε(w_{ℓ+1,m,n} − w_{ℓ−1,m,n}) + ε^{2}(1 + w_{ℓ,m,n})f_{ℓ,m,n} = 0 ,
ℓ, m, n = −q + 1, …, q − 1 ,  (3.1)

w_{ℓ,m,n} = g_{ℓ,m,n} for |m| = q or |n| = q and for ℓ = −q, −q + 1 .

The boundary conditions in (3.1) assume that the field can be measured everywhere on each of the faces of Ω_{j}. The values of g_{ℓ,m,n} = f^{j}(hℓθ_{j}^{1} + hnθ_{j}^{2}) for ℓ = q − 1 have to be computed from the field in |x| ≥ ρ by numerical differentiation. (3.1) is solved in a recursive way, i.e. if w in known on the levels ℓ, ℓ−1 we compute it for ℓ+1 by (3.1). In order to preserve stability we have to filter out the frequencies greater than κ at each stage of this process. This can be done by doing a 2D FFT on each matrix (w_{ℓ,m,n})_{m,n=−q,...,q} as soon as it is computed, zeroing all the entries with \sqrt{m^{2} + n^{2}} > 2ρκ/π, followed by a 2D inverse FFT. The total operator count for solving (2.1) - (2.2) once is O(q^{3} log q). For details see [16]. (2.5) is solved exactly in the same way.

4 Numerical experiments

In a first test we checked the initial value method for accuracy. We solved the forward problem (1.1) for the function

f(x) = \{ \begin{matrix}0.2 & , & |x| ≤ 0.5 \\ 0 & , & otherwise \end{matrix}  (4.1)

and k = 50 analytically and compared the exact solution with the approximate solution of the initial value method for h = 32. We found satisfactory agreement.

In a second test we used the exact data for (4.1) and k = 50 as input to our reconstruction method with p = 25, q = 32. As initial approximation we chose (4.1) with 0.2 replaced by 0.1. After 3 sweeps of our algorithm we obtained a reconstruction very similar to a low pass filtered version of f with cut off 50.

Finally we created a 3D breastphantom, see Fig. 1. It consists of fat, glandular tissue, a tumor and a cyst. The breast is suspended in a cube of sidelength 12 cm which is filled with water. f is given by

f = \frac{c_{0}^{2}}{c^{2}} − 1 − i\frac{2αc_{0}}{kc}  (4.2)

with c_{0} = 1500 m sec^{−1} the speed of sound in water. The values of c and α (at 1MHz) are

tissue  c[\frac{m}{sec}] α[\frac{db}{m}]  Ref  Imf
fat  1458 41  0.058 −9.4/k
glandular tissue 1519 80 −0.025 −18.4/k
tumor  1564 118 −0.080 −27.2/k
cyst  1568 10 −0.084 −2.3/k

Since the top face of the cube is not accessible, we have to modify the finite difference method (3.1). We stipulate the boundary condition ∂w/∂ν = 0 on the top face of Ω_{j}, i.e. we let n run from −q to q and put w_{ℓ,m,q+1} = w_{ℓ,m,q−1}. The boundary value problem (2.5) for z has to be changed accordingly. Of course this procedure is questionable, but at the present state of our work we just don’t know anything better.

We generated data for p = 32 equally spaced directions in [0, 2π] using our initial value method with q = 32, i.e. h = 6cm/32 = 1.875mm. The frequency of the iradiating waves was chosen to be 250 KHz, i.e. k = 10.47cm^{−1}. This corresponds to a wavelength of 6mm. After three sweeps of our algorithm we obtained a reconstruction in which the cyst and the tumor where clearly visible and distinguishable. The computing time per sweep on a SPARC 20 was 10 minutes.

References

[1] Borup, D.T. - Johnson, S.A. - Kim, W.W. - Berggren, M.J.: Nonperturbative diffraction tomography via Gauss-Newton iteration applied to the scattering integral equation, Ultrasonic Imaging 14, 69-85 (1992).

[2] Colton, D. - Monk, P.: A modified dual space method for solving the electromagnetic inverse scattering problem for an infinite cylinder, Inverse Problems 10, 87-108 (1994).

[3] Chew, W.C. - Wang, Y.M.: Reconstruction of Two-Dimensional Permittivity Distribution Using the Distorted Born Iteration Method, IEEE Transactions on Medical Imaging 9, 218-225 (1990).

[4] Devaney, A.J.: A filtered backpropagation algorithm for diffraction tomography, Ultrasonic Imaging 4, 336-350 (1982).

[5] Gutman, S. - Klibanov, M.: Three-Dimensional inhomogeneous media imaging, Inverse Problems 10, 239-246 (1994).

[6] Johnson, S.A.: Progress in applied and inverse scattering for ground penetrating radar, sonar, and medical ultrasound. Lecture on the IMA Workshop. Inverse Problems in Wave Propagation, March 6-17, 1995.

[7] Kleinman, R.E. - van den Berg, P.M.: A modified gradient method for two-dimensional problems in tomography, J. Comp. Appl. Math. 42, 17-35 (1992).

[8] Keller, J.B. - Papadakis, J.S.: Wave Propagation and Underwater Acoustics. Lecture Notes in Physics, vol. 70, Springer 1977.

[9] Knightly, G.H. - Mary, D.F.St.: Stable marching schemes based on elliptic models of wave propagation, J. Acoust. Soc. Am. 93, 1866-1872 (1993).

[10] Natterer, F.: Finite difference methods for inverse problems and Applications to Geophysics, Industry, Medicine and Technology, Proceedings of the International Workshop on Inverse Problems, January 17-19, 1995, HoChiMinh City. Publications of The HoChiMinh City Mathematical Society 2, 1995.

[11] Natterer, F.: An initial value approach to the inverse Helmholtz problem at fixed frequency. Proceedings of the Conference “Inverse Problems in Medical Imaging and Nondestructive Testing”, Oberwolfach, February 4-10, 1996. To appear in Springer 1996.

[12] Natterer, F.: Numerical solution of bilinear inverse Problems, Technical Report 19/96-N, Fachbereich Mathematik der Universität Münster.

[13] Natterer, F. - Wübbeling, F.: A finite difference method for the inverse scattering problem at fixed frequency, Lecture Notes in Physics, Vol. 422, p. 157-166, Springer 1993.

[14] Natterer, F. - Wübbeling, F.: A propagation - backpropagation method for ultrasound tomography, Inverse Problems 11, 1225-1232 (1995).

[15] Stenger, F. - O’Reilly, M.: A new approach to sonic inversion, to appear in Inverse Problems.

[16] Wübbeling, F.: Direktes und inverses Streuproblem bei fester Frequenz, Thesis, Fachbereich Mathematik, Universität Münster (1994).

0LaTeX: berg f berger paperfne Aix-les-Bains.tex

Fig.: Breast phantom at 200 KHz. Real part right, imaginary part left. Note that cyst and tumor can be distinguished clearly in the imaginary part.