SHAPE RECONSTRUCTION IN SCATTERING MEDIA WITH VOIDS USING A TRANSPORT MODEL AND LEVEL SETS

OLIVER DORN

Abstract. A two-step shape reconstruction method for diffuse optical tomography (DOT) is presented which uses adjoint fields and level sets. The propagation of near-infrared photons in tissue is modeled by the timedependent linear transport equation, of which the absorption parameter has to be reconstructed from boundary measurements. In the shape reconstruction approach, it is assumed that the inhomogeneous background absorption parameter and the values inside the obstacles (which typically have a high contrast to the background) are known, but that the number, sizes, shapes, and locations of these obstacles have to be reconstructed from the data. An additional difficulty arises due to the presence of socalled clear regions in the medium. The first step of the reconstruction scheme is a transport-backtransport (TBT) method which provides us with a low-contrast approximation to the sought objects. The second step uses this result as an initial guess for solving the shape reconstruction problem. A key point in this second step is the fusion of the ’level set technique’ for representing the shapes of the reconstructed obstacles, and an ’adjoint-field technique’ for solving the nonlinear inverse problem. Numerical experiments are presented which show that this novel method is able to recover one or more objects very fast and with good accuracy.

1. Introduction

1.1. A novel two-step shape reconstruction scheme. In a recent paper [12, 13] a two-step shape reconstruction method was introduced for the situation of 2D cross-borehole electromagnetic (EM) imaging. In that paper, electromagnetic crossborehole tomography was modeled as an inverse problem for the 2D Helmholtz equation. By reformulating the inverse problem as a shape reconstruction problem, and applying a novel level set based shape reconstruction scheme, it was possible to recover nontrivial shapes from noisy, limited-view, multifrequency data.

The method presented in [12, 13] is quite general. In the present paper, we demonstrate how this novel two-step shape reconstruction scheme can be applied to a completely different imaging situation. We are now interested in the imaging of absorption characteristics in human tissue from near-infrared light data. This imaging technique has become known as Diffuse Optical Tomography (DOT), see for example [1]. In contrast to the majority of approaches in DOT which use a diffusion approximation to the linear transport equation for the inversion task, we consider here the (time-dependent) linear transport equation itself, which is generally accepted as being the correct model for the photon propagation in tissue.

This article was written during the stay of the author at the Mathematical Sciences Research Institute (MSRI), Berkeley, from October to December 2001. The research presented here is supported in part by the NSERC CRD Grant 80357, and research during the stay at MSRI is supported by NSF grant DMS-9810361.

The inversion method is a two-step shape reconstruction scheme which tries to recover the unknown shapes of abnormalities in the absorption coefficient in the tissue. The first step consists of an inverse scattering algorithm for DOT which is called a ’Transport-Backtransport’ (TBT) method and which was originally introduced in [9, 10]. In the situation of DOT, this TBT algorithm yields already during the early iterations a relatively good approximation to the locations of the unknown obstacles, but requires many more iterations in order to get better estimates for the correct extension and shapes of these objects. Therefore, the result of these early sweeps is only used as an initial guess for initiating the second step of the two-step scheme, which is a level set based shape reconstruction method. This second step calculates successive corrections for the shapes of the obstacles which aim to reduce the misfit in the data. Here, the level set representation of the shapes is chosen because it enables us to easily describe and keep track of complicated geometries which might arise during the evolution process.

Both, the initialization step as well as the level set based shape evolution step employ an ’adjoint-field technique’ for the inversion which has the very useful property that the inverse problem can be solved approximately by making two uses of the same forward modeling code. Using a somewhat oversimplified description of our technique, the updates to the level set function are obtained by first making one pass through the code using the parameter distribution (which is here the absorption parameter) corresponding to the latest best guess of the level set function, and then another pass with the adjoint operator applied to the differences in computed and measured data. Then the results of these two calculations are combined to determine updates to the level set function. The resulting procedure is iterative, and can be applied successively to parts of the data, e.g., data associated with one source location can be used to update the model before other source locations are considered. Similar techniques have been successfully applied recently as part of iterative nonlinear inversion schemes in [6, 10, 11, 32, 44]. For alternative inversion schemes in DOT see for example [1, 17, 18, 20, 25, 27, 30, 45] and the references therein.

1.2. The level set method. The level set method was originally developed by Osher and Sethian for describing the motion of curves and surfaces [36, 41]. Since then, it has found applications in a variety of quite different situations. Examples are image enhancement, computer vision, interface problems, crystal growth, or etching and deposition in the microchip fabrication. For an overview see [42].

The idea of using a level set representation as part of a solution scheme for inverse problems involving obstacles was first suggested by Santosa in [40]. In that paper, two linear inverse problems, a deconvolution problem and the problem of reconstructing a diffraction screen, are solved by employing an optimization approach as well as a time evolution approach using the level set technique for describing the shapes. Santosa also outlines in that paper how the two presented methods can be generalized to nonlinear shape recovery problems. More recently, a related method was applied to a nonlinear shape recovery problem by Litman et al. in [28]. In that work, an inverse transmission problem in free space is solved by a controled evolution of a level set function. This evolution is governed by a Hamilton-Jacobi type equation, whose velocity function has to be determined properly in order to minimize a given cost functional.

In Dorn et al. [12, 13] a two-step shape reconstruction algorithm was presented for the situation of electromagnetic cross-borehole tomography in 2D. In contrary to the work mentioned above, this approach is not based on a Hamilton-Jacobi type equation, but uses an optimization approach instead. Moreover, a very specific inversion routine, an adjoint-field technique, is employed which calculates an update for the most recent level set function directly by just solving one forward and one adjoint problem on the computer. The present paper is based on this approach, and a detailed discussion will be given in the following sections.

As further examples for related and very interesting recent approaches we also want to mention here the work described in Burger [2], Elangovan et al [16], Hintermüller et al [21], Ito et al [22], Jackowska-Strumillo et al [23], Osher et al [37], Ramananjaona et al. [38], Sethian et al [43], and Zhao et al [46].

The paper is organized as follows. In section 2 we will present the basic equations of DOT in a form convenient for development of the shape reconstruction technique. In section 3 we formulate the shape reconstruction problem and introduce the level set formulation of this problem. In section 4, we derive the basic shape reconstruction algorithm using level sets and adjoint fields. Section 5 describes shortly the procedure employed here for finding a good starting guess for our shape recontruction algorithm. In Section 6 numerical experiments are presented which demonstrate the performance of the algorithm in different situations. The final section summarizes the results of this paper and indicates some directions for future research.

2. The physical experiment in DOT

2.1. The linear transport equation. We model the propagation of photons in the medium by the time-dependent linear transport equation

\frac{∂u}{∂t} + θ · ∇u(x, θ, t) + (a(x) + b(x))u(x, θ, t) − b(x) ∫_{S^{n−1}} η(θ · θ′)u(x, θ′, t)dθ′

= q(x, θ, t)  in Ω × S^{n−1} × [0, T]  (1)

with initial condition

u(x, θ, 0) = 0  in Ω × S^{n−1}  (2)

and boundary condition

u(x, θ, t) = 0  on Γ_{−}.  (3)

Here,

Γ_{±} := {(x, θ, t) ∈ ∂Ω × S^{n−1} × [0, T], ±ν(x) · θ > 0} .

Ω is a convex, compact domain in IR^{n}, n = 2, 3, with smooth boundary ∂Ω. In our numerical experiments, we will only consider the case n = 2, but the algorithm extends in a straightforward way to n = 3. ν(x) denotes the outward unit normal to ∂Ω at the point x ∈ ∂Ω, and u(x, θ, t) describes the density of particles (photons) which travel in Ω at time t through the point x in the direction θ. The velocity c of the particles is assumed to be normalized to c = 1cms^{−1} and has been neglected in the formulation of (1).

a(x) is the absorption cross-section (in short ’absorption’), b(x) is the scattering cross-section, and μ(x) := a(x)+b(x) is the total cross-section or attenuation. These parameters are assumed to be real, strictly positive functions of the position x. The quantity μ^{−1} is the mean free path of the photons. Typical values in DOT are a ≈ 0.1 − 1.0cm^{−1}, b ≈ 100.0 − 200.0cm^{−1}, μ^{−1} ≈ 0.005 − 0.01cm, see for example [24]. The scattering function η(θ ·θ′) describes the probability for a particle entering a scattering process with the direction of propagation θ′ to leave this process with the direction θ. It is normalized to

∫_{S^{n−1}} η(θ · θ′)dθ = 1,  (4)

which expresses particle conservation in pure scattering events. The dot-product in the argument indicates that η depends only on the cosine of the scattering angle cos ϑ = θ · θ′ and, in particular, is independent of the location of the scattering event x. These latter assumptions simplify the following calculations, but the reconstruction method can be extended to more general scattering functions as well. In our numerical experiments, we will use the following Henyey-Greenstein scattering function [24]

