Inverse Problems 13 (1997) 1341–1361. Printed in the UK PII: S0266-5611(97)83347-3
A fast and accurate imaging algorithm in optical/diffusion
tomography
M V Klibanov, T R Lucas and R M Frank
Department of Mathematics, University of North Carolina at Charlotte, Charlotte, NC 28223,
USA
Received 11 April 1997
Abstract. An n-dimensional (n D 2; 3) inverse problem for the parabolic/diffusion equation
ut D div .D.x/ru/ − a.x/u, u.x; 0/ D .x − x0/, x 2 Rn, t 2 .0; T / is considered. The
problem consists of determining the function a.x/ inside of a bounded domain   Rn given
the values of the solution u.x; t/ for a single source location x0 2 @ on a set of detectors
fxigmiD1  @, where @ is the boundary of . A novel numerical method is derived and
tested. Numerical tests are conducted for n D 2 and for ranges of parameters which are realistic
for applications to early breast cancer diagnosis and the search for mines in murky shallow water
using ultrafast laser pulses. The main innovation of this method lies in a new approach for a
novel linearized problem (LP). Such a LP is derived and reduced to a well-posed boundary-value
problem for a coupled system of elliptic partial differential equations. A principal advantage of
this technique is in its speed and accuracy, since it leads to the factorization of well conditioned,
sparse matrices with non-zero entries clustered in a narrow band near the diagonal. The authors
call this approach the elliptic systems method (ESM). The ESM can be extended to other imaging
modalities.
1. Introduction
In many applications of inverse problems any numerical method with a claim to be
practically valuable must be both fast and accurate. In this paper we present a new imaging
algorithm with the last two properties. At least two applications of this technique are to
early breast cancer diagnosis and the search for mines in murky water using ultrafast laser
pulses. In the recently published report of the National Research Council (NRC) [27] a
call was made to develop effective imaging algorithms for optical/diffusion tomography
(OT/DT), because of its applications to medical imaging. OT is a rapidly growing field.
We consider this publication as our response to this call of the NRC.
The core of the numerical method of this paper consists in a novel approach for a
new LP. Such an LP is derived and reduced to a well-posed boundary-value problem for
a coupled system of elliptic partial differential equations (PDEs). On the other hand, the
classical theory of numerical methods for PDEs implies that, unlike a conventional integral
equation, a well-posed boundary-value problem for a partial differential equation(s) can be
solved by the factorization of a well-conditioned , sparse matrix, with non-zero elements
clustered in a narrow band near the diagonal, which can be done rapidly. We call our
technique the elliptic systems method (ESM). The only pre-computation required is the
solution of the parabolic forward problem (2.1) for a background medium for a single
source location, which is done rapidly.
0266-5611/97/051341+21$19.50 c© 1997 IOP Publishing Ltd 1341

1342 M V Klibanov et al
Bukhgeim and Kazantsev [5] have demonstrated (though without numerical examples)
that the problem of computed tomography can be reduced to an elliptic system of the
Beltrami type. Note that the idea of [5] is different from ours due to the difference of the
problems.
In the conventional approach to numerical methods for OT and, more generally, to
numerical methods for multidimensional inverse problems, these problems are reduced to
ill-posed integral equations, which are generated by Green’s functions, cf Arridge [2],
Barbour et al [4], Cai et al [7], O’Leary et al [29], and references cited therein. This
approach leads to two time-consuming procedures.
(1) The solution of an ill-posed integral equation. This solution is especially time
consuming in the case of imaging inclusions whose sizes are small compared with the sizes
of the media of interest (tumours and mines, for example), because of the necessity of using
a fine mesh. This usually leads to the factorization of dense, ill-conditioned matrices with
millions of non-zero entries. Perhaps the most time-consuming procedure here is a choice
of the regularization parameter, which differs essentially for different matrix structures. The
matrix structure, in turn, depends critically on the background medium.
(2) Pre-computation and storing of the Green’s function for the background medium,
which is equivalent to the solution of many forward problems for many source locations.
We note that in medical applications, for example, the background media are heterogeneous
and vary essentially from one patient to another, cf Barbour et al [4].
Some recently developed numerical methods have suggested some strong approaches
for fast inverse solvers, see Colton and Kirsch [10], Elayaan and Isakov [13], Souriau et al
[31], and van den Berg and Kleinman [32]. However, these algorithms still rely either on
lengthly pre-computations of Green’s functions, or on explicit formulae in the case of a
homogeneous background medium. Recently Colak et al [11] suggested a semi-empirical
approach, which reduces OT to conventional computed tomography [8]. However, this
approach works only for a homogeneous background medium, whereas breast tissues are
essentially heterogeneous and are different for different patients, cf Barbour et al [4]. A fast
algorithm for ultrasound imaging has been developed by Natterer and Wu¨bbeling [28]. The
numerical method of [28] works for a heterogeneous background without pre-computation
of the Green’s function. However, this algorithm works only with high frequencies, which
usually do not propagate far into the media.
In section 2 we state the inverse problem under consideration and discuss its applications.
In section 3 the ESM is derived. In section 4 we prove the well-posedness of the resulting
elliptic boundary-value problem. In section 5 we discuss some details of the numerical
implementation. Section 6 is devoted to numerical experiments. In section 7 we discuss
the obtained results. In section 8 we outline extensions of the ESM to other data collection
schemes. Finally, theorems A2 and A3 in the appendix establish an asymptotic behaviour
of the function u.x; t/ as t ! 0, which is an important part of the ESM.
2. Inverse problem under consideration and its applications
Consider the Cauchy problem for the parabolic/diffusion equation
ut D div.D.x/ru/− a.x/u x 2 Rn; n D 2; 3I t 2 .0; T /
ujtD0 D .x − x0/: (2.1)
We will deal with a single source location fx0g. Let   Rn be a bounded domain with a
piecewise smooth boundary @. Let fxigmiD1  @ be a set of detectors placed around @.