η(θ · θ′) = \frac{1 − g^{2}}{2(1 + g^{2} − 2g cos ϑ)^{3/2}} ,  (5)

with −1 < g < 1. The parameter g in (5) is the mean cosine of the scattering function. Values of g close to one indicate that the scattering is primarily forward directed, whereas values close to zero indicate that scattering is almost isotropic. In our numerical experiments, we will choose g to be 0.9, which is a typical value for DOT.

The initial condition (2) indicates that there are no photons moving inside of Ω at the starting time of our experiment. The boundary condition (3) indicates that during the experiment no photons enter the domain Ω from the outside. All photons inside of Ω originate from the source q which, however, can be situated at the boundary ∂Ω. We mention that an alternative way of describing laser sources at the boundary is to use an inhomogeneous boundary condition instead of (3) and a zero source q = 0. Both choices are equivalent, and the derivation of the TBTmethod as well as the shape reconstruction method are almost identical for these two possibilities.

We consider the problem (1)-(3) for p different sources q_{j}, j = 1, …, p, positioned at ∂Ω. Typical sources in applications are delta-like pulses transmitted at time t = 0 at the position s_{j} ∈ ∂Ω into the direction θ_{j}, which can be described by the distributional expressions

q̃_{j}(x, θ, t) = δ_{x}(s_{j})δ_{t}(0)δ_{θ}(θ_{j}),  j = 1, …, p,  (6)

with ν(s_{j}) · θ_{j} < 0. We assume that a given source q_{j} gives rise to the physical fields ũ_{j}(x, θ, t) which are solutions of (1)-(3). Our measurements consist of the outgoing flux across the boundary ∂Ω which has with (3) the form

G̃_{j}(x, t) = ∫_{S^{n−1}} ν(x) · θũ_{j}(x, θ, t)dθ on ∂Ω × [0, T],  (7)

for j = 1, …, p. With these definitions, we are now able to formulate the inverse problem which we will consider in the following.

Definition 2.1. Inverse problem of DOT: Assume that we measure for p different sources q_{j}, j = 1, …, p, the data (7). Provided with this information and knowing the scattering functions η and b, how can we reconstruct the coefficient a(x) of (1) inside of Ω?

In general, also the scattering parameter b(x) is part of the inversion problem. However, for simplicity, we will restrict ourselves in this paper to the determination of only one parameter, namely a(x), and will assume that the scattering parameter b(x) is given. Extensions of this method will be addressed in our future research.

2.2. Clear regions. An additional difficulty arises in our application of DOT due to the presence of so-called ’clear regions’ in the imaging domain. These are regions where the scattering coefficient b(x) is very small or zero, such that the dominant behavior of the traveling photons inside these regions is transport rather than diffusion. Most of the imaging methods in DOT are purely based on the diffusion approximation to the linear transport equation, such that they cannot be used in applications where clear regions are present. Clear regions play an important role for example when imaging the head of neonates or in functional brain imaging. See [1, 26, 35, 39].

Although we assume in this paper that the locations of the clear regions are known a priori (which in realistic applications is not always the case such that the shapes of the clear regions become part of the imaging task themselves), the presence of these regions complicates the numerical inversion process.

One problem caused by the clear layers is their effect on the sensitivity structure of the data with respect to the unknown parameters. The clear regions (and, in particular, the clear layers which are surrounding the region of interest, see figure 1) act like ’wave guides’ for the photons, since the propagation is much faster inside these layers than it is in the diffusive regions. Therefore, the photons which have travelled through these clear layers dominate in particular the early parts of the time-dependent measurements. Unfortunately, these photons do not carry much information related to the internal structure of the medium. In absence of clear layers, the early photons are usually assumed to yield the reconstructions with the best resolution. If clear layers are present, however, these early photons are now mixed with those photons that travelled through the clear layers, such that it is difficult to extract the necessary information about the interior domain from the early data. Therefore, we will in our numerical reconstructions only use photons which arrived at later times at the receiver locations. Our sensitivity analysis indicates that these later photons are less affected by the clear layers. We will show that the information of the later arriving photons is sufficient for achieving good reconstructions of the absorption distribution inside the regions of interest, which might however have a slightly reduced resolution. See also the more detailed discussion of sensitivity functions in section 4.7.

3. The shape reconstruction problem

In this section we formulate the shape reconstruction problem which we want to solve, and cast it in a form which makes use of the level set representation of the domains.

3.1. The shape reconstruction problem in DOT. To start with we introduce some terminology which we will use throughout the paper.

Definition 3.1. Let us assume that we are given a constant â > 0, and a bounded function a_{b} : Ω → IR. We call a pair (D, a), which consists of a compact domain D ⊂⊂ Ω and a bounded function a : Ω → IR, admissible if we have

a|_{D} = â, a|_{Ω\D} = a_{b}|_{Ω\D}.  (8)

In other words, a pair (D, a) is admissible if a is equal to a preassigned constant value â inside of D, and equal to the preassigned background absorption a_{b} outside of D. The domain D is called the scattering domain.

Remark 3.1. For an admissible pair (D, a), and for given â, a_{b}, the absorption a is uniquely determined by D.

With this definition, we can now formulate the shape reconstruction problem.

Shape reconstruction problem. Let us assume that we are given a constant â > 0, a bounded function a_{b} : Ω → IR, and some data G̃_{j} as in (7). Find a domain D̃ such that the admissible pair (D̃, ã) reproduces the data, i.e. (7) holds with ũ_{j} being a solution of (1)-(3) with a = ã for j = 1, …, p.

Using the same notation and assumptions as in definition 3.1, we want to formulate another inverse problem which we will call the inverse scattering problem and which will play an important part when solving the shape reconstruction problem.

Inverse Scattering Problem. Let us assume that we are given a bounded function a_{b} : Ω → IR, and some data G̃_{j} as in (7). Find a bounded function ã_{s} : Ω → IR such that ã = a_{b} + ã_{s} reproduces the data, i.e. (6) holds with ũ_{j} given by (1)-(3) and with a = ã for j = 1, …, p.

The inverse scattering problem gives rise to the decomposition

a = a_{b} + a_{s} in Ω  (9)

In other words, the absorption distribution a is decomposed into the background distribution a_{b} and the perturbation a_{s} which we will refer to in the following as the obstacle.

Solving the shape reconstruction problem requires only to find the shape of the domain D̃, since the function ã is then uniquely determined by (8). Solving the inverse scattering problem, on the other hand, amounts to finding the entire function ã_{s} from the given data, which is much harder to do. However, it will turn out that finding a good approximate solution of the inverse scattering problem is much easier and faster to achieve and will provide us with a very good initial guess for starting our shape reconstruction routine.

Definition 3.1 allows us to formulate a first version of the strategy which we want to use for solving the shape reconstruction problem.

Strategy for solving the shape reconstruction problem. Construct a series of admissible pairs (D^{(n)}, a^{(n)}), n = 0, 1, 2, …, such that the misfit between the physically measured data and the calculated data corresponding to (D^{(n)}, a^{(n)}) decreases with increasing n, and ideally, i.e. in absence of noise, tends to zero in the limit n → ∞. Use the approximate solution of the inverse scattering problem to initialize this series by determining a good starting element (D^{(0)}, a^{(0)}).

3.2. The domains D^{(n)}. In our numerical examples, each of the domains D^{(n)} which we are looking for can be given as a collection of a finite number L_{n} of disjoint, compact subdomains D_{l}^{(n)}, l = 1, …, L_{n}, with

D^{(n)} = ⋃_{l=1}^{L_{n}} D_{l}^{(n)},  D_{l}^{(n)} ∩ D_{l′}^{(n)} = ∅ for l ≠ l′.  (10)

The shapes of these subdomains D_{l}^{(n)} can in principle be arbitrary. In particular, they are allowed to be multiply connected, and to enclose some ’cavities’ or ’holes’ filled with background material, and moreover the number L_{n} of these subdomains might (and usually does) vary with the iteration number n. For the derivation of the inversion method, we assume that the boundaries ∂D_{l}^{(n)} of these domains are sufficiently smooth.

It is essential for the success and the efficiency of the reconstruction scheme to have a good and flexible way of keeping track of the shape evolution during the reconstruction process. The method we have chosen in our reconstruction algorithm is a level set representation of the shapes as it was suggested by Santosa [40]. This representation has the advantage that the level set functions, which are in principle only used for representing the shapes, can in a natural way be made part of the reconstruction scheme itself. Doing so, it is not necessary anymore to refer to the shapes of the domains until the reconstruction process is completed. The final shape is then recovered from the representing level set function easily. In the following we will discuss in a more formal way how this can be achieved.

3.3. Level set representation of the domains D^{(n)}. Assume that we are given a domain D ⊂⊂ Ω. The characteristic function χ_{D} : Ω → {0, 1} is defined in the usual way as

χ_{D}(x) = \{ \begin{matrix}1 , & x ∈ D \\ 0 , & x ∈ Ω\D. \end{matrix}  (11)

Definition 3.2. We call a function ϕ : Ω → IR a level set representation of D if

χ_{D}(x) = Ψ_{ϕ}(x)  on Ω  (12)

where Ψ_{ϕ} : Ω → {0, 1} is defined as

Ψ_{ϕ}(x) = \{ \begin{matrix}1 , & ϕ(x) ≤ 0 \\ 0 , & ϕ(x) > 0. \end{matrix}  (13)

For each function ϕ : Ω → IR there is a domain D associated with ϕ by (12),(13) which we call D[ϕ]. It is clear that different functions ϕ_{1}, ϕ_{2}, ϕ_{1} ≠ ϕ_{2}, can be associated with the same domain D[ϕ_{1}] = D[ϕ_{2}], but that different domains cannot have the same level set representation. Therefore, we can use the level set representation for unambiguously specifying a domain D by any one of its associated level set functions.

The boundary Γ = ∂D[ϕ] of a domain D[ϕ], represented by the level set function ϕ, is defined as

Γ = {x ∈ Ω : for all ρ > 0 we can find x_{1}, x_{2} ∈ B_{ρ}(x)  (14)
with ϕ(x_{1}) > 0 and ϕ(x_{2}) < 0 }

Definition 3.3. We call a triple (D, a, ϕ), which consists of a domain D ⊂⊂ Ω and bounded functions a, ϕ : Ω → IR, admissible if the pair (D, a) is admissible in the sense of definition 3.1, and ϕ is a valid level set representation of D.

Remark 3.2. For an admissible triple (D, a, ϕ), and for given â, a_{b}, the pair (D, a) is uniquely determined by ϕ.

We use these definitions to reformulate our shape reconstruction problem.

Level set formulation of the shape reconstruction problem. Given a constant â > 0, a background distribution a_{b}, and some data G̃ as in (7). Find a level set function ϕ̃ such that the corresponding admissible triple (D̃, ã, ϕ̃) reproduces the data, i.e. (7) holds with ũ_{j} being a solution of (1)-(3) with a = ã for j = 1, …, p.

The strategy for solving this shape reconstruction problem has to be reformulated, too. It reads now as follows.

Strategy for solving the reformulated shape reconstruction problem. Construct a series of admissible triples (D^{(n)}, a^{(n)}, ϕ^{(n)}), n = 0, 1, 2, …, such that the misfit between the data (7) and the calculated data corresponding to (D^{(n)}, a^{(n)}, ϕ^{(n)}) decreases with increasing n, and ideally, i.e. in absence of noise, tends to zero in the limit n → ∞. For finding this series we only have to keep track of a^{(n)} and ϕ^{(n)}, but not of D^{(n)}. The function a^{(n)} is needed in each step for solving a forward problem (1)-(3), and a corresponding adjoint problem. The knowledge of ϕ^{(n)} is used in each step to determine a^{(n)}. The final level set function ϕ^{(N)}, which satisfies some stopping criterion, is used to recover the final shape D^{(N)} via (12).

4. Solving the shape reconstruction problem

In this section we derive the basic shape reconstruction method which uses adjoint fields and the level set representation introduced above. The initializing procedure for this reconstruction routine will be discussed in section 5.

4.1. Function spaces. We want to specify now the function spaces which we will be working with. The main objective of this section is to introduce the inner products on these function spaces, which will become important when defining the adjoint linearized operators in the following sections.

The natural spaces in the framework of the linear transport equation are L^{1}- spaces for the fields u, and L^{∞}-spaces for the parameters a and b, due to the physical significance of the corresponding norms. This is because the L^{1}-norm of a function u(x, θ, t_{0}) describes physically the number of particles which are in the system at a given time t_{0}. However, in our derivation of the shape reconstruction method it will be more convenient to work with L^{2}-spaces instead, since deriving the adjoint operators used in the inversion routine is more straightforward when using these spaces. See for example [5, 7, 9, 10] for more information concerning L^{1}-spaces in linear transport theory. The spaces which we will consider are

Y = L^{2}(Ω × S^{n−1} × [0, T]), U = L^{2}(Ω × S^{n−1} × [0, T]),  (15)
Z = L^{2}(∂Ω × [0, T]), F = L^{2}(Ω), Φ = L^{2}(Ω),  (16)

equipped with the standard inner products. Here, Y is the space of sources q, U is the space of states u, Z is the space of measurements or data G, F is the space of parameters a, and Φ is the space of level set functions ϕ. We mention that in fact our level set functions ϕ will be assumed to have a certain regularity, since in the derivation of the reconstruction scheme we will make use of first order derivatives of these functions.

4.2. Operators. In the following, we will introduce some operators which will enable us to formulate the shape reconstruction problem in a way suitable for deriving the inversion algorithm.

Given a constant â and a bounded function a_{b} : Ω → IR. Then, with each level set function ϕ ∈ Φ a uniquely determined ’absorption perturbation’ or ’obstacle’ Λ(ϕ) is associated by putting

Λ(ϕ)(x) = \{ \begin{matrix}â − a_{b}(x) , & ϕ(x) ≤ 0 \\ 0  , & ϕ(x) > 0. \end{matrix}  (17)

With (13) we can write this also as

Λ(ϕ)(x) = Ψ_{ϕ}(x)(â − a_{b}(x)) , x ∈ Ω.  (18)

Notice that the operator Λ is chosen such that the triple (D, a, ϕ) with a = a_{b} +Λ(ϕ) and domain D[ϕ] forms an admissible triple (D, a, ϕ) in the sense of definition 3.3. Moreover, for (D, a, ϕ) an admissible triple, we see that Λ(ϕ) is just the absorption perturbation a_{s}

Λ(ϕ)(x) = a_{s}(x) = χ_{D}(x)(â − a_{b}(x)),  x ∈ Ω.  (19)

Let us assume now that we are given a background absorption a_{b} and that we have collected some data G̃_{j} which correspond to the ’true’ absorption distribution

ã = a_{b} + ã_{s},  (20)

where ã_{s} is the ’true’ obstacle. We consider the following measurement operators

G_{j} : F ⟶ Z, G_{j}(a_{s})(x, t) := ∫_{S^{n−1}} ν(x) · θ u_{j}(x, θ, t) dθ,  (21)

j = 1, …, p, where u_{j} solves (1)-(3) with the source q_{j} ∈ Y and the parameters a = a_{b} + a_{s}. For the ’true’ parameters ã_{s} we ask (21) to coincide with the data G̃_{j}

G_{j}(ã_{s})(x, t) = G̃_{j}(x, t) for j = 1, …, p  (22)

on ∂Ω × [0, T]. Assuming G̃_{j} ∈ Z we want to determine ã_{s} such that (22) is valid. For given data G̃_{j}, j = 1, …, p, we define furthermore the residual operators

R_{j} : F ⟶ Z, R_{j}(a_{s})(x, t) = G_{j}(a_{s})(x, t) − G̃_{j}(x, t).  (23)

From (22) it follows that for the ’true’ obstacle the residuals vanish,

R_{j}(ã_{s}) = 0  for j = 1, …, p,  (24)

if the data are noise-free.

The forward operators T_{j} which map a given level set function ϕ ∈ Φ into the corresponding mismatch in the data are defined by

T_{j} : Φ ⟶ Z_{j} , T_{j}(ϕ) = R_{j}(Λ(ϕ))  (25)

for j = 1, …, p. The goal is to find a level set function ϕ̃ ∈ Φ such that

T_{j}(ϕ̃) = 0  for j = 1, …, p.  (26)

We mention that all three operators Λ, R_{j} and T_{j} are nonlinear.

4.3. Linearized operators. For the derivation of the shape reconstruction algorithm, we will need expressions for the linearized operators corresponding to the nonlinear operators introduced above, and for their adjoints with respect to the given inner products. In this section, we define the linearized operators, and expressions for their adjoints are derived in section 4.5.

Analogous to Santosa [40] it is shown that, for a homogeneous background a_{b}, the infinitesimal response δa(x) in the obstacle a_{s}(x) to an infinitesimal change δϕ(x) of the level set function ϕ(x) has the form

δa(x) = − [â − a_{b}] \frac{δϕ(x)}{|∇ϕ(x)|}|_{x∈∂D[ϕ]} .  (27)

The function δa in (27) can be interpreted as a ’surface measure’ on the boundary Γ = ∂D[ϕ]. Similar to (27), we define the linearized operator Λ̃′[ϕ] by

(Λ̃′[ϕ]δϕ) (x) = − [â − a_{b}(x)] \frac{δϕ(x)}{|∇ϕ(x)|}δ̂_{Γ}(x)  (28)

where δ̂_{Γ}(x) denotes the Dirac delta distribution concentrated on Γ = ∂Ω[ϕ]. In this interpretation, (27) describes an infinitesimal ’surface load’ of absorption on Γ which has to be recovered from the mismatch in the data.

However, the expression on the right hand side of (28) is not an element of F which causes problems when we want to calculate the inner products defined in section 4.1. Therefore, we will replace the operator (28) by an approximation which maps from Φ into F and which will be more convenient for the derivation of the reconstruction method.

For a given level set function ϕ ∈ Φ, let Γ = ∂D[ϕ] and B_{ρ}(Γ) = ∪_{y∈Γ}B_{ρ}(y) a small finite width neighborhood of Γ with some given constant 0 < ρ ≪ 1. The (approximated) linearized operator Λ′[ϕ] is defined as

Λ′[ϕ] : Φ ⟶ F, (Λ′[ϕ]δϕ) (x) = − [â − a_{b}(x)] \frac{δϕ(x)}{|∇ϕ(x)|}C_{ρ}(Γ)χ_{B_{ρ}(Γ)}(x) (29)

with χ_{B_{ρ}(Γ)}(x) given by (11). Here, C_{ρ}(Γ) = L(Γ)/Vol(B_{ρ}(Γ)) where L(Γ) = ∫_{Ω} δ̂_{Γ}(x)dx is the length of the boundary Γ, and Vol(B_{ρ}(Γ)) = ∫_{Ω} χ_{B_{ρ}(Γ)}(x)dx is the volume of B_{ρ}(Γ). For a very small ρ we will get a very large weight C_{ρ}(Γ), whereas for increasing ρ this weight C_{ρ}(Γ) decreases accordingly. The operator defined in (29) maps now from Φ into F such that we can make use of the inner products defined on these spaces.

We already mentioned that the term |∇ϕ(x)| in (27), (28), as well as the derivation of these expressions, implies some regularity constraint on ϕ. For example, ϕ ∈ C^{1} would be possible. Another possibility would be to use a ’signed distance function’ as a standard representation of the boundary [42]. We do not want to specify the regularity of ϕ at this point, but assume instead that it is ’sufficiently smooth’ for our purposes.

The linearized residual operator R′_{j}[a_{s}] is defined by

Rj^{′} [a_{s}] : F ⟶ Z,  Rj^{′} [a_{s}] δa = ∫_{S^{n−1}} ν · θ v_{j}(x, θ, t) dθ  (30)

where v_{j} solves the problem

\frac{∂v_{j}}{∂t} + θ · ∇v_{j}(x, θ, t) + (a(x) + b(x)) v_{j}(x, θ, t) − b(x) ∫_{S^{n−1}} η(θ · θ′)v_{j}(x, θ′, t)dθ′

= −δa(x) u_{j}(x, θ, t) in Ω × S^{n−1} × [0, T]  (31)

with initial condition

v_{j}(x, θ, 0) = 0  in Ω × S^{n−1}  (32)

and boundary condition

v_{j}(x, θ, t) = 0  on Γ_{−}.  (33)

Here, u_{j} is the solution of (1)–(3) with source q_{j} and with a = a_{b} + a_{s}. Note that (31) has the same structure as (1). The ‘primary sources’ q of (1) are now replaced by ‘secondary sources’ which are caused by the parameter perturbation δa. The representation (30) can be derived by perturbing

a_{s} → a_{s} + δa , u_{j} → u_{j} + v_{j},  (34)

plugging this into (1) and neglecting terms which are of higher than linear order in the perturbations δa, v_{j}. A rigorous derivation of (30) in the framework of L^{1}-spaces can be found in [7].

As our third linearized operator, we introduce the linearized forward operator T′_{j}[ϕ] by putting

T′_{j}[ϕ] : Φ ⟶ Z_{j} , T′_{j}[ϕ]δϕ = R′_{j} [Λ(ϕ)] Λ′[ϕ]δϕ.  (35)

All three operators Λ′[ϕ], R′_{j}[a_{s}], and T′_{j}[ϕ] are linear.

4.4. A nonlinear Kaczmarz-type approach. The algorithm works in a ’singlestep fashion’ as follows. Instead of using the data (22) for all sources simultaneously, we only use the data for one source at a time while updating the linearized residual operator after each determination of the corresponding incremental correction δϕ.

To be more specific, let us assume that we are given a level set function ϕ^{(n)}(x) and an obstacle a(n)_{s}(x) such that (D^{(n)}, a_{b} + a(n)_{s}, ϕ^{(n)}) forms an admissible triple in the sense of definition 3.3. Using a data set G̃_{j} corresponding to the fixed source position q_{j}, we want to find an update δϕ^{(n)} to ϕ^{(n)} such that for the admissible triple

(D^{(n+1)}, a_{b} + a(n+1)_{s}, ϕ^{(n+1)}) :=  (36)
(D[ϕ^{(n)} + δϕ^{(n)}], a_{b} + Λ(ϕ^{(n)} + δϕ^{(n)}), ϕ^{(n)} + δϕ^{(n)})

the residuals in the data corresponding to this source vanish

T_{j}(ϕ^{(n+1)}) = T_{j}(ϕ^{(n)} + δϕ^{(n)}) = 0.  (37)

Applying a Newton-type approach, we get from (37) a correction δϕ^{(n)} for ϕ^{(n)} by solving

T′_{j}[ϕ^{(n)}]δϕ^{(n)} = − T_{j}(ϕ^{(n)}) = − (M_{j}u_{j} − G̃_{j})  (38)

where u_{j} satisfies (1)–(3) with

a(x) = a_{b}(x) + Λ(ϕ^{(n)})(x) = \{ \begin{matrix}â , & x ∈ D[ϕ^{(n)}] \\ a_{b}(x) , & x ∈ Ω \ D[ϕ^{(n)}]. \end{matrix}  (39)

Since we have only few data given for one source and one frequency, equation (38) usually will have many solutions (in the absence of noise), such that we have to pick one according to some criterion. We choose to take that solution which minimizes the L^{2}-norm of δϕ^{(n)}

Min ∥δϕ^{(n)}∥_{2} subject to T′_{j}(ϕ^{(n)})δϕ^{(n)} = − (M_{j}u_{j} − G̃_{j}) .  (40)

This solution can be formulated explicitly. It is

δϕ(n)_{MN} = − T′_{j}[ϕ^{(n)}]^{∗} (T′_{j}[ϕ^{(n)}]T′_{j}[ϕ^{(n)}]^{∗})^{−1} (M_{j}u_{j} − G̃_{j}) ,  (41)

where T′_{j}[ϕ^{(n)}]^{∗} denotes the adjoint operator to T′_{j}[ϕ^{(n)}].

After correcting ϕ by ϕ → ϕ+δϕ_{j}, where δϕ_{j} is given by (41), we use the updated residual equation (38) to compute the next correction δϕ_{j′}. Doing this for one equation after the other, until each of the sources q_{j} has been considered exactly once, will yield one complete sweep of the algorithm. This procedure is similar to the Kaczmarz method for solving linear systems, or the algebraic reconstruction technique (ART) in x-ray tomography [31] and the simultaneous iterative reconstruction technique (SIRT) as presented in [8]. For related approaches in a variety of imaging problems we refer to [6, 10, 11, 32, 33, 34, 44].

4.5. The adjoint linearized operators. In order to calculate the minimal norm solution (41), we will need practically useful expressions for the adjoints of the linearized operators of section (4.3). We will present such expressions in this section. The calculation of the actions of these operators will typically require us to solve an adjoint transport problem. This explains the name ’adjoint-field method’ of the inversion method employed here.

To start with, a simple calculation gives us the following theorem.

Theorem 4.1. The adjoint operator Λ′[ϕ]^{∗} which corresponds to the linearized operator Λ′[ϕ] is given by

Λ′[ϕ]^{∗} : F ⟶ Φ , (Λ′[ϕ]^{∗}δa_{s}) (x) = − [â − a_{b}(x)] \frac{δa_{s}(x)}{|∇ϕ(x)|}C_{ρ}(Γ)χ_{B_{ρ}(Γ)}(x).

(42)

The next theorem describes the adjoint operator R′_{j}[a_{s}]^{∗} which corresponds to R′_{j}[a_{s}]. Its proof can be found in [9, 10], and is therefore omitted here.

Theorem 4.2. Let z ∈ U^{∗} be a solution of the (adjoint) transport equation

−\frac{∂z}{∂t} − θ · ∇z(x, θ, t) + (a(x) + b(x)) z(x, θ, t) − b(x) ∫_{S^{n−1}} η(θ · θ′) z(x, θ′, t)dθ′

= 0  in Ω × S^{n−1} × [0, T]  (43)

with ‘initial’ condition

z(x, θ, T) = 0  on Ω × S^{n−1}  (44)

and boundary condition

z|_{Γ_{+}} = z̃(x, t) := ζ(x, t).  (45)

Then the adjoint operator Rj^{′}[a_{s}]^{∗} is given by

(Rj^{′}[a_{s}]^{∗}ζ)(x) = I_{j}(a_{s}) := ∫_{[0,T]} ∫_{S^{n−1}} u_{j}(x, θ, t)z_{j}(x, θ, t) dθdt,  (46)

with u_{j} being a solution of (1)–(3) with source q_{j} and with a = a_{b} + a_{s}. Notice that ζ(x, t) is applied uniformly into all directions θ with ν · θ > 0. □

Finally, by combining theorems 4.1 and 4.2, we get an expression for the adjoint operator T′_{j}[ϕ]^{∗} which corresponds to the linearized forward operator T′_{j}[ϕ]. It is described in the following theorem.

Theorem 4.3. Let ζ ∈ Z. The adjoint operator T′_{j}[ϕ]^{∗} acts on ζ in the following way

T′_{j}[ϕ]^{∗}ζ = Λ′[ϕ]^{∗}R′_{j}[Λ(ϕ)]^{∗}ζ = \frac{[â − a_{b}(x)]}{|∇ϕ(x)|}C_{ρ}(Γ)χ_{B_{ρ}(Γ)}(x)I_{j}(Λ(ϕ)) (47)

where I_{j}(Λ(ϕ)) is given by (46) and where u_{j} solves (1)-(3) and z_{j} solves (43)–(45) with a_{s} replaced by Λ(ϕ).

4.6. Updating the level set function. Let us consider the operator C_{j,n} := T′_{j}[ϕ^{(n)}]T′_{j}[ϕ^{(n)}]^{∗} in (41) more closely. In order to calculate a correction δϕ(n)_{MN} by (41) we first have to apply the operator C_{j,n}^{−1} to the expression M_{j}u_{j} − G̃_{j}. This operator is difficult to calculate, since typically a large number of forward problems has to be solved for this purpose if the background medium is inhomogeneous and the corresponding Green functions are not easily available. Moreover, the inversion of this operator is usually unstable due to the ill-posedness of the underlying inverse problem. Therefore, in our numerical experiments, we will usually just replace this operator by a suitably chosen approximate operator Ĉ_{j} which is easier to calculate and whose inversion is more stable. One possibility for such an approximation (which is adopted here) is to simply use a multiple of the identity, combined with an appropriate selection of receivers in the current step of the nonlinear Kaczmarz-type approach. We refer here for a more thorough discussion of this operator and its approximations to [9, 10, 11, 12, 13, 34]. With this approximation, we define

ζ_{j} := Ĉ_{j}^{−1} (M_{j}u_{j} − G̃_{j}) .  (48)

Next, in order to find δϕ(n)_{MN}, we have to calculate T′_{j}[ϕ^{(n)}]^{∗}ζ_{j}. Using (47) we get

T′_{j}[ϕ^{(n)}]^{∗}ζ_{j} = \frac{(â − a_{b}(x))}{|∇ϕ^{(n)}(x)|} C_{ρ}(Γ)χ_{B_{ρ}(Γ)}(x)I_{j}(Λ(ϕ^{(n)})),  (49)

where u_{j} solves (1)-(3) with source q_{j} and z_{j} solves (43)–(45) with ζ = ζ_{j} and a_{s} = Λ(ϕ^{(n)}).

In the following we want to introduce an approximation to the update formula (49) which in our numerical experiments so far has shown to simplify and stabilize the reconstruction process. We mentioned already earlier that the level set representation of a given contour is not uniquely defined. Therefore, it would be desirable to pick from all those possible level set functions representing a given contour that one with the most convenient properties.

Since we are updating the level set function in each step of our algorithm, it might happen that the denominator in (49) becomes very small, and numerical noise and roundoff-errors can destabilize the inversion routine. In order to avoid this problem, we rescale the level set function ϕ after each update with the aim to keep |∇ϕ(x)| near the boundaries of the inclusions in a certain range of values. As a simple scaling rule, we multiply the level set function with a scaling factor such that the global minimum (or alternatively maximum) of the level set function is always equal to a predefined constant value. Notice that this rescaling procedure does not change the domains which are represented by the given level set function. In other words, we try to pick from the family of level set functions representing a given domain that one which has a convenient behavior of |∇ϕ(x)| in the neighborhood of the zero level set.

We go one step further and make the assumption that, with this rescaling procedure, the gradient |∇ϕ(x)| along the boundary can be approximated in (49) by some constant c_{1}. Although, with this approximation, we are not marching anymore into the ’optimal’ direction when calculating the updates, we will gain some stability in the inversion routine since we avoid dividing by small numerical derivatives (these are still possible even after the rescaling) in (49) which are strongly influenced by unavoidable numerical noise and roundoff errors. Moreover, our numerical experiments show that these small deviations from the ’optimal’ marching direction are easily corrected by the succeeding updates.

With this modification, we get the following update formula for the level set function

δϕ̂^{(n)}(x) = − (â − a_{b}(x))c_{1}^{−1} C_{ρ}(Γ)χ_{B_{ρ}(Γ)}(x)I_{j}(Λ(ϕ^{(n)})).  (50)

In (50) we did not incorporate the rescaling procedure described above.

Notice that, although we did not explicitly impose any regularity constraints on the updates (50), they are now in the range of R′_{j}[Λ(ϕ^{(n)})]^{∗} (up to the factor â − a_{b}(x)) which implicitly gives us some information about the regularity we can expect. Our numerical experiments so far indicate that the degree of regularity which is achieved by applying (50) is typically sufficient ’for practical purposes’ in our numerical examples.

4.7. Sensitivity analysis. A valuable tool for designing an imaging experiment is a careful sensitivity analysis of the underlying inverse problem. The sensitivity functions which correspond to a nonlinear inverse problem are functions which decompose the linearized residual operator (which is often called ’Fréchet derivative’) as well as its adjoint operator, as described in more details in [14, 15]. The sensitivity functions corresponding to one given source position and one given receiver location describe the sensitivity of this source-receiver pair to small variations of the parameters in the medium. Locations in the medium where the sensitivities are small have almost no influence on the measurements at this specific detector location. Therefore, hidden objects at these locations are difficult to recover unless other source-receiver pairs have larger sensitivity values at these specific locations. Ideally, we would like to design an experiment where certain source-receiver configurations have a sensitivity structure which is large on a very small area in the domain, and very small elsewhere. The corresponding data can then be used for determining the coefficients at these specific locations. If we find a set of sensitivity functions such that each location in the medium is ’probed’ by such an ideal configuration, we expect to get good and fast reconstructions. Unfortunately, in diffusion-type inverse problems like ours, the sensitivity functions are typically spread out over a relatively large area. In addition, they have a large dynamical range since they have very large values close to the source and receiver locations, and very small values in the regions of interest. Therefore, it is very difficult to probe localized areas inside the domain with suitable superpositions of sensitivity functions. The role of the sensitivity functions in our TBT-reconstruction scheme is described in the following theorem.

Theorem 4.4. Transport sensitivity functions. Given a source (6). We can find functions Ψ_{a}(x_{r}, t_{r}; x_{sc}) with the following property. 1.)’Projection step’: The action of the linearized residual operator R′[a_{s}] on the perturbation δa in parameter space can be described as follows

(R′[a_{s}]δa) (x_{r}, t_{r}) = ∫_{Ω} Ψ_{a}(x_{r}, t_{r}; x_{sc})δa(x_{sc})dx_{sc}.  (51)

2.) ’Backprojection step’: The action of the adjoint linearized residual operator R′[a_{s}]^{∗} on the vector ζ in data space can be described as follows

(R_{a}^{′∗}ζ) (x_{sc}) = ∫_{∂Ω} ∫_{[0,T]} Ψ_{a}(x_{r}, t_{r}; x_{sc})ζ(x_{r}, t_{r}) dt_{r}dx_{r}.  (52)

The functions Ψ_{a}(x_{r}, t_{r}; x_{sc}) are called ’absorption transport sensitivity functions’. Here, x_{r} and t_{r} denote receiver location and receiver time, respectively, and x_{sc} is an arbitrary scattering point in the medium.

A derivation of this theorem, and a more detailed discussion of these sensitivity functions, can be found in [14]. We have displayed typical sensitivity functions for our imaging configuration in figure 1. The upper row of the figure shows the distribution of the scattering parameter for two different imaging situations. On the left, we have no clear layers, on the right we have a clear layer close to the boundary of the imaging domain which might represent for example the Cerebro-Spinal Fluid Layer (CFL) of the human head [1, 35, 39]. The source position considered here is indicated as point ’P1’ in the figure on the left hand side. The receiver is located at the point ’P3’. (We are not considering here the also marked point ’P2’.) These points are the same in the situation shown on the right hand side. The second row displays the sensitivity functions with respect to the absorption coefficient for these two situations and for those data which are due to a source position at P1 and taken at the receiver position P3 at an early time, namely t = 10 s (recall the normalization c = 1 cm/s). It can be seen that these early times are very sensitive to the clear layer. Our interpretation of this observation is that the early photons are very likely to have been channeled through this clear layer, and are therefore very sensitive to small variations of the absorption parameter inside the clear layer.

The third row shows the same sensitivity functions, but at a later time t = 24 s. We see that these later times are more sensitive to variations in the absorption parameter inside the domain and less to those inside the clear layer. Photons which arrive that late at the given receiver, are very likely to have travelled a long time in the scattering medium, and are therefore more affected by parameter variations inside the medium. However, we also observe that these later sensitivity functions are more ’spread-out’ than those corresponding to earlier times, since the photons had more time to reach remote areas inside the imaging domain.

Figure 1. Transport absorption sensitivity functions −Ψ_{a}(x_{r}, t_{r}; x_{sc}) for a source located at P 1 and a receiver at P 3. Left column: homogeneous background. Right column: homogeneous background with clear region close to boundary. Upper row: the domains. Middle row: at an early time step (t_{r} = 10 s). Bottom row: at a later time step (t_{r} = 24 s).

Since the sensitivity structure of these later photons is more focused onto the domain of interest, we will use only these data for solving our inversion task. The inversion process incorporates an iterative backprojection (’backtransport’) step along these sensitivity functions, as indicated in theorem 4.4 and explained in [14, 15]. Since the functions along which the data are backprojected are very spread-out, we expect a loss in resolution compared to those reconstructions which can make use of the more focused earlier-time sensitivity functions.

4.8. The algorithm. The nonlinear Kaczmarz-type method for shape reconstruction using level sets is described in brief algorithmic form in table 1. Here, η is a relaxation parameter for the update of the level set function which is determined empirically. The operator Ĉ_{j}^{−1} is taken to be simply a multiple of the identity, such that it can be considered as part of η. The constant C_{ρ}(Γ) could be calculated Figure 2. A typical data curve at one receiver in DOT.

Table 1: The level set reconstruction algorithm.

Initialization.

n = 0;
(D^{(0)}, a^{(0)}, ϕ^{(0)}) given from TBT.

Reconstruction loop.

FOR i = 1 : I_{m} perform I_{m} sweeps FOR j = 1 : p march over source positions

ζ_{j} = Ĉ_{j}^{−1}(M_{j}u_{j} − G̃_{j}); u_{j} solves (1)–(3) with a^{(n)}, q_{j}

δϕ^{(n)} = −(â − a_{b}(x))I_{j}(Λ(ϕ^{(n)}))χ_{B_{ρ}(Γ)}; z_{j} solves (43)-(45) with a^{(n)} and ζ_{j}

ϕ^{(n+1)} = C(n)_{LS}(ϕ^{(n)} + η \frac{C_{ρ}(Γ)}{c_{1}}δϕ^{(n)}); update level set function

a^{(n+1)} = a_{b} + Λ(ϕ^{(n+1)}); n = n + 1; reinitialization n → n + 1 END END (D^{(N)}, a^{(N)}, ϕ^{(N)}) = (D[ϕ^{(n)}], a_{b} + Λ(ϕ^{(n)}), ϕ^{(n)}); final reconstruction.

explicitly for the actual curve Γ^{(n)}, or it could be approximated by some value corresponding to a simple geometrical object (to give an example, in case of a single circle it would be C_{ρ}(Γ) = (2ρ)^{−1}). In our numerical experiments so far, however, it is simply considered as part of η. The same holds true for c_{1}. The constant ρ is in our numerical experiments chosen between one and two grid cells. The scaling factor C(n)_{LS} is determined after each update to keep the global minimum (or maximum) of the level set function at a constant value.

Notice that we typically apply a very small relaxation parameter η, which has the effect that the accumulated contributions of the earlier updates are not totally destroyed by the new correction, but instead have some ’inertia’ which is the larger the smaller we choose η. This slows down the whole reconstruction process to some degree, but, on the other hand, stabilizes it significantly. In our numerical experiments, we typically choose η such that the corrections to the boundaries of the shapes in each update are not much larger than one or two pixels in the normal direction to these boundaries.

5. The TBT-method and the first guess

For starting our shape reconstruction method using level sets we will need an initial guess (D^{(0)}, a^{(0)}, ϕ^{(0)}). Although it is possible just to create an arbitrary initial guess without using any data at all, there are several reasons for employing a preprocessing step which is relatively inexpensive and which yields a good initial guess for the level set reconstruction method. This preprocessing routine usually indicates already the approximate locations of possible inclusions, such that the amount of work which has to be done by the level set method is significantly reduced, and the algorithm converges much faster.

We are employing an inverse scattering algorithm (the TBT-algorithm) for this purpose which was introduced in [9, 10] and which shares basic features with the level set reconstruction scheme. Typically, in the diffusive regime, this inverse scattering algorithm yields already in the early sweeps good low-contrast approximations to the high-contrast objects. These low-contrast approximations usually indicate already the approximate locations of the unknown objects, but typically do not capture the correct number, sizes, shapes or contrasts of these objects. For this, many more iterations of the TBT scheme would be necessary. Therefore, it is advantageous to apply only a few steps of this inverse scattering method, and to switch then to a shape reconstruction formulation which explicitly takes into account the highcontrast assumption. We will show that in fact much better reconstructions can be achieved rapidly in these situations by the shape reconstruction method compared to the inverse scattering algorithm. We refer for a detailed description of the TBTmethod to the publications [9, 10], and describe in the following only the scheme which we use for extracting a suitable initial triple (D^{(0)}, a^{(0)}, ϕ^{(0)}) from the TBT reconstructions a_{TBT}(x).

Let as before D̃ denote the sought domain. We assume that we are working in a high-contrast situation, such that exactly one of the following conditions is satisfied

â − a_{b}(x) ≫ 0 for all x ∈ D̃  (53)
â − a_{b}(x) ≪ 0 for all x ∈ D̃.  (54)

Since we know â and a_{b}(x), we know the constant

sign(D̃) := \{ \begin{matrix}1 , & if (53) holds \\ −1 , & if (54) holds. \end{matrix}  (55)

Choose a threshold value 0 < γ_{LS} < 1 (in our numerical examples presented in section 6 we use γ_{LS} = 0.9) and define

a_{LS} := γ_{LS} max_{x∈Ω} a_{TBT}(x).  (56)

For the level set zero L_{0}^{(0)} of ϕ^{(0)} we require that

L_{0}^{(0)} = {x ∈ Ω : a_{TBT}(x) = a_{LS}} .  (57)

This means that we want all points of Ω where the reconstruction a_{TBT}(x) has exactly the value a_{LS} to be mapped to zero by the level set function ϕ^{(0)}

ϕ^{(0)}(x) = 0 for all x ∈ L_{0}^{(0)}.  (58)

The level set function is now defined as

ϕ^{(0)}(x) = C(0)_{LS}sign(D̃) (a_{LS} − a_{TBT}(x)) ,  (59)

where C(0)_{LS} is some suitably chosen scaling factor. Notice that (58) and (59) are consistent. The initial scattering domain D^{(0)} and the absorption a^{(0)} are defined as

D^{(0)} = D[ϕ^{(0)}],  a^{(0)} = a_{b} + Λ(ϕ^{(0)}).  (60)

Together with ϕ^{(0)} they form an admissible triple (D^{(0)}, a^{(0)}, ϕ^{(0)}).

6. Numerical Experiments

6.1. The computational domain. In our numerical experiments we consider a domain of the size 5×5 cm^{2} which is divided into 50×50 pixels each having a size of 1×1 mm^{2}. One time step lasts 0.2 seconds ([s]), and measurements are taken during the first 100 time steps or 20.0 s. Since the velocity of the light in the medium is normalized to c = 1.0 cm/s, each particle can travel 20 cm in that time. The directional variable is discretized into 12 direction vectors which are equidistantly distributed over the unit circle. The scattering law is a Henyey-Greenstein function (5) with g = 0.9. The space and time derivatives of (1) are discretized by using a finite difference scheme. A more detailed description of our numerical scheme for solving (1)–(3) can be found in [9, 10].

Our sources have a width of 5 pixels or 5 mm along the boundary. Using these slightly extended sources reduces artifacts which tend to appear close to the boundaries, but is not essential for the presented method. In total, there are 16 of these sources, 4 at each side of the domain. These sources are emitting an ultrashort (one time step) laser pulse at time t = 0 in the direction perpendicular to the boundary into the medium. They are arranged in a way such that only the central 20 mm of each side are illuminated during the whole imaging process. In other words, those parts of the boundary whose distance to any of the corners is less than 15 mm are not illuminated in our experiments. The reason for this source arrangement is the location of the clear layer close to the boundary, combined with the fact that our domain is rectangular rather than circular. By this source arrangement, we avoid that ballistic (i.e. unscattered) photons travel along the clear regions.

Receivers are positioned at those parts of the boundary which have a distance (along the boundary) of at least 5 cm or 50 pixels to the source position. Therefore, only half of the boundary is actually covered by the receivers in a given experiment. This arrangement has been chosen since the source-near receivers are not very sensitive to the regions of interest inside the domain. Therefore, they are not supposed to carry much information about the sought objects. In addition, we believe that the approximation of the matrices C_{j,n} by a multiple of the identity is easier to justify with this source-receiver geometry, because of the reduced ’dynamic range’ of the corresponding data. However, we do not have any proof for this assumption so far.

Another important detail is that we are only using those photons as data which arrived at the receiver positions at times between t = 8 s and t = 20 s. In other words, the early photons are not used for the inversion task, since they are more likely to have been channeled through the clear region close to the boundary. We plan to use the information of these early photons in our future work for the additional determination of the structure of the clear layers from the measured data, which is here simply assumed to be known.

The background of our imaging medium is homogeneous except of the embedded clear layers. These homogeneous regions have a scattering parameter of b = 100.0 cm^{−1}, and an absorption parameter of a_{b} = 0.1 cm^{−1}. Inside the clear regions, the parameters are b = 0.01 cm^{−1} and a_{b} = 0.01 cm^{−1}. In the first two numerical experiments, we have only one clear layer close to the boundary and surrounding the interior imaging domain, see figures 1, 3, 4. In the third numerical experiment, we have two additional disc-shaped clear regions imbedded in the interior imaging domain, see figure 5.

6.2. Three numerical experiments. The results of our first numerical experiment are shown in figure 3. We want to recover three disc-shaped objects in the interior of the domain from time-dependent data. The true configuration is shown in the bottom right image of the figure. The values inside the obstacles are a = 0.5 cm^{−1} and b = 100.0 cm^{−1}. The data are created by running our finite differences discretization of the linear transport equation (1)-(3) on the correct model. So far, no noise is added to the data. We plan to investigate the effect of noise in the data more carefully in our future research.

First, we run our preprocessing step (a TBT algorithm) in order to find a good estimate of the number and locations of the obstacles. The results of these runs after 5 sweeps and after 20 sweeps of the TBT method are shown in the first two images of the figure. Notice that one sweep of the TBT method consists of using the data for each of the sources exactly once, in a nonlinear Kaczmarz-type scheme. For more details of that scheme see section 4.4 and [9, 10, 34]. We see that the locations and even the number of the unknown obstacles are already indicated in the reconstructions after 20 sweeps, but that the contrast and the shapes are far away from being correct. We use the result in order to start our level-set based shape reconstruction scheme. The figure shows the initial shape which is extracted from the TBT reconstruction, the reconstruction after six steps of the algorithm (which means after only using the information corresponding to six of the 16 source positions once), and after 16 steps or one sweep of the algorithm. We see that after only one additional sweep, now with the level set based shape reconstruction scheme, we have arrived at a very good estimate of the shapes and locations of the unknown obstacles. Additional iterations (sweeps) do not improve the reconstruction significantly, such that they are not displayed here. In this second step we have assumed that we know the correct contrast of the obstacles to the background. We conclude that the explicit incorporation of this additional information of the contrast in our reconstruction scheme enables us to get rapidly a much better estimate of the shapes of the obstacles than it would have been possible by just continuing the pixel-based TBT reconstruction scheme.

Certainly, the final goal of the reconstruction algorithm will be to find the shapes as well as the contrast of the obstacles simultaneously from the given data. We will address this question in our future research. Notice that the shapes of the Figure 3. The correct value a = 0.5 cm^{−1} is used for the reconstruction. Top row from left to right: a_{TBT} after 5 sweeps of TBT; after 20 sweeps of TBT; initial shape for level set algorithm. Bottom row from left to right: shape after 6 steps of level set algorithm; after 16 steps (i.e. 1 sweep); true object.

obstacles change topology during the evolution process. The level set method has no difficulties in keeping track of these topology changes.

As mentioned, in practical applications we usually do not know the correct values of the absorption parameter inside the shapes. In our second numerical experiment, we investigate the behavior of the level set scheme in the case that we prescribe a slightly incorrect value for the contrast. Instead of a_{s} = 0.5 cm^{−1} inside the obstacle we assume a_{s} = 0.55 cm^{−1}. The new configuration is displayed in the bottom right image of figure 4. The other images of the figure show the results of the algorithm. We see from these results that the effect of this model error on the reconstruction is very small. After only one additional sweep with the level set based reconstruction scheme, we arrive at a very good estimate for the shapes of the obstacles. These shapes also do not differ much from those where the correct value is used. That indicates that the data which are used here are not very sensitive to small changes in the parameter values.

In our third numerical experiment we investigate a slightly more complex situation. Now, there are two additional disc-shaped clear regions in the domain as shown in Figure 5. The bottom right image of the figure shows the true configuration. Although we assume that we know the boundaries of these additional clear regions, such that we do not have to recover them simultaneously with the unknown obstacles, these regions might have a negative impact on our ability to recover these obstacles due to the changed sensitivity structure. We refer to [14] for visualizations of the sensitivity functions in a similar situation. As an additional difficulty, we assume now that each of the three objects which we are looking for has a different contrast to the background. The three values of the absorption parameter inside the objects are a = 0.4 cm^{−1}, a = 0.5 cm^{−1}, and a = 0.6 cm^{−1}. In our reconstruction process, we apply the constant value a = 0.4 cm^{−1} for all three unknown obstacles, Figure 4. An approximate value a = 0.55 cm^{−1} (instead of a = 0.5 cm^{−1}) is used for the reconstruction. Top row from left to right: a_{TBT} after 5 sweeps of TBT; after 20 sweeps of TBT; initial shape for level set algorithm. Bottom row from left to right: shape after 6 steps of level set algorithm; after 16 steps (i.e. 1 sweep); true object.

such that two of the objects are reconstructed with an incorrect value. Figure 5 shows first the reconstruction of the TBT-method after 20 sweeps. This is used for extracting a first guess for our shape reconstruction scheme, as it is shown in the second image on the top row. The third image on the top row shows the reconstruction after 6 steps of the shape reconstruction scheme (which means that only the information corresponding to the first six source position has been used exactly once), after three sweeps, and after 10 sweeps of the scheme. As in our first two numerical examples, the algorithm delivered very good estimates of the shapes after only a few sweeps. The evolution of the norm of the residuals during the reconstruction is displayed in figure 6. The residuals converge now more slowly than in the first two examples, and the later sweeps still yield small improvements to the reconstructions. Notice also that the final reconstructions of the obstacles tend to be slightly more extended than the original objects in the case that a smaller contrast value was used for the inversion. This, of course, might have been expected since a larger area of absorption material is now necessary in order to have the same effect on the data as the original object with a higher absorption value. Finally, we observe also here that the shapes do change topology in order to arrive at the final solution. This is handled easily by the level set formulation.

7. Summary and future research directions

We have presented a novel two-step shape reconstruction method for DOT which uses level sets and adjoint fields. The first step of the method employes the TBT method in order to determine a first guess for the shapes. Then, the second step uses this first guess to start a shape evolution scheme in order to complete the inversion task. The combination of these two schemes makes use of the fact that the TBT method typically finds a good guess for the locations of hidden obstacles in the Figure 5. The smallest (a = 0.4 cm^{−1}) of the three values (a = 0.4 cm^{−1}, a = 0.5 cm^{−1} and a = 0.6 cm^{−1}) is used for the reconstruction. Top row from left to right: a_{TBT} after 20 sweeps of TBT; initial shape for level set algorithm; shape after 6 steps of level set algorithm; Bottom row from left to right: shape after 3 sweeps of level set algorithm; after 10 sweeps; true object.

Figure 6. Norm of the residuals during the level set reconstruction process shown in figure 5

early sweeps, and that the level set based shape reconstruction method is particularly fast once a good approximation for the shapes has been found. Therefore, the combination of these two methods enables us to get very fast and accurate reconstructions of the unknown objects. The reformulation of the inverse problem as a shape reconstruction problem can be considered as a regularization technique. The incorporated prior information about the high contrast of the inclusions and the approximately known values of the parameters stabilizes and speeds up the inversion task. We believe that also the simultaneous reconstruction of shapes and contrast will be possible with a slight modification of the algorithm, and plan to investigate this interesting problem in our future work.

So far, we have been using mainly the information from later arriving photons for the inversion task, in order to cope with the specific structure of the sensitivity functions when clear layers are present. In our future work, the additional information of the early photons will be used in order to try to determine also the structure of the clear layers from the given data. Moreover, we plan to extend the imaging problem to finding the scattering parameter b(x) as well. An important question to be addressed is also the proper incorporation of additional regularization techniques, and the correct choice of the underlying function spaces. Especially the choice of the level set space is not obvious from the beginning. The proper choice of this space will ideally also be connected to our regularization strategy. We conclude that the novel two-step shape reconstruction technique presented here is extremely promising for the application in DOT, and that it poses many interesting and mathematically challenging questions which are open for future research.

References

[1] Arridge S R (1999) Optical tomography in medical imaging Inverse Problems 15 (2), R41–R93 [2] Burger M (2001) A level set method for inverse problems Inverse Problems 17 1327–1355 [3] Case K M and Zweifel P F (1967) Linear Transport Theory (New York: Plenum Press)

[4] Claerbout J F 1976 Fundamentals of Geophysical Data Processing: With Applications to Petroleum Prospecting McGraw-Hill, New York

[5] Dautray R and Lions J L (1993) Mathematical Analysis and Numerical Methods for Science and Technology Vol.6 (Berlin: Springer)

[6] Dierkes T (2000) Rekonstruktionsverfahren zur optischen Tomographie Thesis Preprints Angewandte Mathematik und Informatik 16/00-N, Münster (in german)

[7] Dierkes T, Dorn O, Natterer F, Palamodov V and Sielschott H (2001) Fréchet derivatives for some bilinear inverse problems to appear in ’SIAM Journal on Applied Mathematics’.

[8] Dines K A and Little R J 1979 Computerized geophysical tomography Proc IEEE 67 1065–73 [9] Dorn O (1997) Das inverse Transportproblem in der Lasertomographie Thesis Preprints “Ange-

wandte Mathematik und Informatik” 7/97-N, Münster (in german)

[10] Dorn O (1998) A transport-backtransport method for optical tomography Inverse Problems 14 1107–1130

[11] Dorn O, Bertete-Aguirre H, Berryman J G and Papanicolaou G C 1999 A nonlinear inversion method for 3D electromagnetic imaging using adjoint fields Inverse Problems 15 1523-1558

[12] Dorn O, Miller E L and Rappaport C M 2000 A shape reconstruction method for electromagnetic tomography using adjoint fields and level sets, Inverse Problems 16 1119-1156

[13] Dorn O, Miller E L and Rappaport C M 2001 Shape reconstruction in 2D from limited-view multifrequency electromagnetic data, AMS series Contemporary Mathematics, 278, ’Radon Transform and Tomography’, 97-122

[14] Dorn O 2000 Scattering and absorption transport sensitivity functions for optical tomography Optics Express 7 (13) (Focus Issue: Diffuse Optical Tomography) 492–506, electronically available at http://www.opticsexpress.org/opticsexpress/tocv7n13.htm

[15] Dorn O, Bertete-Aguirre H, Berryman J G and Papanicolaou G C 2001 Sensitivity analysis of a nonlinear inversion method for 3D electromagnetic imaging in anisotropic media to appear in ’Inverse Problems’.

[16] Elangovan V and Whitaker R (2001) From Sinograms to surfaces: a direct approach to the segmentation of tomographic data International Conference on Medical Image Computing and Computer-Assisted Intervention (MICCAI) 2001 Utrecht, The Netherlands

[17] Gaudette J J, Brooks C A, Kilmer M E, Miller E L, Gaudette T and Boas D (2000) A comparison study of linear reconstruction techniques for diffuse optical tomographic imaging of absorption coefficient, Phys.Med.Biol. 45 (4) 1051–1070.

[18] Grünbaum F A (2001) A nonlinear inverse problem inspired by three-dimensional diffuse tomography, Inverse Problems 17 (6) 1907-1922

[19] Harabetian E and Osher S (1998) Regularization of ill-posed problems via the level set approach SIAM J.Appl.Math. 58 (6) 1689–1706

[20] Hielscher A, Klose A and Beuthan J Evolution strategies for optical tomographic characterization of homogeneous media, Optics Express 7 (13) (Focus Issue: Diffuse Optical Tomography) 507

[21] Hintermüller M and Ring W (2001) A level set approach for the solution of a state-constraint optimal control problem Report No. 212 University of Graz and Technical University of Graz, February 2001

[22] Ito K, Kunisch K and Li Z (2001) Level-set function approach to an inverse interface problem Inverse Problems 17 1225–1242

[23] Jackowska-Strumillo L, Sokolowski J, Zochowski A and Henrot A (2001) On numerical solution of shape inverse problems, preprint.

[24] Kaltenbach J M and Kaschke M (1993) Frequency- and Time-Domain Modelling of Light transport in Random Media, in: Medical Optical Tomography ed. Potter (1993) SPIE Optical Engineering Press, Vol. IS11 65-86

[25] Klibanov M V, Lucas T R and Frank R M (1997) A fast and accurate imaging algorithm in optical/diffusion tomography Inverse Problems 13 1341-1361

[26] Kolehmainen V, Arridge S R, Lionheart W R B, Vauhkonen M and Kaipio J P 1999 Recovery of region boundaries of piecewise constant coefficients of an elliptic PDE from boundary data Inverse Problems 15 1375-1391

[27] Kolehmainen V (2001) Novel approaches to image reconstruction in diffusion tomography thesis, Kuopio University Publications C. Natural and Environmental Sciences 125

[28] Litman A, Lesselier D and Santosa F 1998 Reconstruction of a two-dimensional binary obstacle by controlled evolution of a level-set Inverse Problems 14 685-706

[29] Miller E, Kilmer M and Rappaport C 1999 A New Shape-Based Method for Object Localization and Characterization from Scattered Field Data, submitted to IEEE Trans. Geoscience and Remote Sensing

[30] Moscoso M, Keller J B, and Papanicolaou G (2001) Depolarization and blurring of optical images by biological tissues, J. Opt. Soc. Am. A 18 948-960

[31] Natterer F 1986 The Mathematics of Computerized Tomography (Stuttgart: Teubner)

[32] Natterer F and Wübbeling F 1995 A propagation-backpropagation method for ultrasound tomography Inverse Problems 11 1225–1232

[33] Natterer F 1996 Numerical Solution of Bilinear Inverse Problems Preprints “Angewandte Mathematik und Informatik” 19/96-N Münster. This preprint is electronically available at http://wwwmath.uni-muenster.de/math/num/Preprints/1997/natterer 3.

[34] Natterer F and Wübbeling F, Mathematical Methods in Image Reconstruction, Monographs on Mathematical Modeling and Computation 5, SIAM 2001.

[35] E. Okada, M. Firbank, M. Schweiger, S. R. Arridge, M. Cope and D. T. Delpy, “Theoretical and experimental investigation of near-infrared light propagation in a model of the adult head,” Appl. Opt. 36 (1), 21–31 (1997)

[36] Osher S and Sethian J 1988 Fronts propagation with curvature dependent speed: Algorithms based on Hamilton-Jacobi formulations Journal of Computational Physics 56 12-49

[37] Osher J S and Santosa F (2001) Level set methods for optimization problems involving geometry and constraints I. Frequencies of a two-density inhomogeneous drum J. Comp. Phys. 171 272-288

[38] Ramananjaona C, Lambert M, Lesselier D and Zolesio J-P (2001) Shape reconstruction of buried obstacles by controlled evolution of a level set: from a min-max formulation to numerical experiments Inverse Problems 17 1087–1111

[39] J. Riley, H. Dehghani, M. Schweiger, S. R. Arridge, J. Ripoll and M. Nieto-Vesperinas (2000) “3D Optical Tomography in the Presence of Void Regions,” Optics Express 7 (13) (Focus Issue: Diffuse Optical Tomography) 462

[40] Santosa F 1996 A Level-Set Approach for Inverse Problems Involving Obstacles ESAIM: Control, Optimization and Calculus of Variations 1 17-33 [41] Sethian J A 1990 Numerical algorithms for propagating interfaces: Hamilton-Jacobi equations

and conservation laws J. Diff. Geom. 31 131-61

[42] Sethian J A 1999 Level Set Methods and Fast Marching Methods (2nd ed) Cambridge University Press

[43] Sethian J A and Wiegmann A (1999) Structural boundary design via level set and immersed interface methods J. Comp. Phys. 163 (2) 489-528

[44] Sielschott H (2001) Rückpropagationsverfahren für die Wellengleichung in bewegtem Medium thesis, Preprints “Angewandte Mathematik und Informatik” University of Münster, Germany (in german)

[45] Svaiter B F and Zubelli J P (2001) Convexity of the diffuse tomography model Inverse Problems 17 729–738

[46] Zhao H K, Osher S, Merriman B, and Kang M (2000) Implicit and Non-parametric Shape Reconstruction from Unorganized Points Using Variational Level Set Method Computer Vision and Image Understanding 80 295-319

Oliver Dorn, Department of Computer Science, University of British Columbia, 201-2366 Main Mall, Vancouver, B.C. V6T1Z4, Canada

   E-mail address: dorn@cs.ubc.ca