D-76128 KarlsruheGEOMETRIC RECONSTRUCTION IN BIOLUMINESCENCE TOMOGRAPHY

In bioluminescence tomography the location as well as the radiation intensity 
of a photon source (marked cell clusters) inside an organism have to be determined given the outside photon count. This inverse source problem is ill-posed: it suffers not only from strong instability but also from non-uniqueness. 
To cope with these difficulties the source is modeled as a linear combination of indicator functions of measurable domains leading to a nonlinear operator equation. The solution process is stabilized by a Tikhonov like functional which penalizes the perimeter of the domains. For the resulting minimization problem 
existence of a minimizer, stability, and regularization property are shown. 
Moreover, an approximate variational principle is developed based on the calculated domain derivatives which states that there exist smooth almost stationary points of the Tikhonov like functional near to any of its minimizers. This is a crucial property from a numerical point of view 
as it allows to approximate the searched-for domain by smooth domains. Based on the theoretical 
findings numerical schemes are proposed and tested for star-shaped sources in 2D: 
computational experiments illustrate performance and limitations of the considered approach.


Introduction
Bioluminescence tomography (BLT) is a novel technique to image cells in a living organism (in vivo). To this end DNA of a luminescent protein (so-called luciferase) is infiltrated into the target cells (e.g. tumor cells). These cells will emit photons triggered by luciferin which has to be injected prior to imaging, see [5,22]. From the observed photon flux over the organism's surface one has to recover location and intensity of the photon source [4]. Thus, BLT is an inverse source problem.
In this article we work with the simplest mathematical model for BLT which is the diffusion approximation of the radiative transport equation [11]: Let Ω ⊂ R d , d ∈ {2, 3}, be the object (organism) and let u : Ω → R denote the photon density. Then −div(D∇u) + µu = q in Ω, 2D ∂u ∂ν The measurements are described by the boundary condition (2) D ∂u ∂ν =g on ∂Ω.
If not otherwise required, we assume in the following that the absorption coefficient µ ∈ L ∞ (Ω) as well as the diffusion coefficient D ∈ L ∞ (Ω) are bounded away from zero by µ 0 and D 0 , respectively: µ ≥ µ 0 > 0 and D ≥ D 0 > 0 almost everywhere in Ω. The domain Ω is assumed to be convex with a sufficiently smooth boundary ∂Ω (Lipschitz continuous at least). The term g − describes the photon flux penetrating the object and is known in advance. For the sake of simplicity we assume that it vanishes which is the case in most applications.
Date: March 29, 2012. Part of this work was done while the first author was visiting the numerical analysis group of the Department of Mathematics at the University of Iowa. This visit was partly funded by the Karlsruhe House of Young Scientists (KHYS) whose support is greatly acknowledged. The authors thank Weimin Han (University of Iowa) for fruitful discussions.
Subtracting twice the Neumann values in (2) from the boundary condition in (1) we obtain another possibility to model the measurements, namely by the Dirichlet data (3) u = g on ∂Ω with g = g − − 2g. Since it is numerically more stable to evaluate the Dirichlet data than the Neumann values, we use (3) subsequently. As briefly explained above the bioluminescence sources are marked cells. The light intensity of every living cell is determined by the used marker, more precisely by the luciferase, and constant over the cell. Surely we are not able to resolve every cell, but still on a structure, e.g. a tumor, we may assume a constant intensity. Due to dead cells in this structure we do not know the exact strength over it, but it will lie near the intensity of the used cell line. Additionally, the source function vanishes outside of the cell structure. Consequently, we assume that the source function can be modeled by where χ G i is the characteristic function of a measurable domain G i ⊂ Ω and λ i ∈ [λ i , λ i ] = Λ i . The number I is fixed and has to be set in advance. Moreover, we assume G i ⊂ Ω i for an open subset Ω i ⊂ Ω since an a priori knowledge about the location of the sources may be available. Let us use the notations λ = (λ 1 , . . . , λ I ), G = (G 1 , . . . , G I ) and Λ = Λ 1 ×· · ·×Λ I . In order to analyze the BLT problem and develop some reconstruction algorithms in the following chapters we will write it as a nonlinear operator equation F (λ, G) = g. Here the forward operator F is given by F : Λ × L → L 2 (∂Ω), (λ, G) → u| ∂Ω (5) with u denoting the solution of the BVP (1) and L is some appropriate set of I-tuples of subdomains of Ω. Defining the linear and bounded operator A : L 2 (Ω) → L 2 (∂Ω) that maps the source term q to the Dirichlet values of the solution of the BVP (1), we can rewrite F as F (λ, G) = λ i Aχ G i . It can be shown that the operator A is even compact.
So the inverse problem of the bioluminescence tomography under these assumptions can be written as: Problem 1.1. Given the measurements g, find an intensity vector λ ∈ Λ and a tuple of domains G ∈ L such that F (λ, G) = g.
The ill-posedness of Problem 1.1 originates in the compactness of the operator A. Furthermore, it is not uniquely solvable, even for ball-shaped sources, see [21]. Therefore, the problem needs to be regularized to get stable and in special cases unique solutions. In this work we will consider regularization with a total variation (TV) penalty term which will result in smoothing the boundary of the domains G i . In other words we want to minimize the Tikhonov like functional where we write |Dv| for the BV semi-norm for v ∈ BV (R d ) (see e.g. [2] for details). The regularization term in (6) is identical with the perimeter of the domains G i , see e.g. [2], and will be denoted by In the case of a Lipschitz domain G i the perimeter coincides with the (d − 1)-dimensional Hausdorff measure of ∂G i . In [18] a similar approach was used by Ramlau and Ring for Xray computerized tomography and they called the functional of type J α a Mumford-Shah like functional. So we will do in the following.
Let us point out that in the stated framework the source q is essentially the same under changes on a set of measure zero. Also the perimeter is invariant under such alterations [10]. Therefore, it is reasonable to consider equivalence classes of the measurable domains G i , i.e. domains that coincide but on a set of measure zero, rather than an explicit representative.
The first step is to analyze the above stated minimization problem in detail. Existence, stability and regularization results are presented in Section 2. Since the minimization functional is not differentiable with respect to arbitrary domains, we approximate it by smooth domains and develop an approximate variational principle in Section 3. The required domain derivatives are calculated for both the operator F and the perimeter with respect to G in Section 3.1. Since the first operator is linear with respect to λ and the second does not depend on λ, we obtain the derivative of J α immediately.
In Section 4 we develop the theory for star-shaped domains where we can act on the linear space of parameterizations rather than on a set of domains. Similar results as in the case of general domains are presented. Based on these results we will develop a minimization method in Section 5 and present numerical results in Section 6.

Analysis of the Minimization Problem
Let us study in this section the problem of minimizing J α defined in (6) with G i being general measurable subsets of Ω i , i.e.
with L = L Ω 1 × · · · × L Ω I and L Ω i denoting the set of all measurable subsets of Ω i . We proceed similar to [19] where existence, stability and regularization results for the Mumford-Shah approach under an injectivity assumption was proven. However, the BLT forward operator does not satisfy this property (Assumption 3 in [19]). But we think that this assumption can be weakened such that only injectivity needs to hold with respect to span{χ G i | i = 1, . . . , I} for fixed but arbitrary G i = Ω . In this framework the BLT forward operator would fit if I = 1. But in the case I ≥ 2 even this assumption is not satisfied. A counterexample can be constructed in a ball using a ball-shaped source enclosed by a ring-shaped source. Therefore, we present a different way to show existence and stability of the solution of the minimization problem (7). In contrast to [19], we use the constraint on λ to obtain a compactness result in that variable and based on this similar results as in the cited article.
We point out that the following analysis is valid for any operator F : Λ × L → Y that can be written in the form A λ i χ G i , where Y is a Banach space and A a linear and bounded operator from L 2 to Y .

2.2.
Stability. The aim of introducing the regularization term is to stabilize the reconstruction. This is indeed the case for our approach as we will validate in the sequel where we rely on the following lemma taken from [19].
Lemma 2.2. Let g n → g in L 2 as n → ∞ and denote by J n α the functional J α with g replaced by g n . Further, let (λ n , G n ) be a minimizer of J n α over Λ × L. Then there exists a constant C > 0 with Per(G n ) ≤ C for all n.
Theorem 2.3 (Stability). Let g n → g in L 2 as n → ∞ and let (λ n , G n ) minimize Then there exists a subsequence {(λ n k , G n k )} k converging to a minimizer (λ * , G * ) ∈ Λ × L of J α in the sense that Furthermore, every convergent subsequence of {(λ n , G n )} n converges as defined by (8) to a minimizer of J α .
Proof. From Lemma 2.2 we derive the uniform boundedness of the perimeter of G n . As in the proof of Theorem 2.1 we find a subsequence {(λ n k , G n k )} k and a pair (λ * , G * ) such that χ G n k i converges to χ G * i in L 1 as well as λ n k i χ G n k i to λ * i χ G * i in L 2 for every i. It remains to show that the limit is indeed a minimizer of J α . Since the operator A is bounded, we have as k → ∞. Using this convergence, the lower semicontinuity of the perimeter and the minimal property of (λ n k , G n k ) we conclude for any (λ, G) ∈ Λ × L, i.e. the limit (λ * , G * ) is a minimizer of J α .
2.3. Regularization Property. Combining the above ideas of constructing a convergent subsequence with the regularization result from [19] in a straightforward manner we get that the Mumford-Shah like approach is indeed a regularization method.
In addition, let {δ n } n be a positive null sequence and {g n } n such that Then, with the notation of Theorem 2.3, the sequence {(λ n , G n )} of minimizers of J n α(δn) possesses a subsequence converging to (λ + , G + ) which satisfies Furthermore, every convergent subsequence of {(λ n , G n )} n converges in terms of (8) to a pair (λ † , G † ) with property (9).

Approximation by Smooth Domains
To calculate the derivative of J α with respect to the domain, which is essential for a variational principle, we need some smoothness assumptions. These assumptions may be weakened, but to avoid technical difficulties we suppose throughout the following analysis that the coefficients D, µ are continuously differentiable, Ω ⊂ R d , d ∈ {2, 3}, is an open domain with a C 2 −boundary ∂Ω and that We introduce the shorthand notation of the latter relation G ∈ G = G 1 × · · · × G I .
In view of the following lemma, cf. [6, Theorem 5.5, Chapter 3], our smoothness assumption on G appears not to be too restrictive.    [14,20] we consider variations Γ h of the domain Γ ∈ G 0 = {Γ ⊂ Ω | ∂Γ ∈ C 2 } caused by a vector field h ∈ C 2 0 (Ω, R d ): If h is small enough, say if h C 2 < 1/2, then the vector field h is a contraction and thus ϕ = id + h a diffeomorphism on Ω, where id is the identity map. In this case, Γ h ∈ G 0 . By the domain derivative of a mapping Φ : Since the mappings we want to differentiate depend on the intensity λ as well, we will write ∂ Γ Φ := Φ (Γ) for the domain derivative and will replace Γ by the respective component of G.
In [14] the domain derivative of operators involving general boundary value problems were considered. As a special case we obtain the domain derivative of the operator F . For that purpose we introduce some notation: The jump of a function u at ∂Γ is denoted by where the symbols | + and | − indicate the trace of u approaching ∂Γ from the exterior Ω\Γ and the interior Γ, respectively. The term h ν symbolizes the normal component of h, i.e.
Lemma 3.2 (Domain derivative of F ). The derivative of the operator F defined in (5) with respect to the ith domain in direction h ∈ C 2 0 (Ω, R d ) about (λ, G) is given by Proof. See [14, Theorem 2.9].
Remark 3.3. For later use in Section 5 we give the weak formulation of the transmission boundary problem (10): The next step is to calculate the domain derivative of the penalty term, i.e. of the perimeter operator Per : G 0 → R given by (12) Per(Γ) = |D(χ Γ )|.
Since the boundary of Γ is in particular Lipschitz, we obtain by Remark 10 Based on the right identity of (13) and the explanations in [20] we are able to calculate the derivative of the perimeter function Per with respect to the domain.  (12) with respect to the ith domain in direction h ∈ C 2 0 (Ω, R d ) about G is given by where H ∂G i denotes the mean curvature of ∂G i .
3.1.2. The Combined Derivative. As F , see (5), is linear in λ, its partial Fréchet derivative with respect to the intensity in direction k ∈ R I about (λ, G) ∈ Λ × G is given by Combining this with the domain derivative we are able to differentiate the regularization functional J α .
The term u i is the solution of the transmission boundary value problem (10).

An Approximate Variational
Principle. Based on the derivative of the Mumford-Shah like functional on a dense subset, namely the sets of C 2 -domains, we will now present an approximate variational principle which states that near the minimizer the derivative becomes arbitrarily small. For a rigorous formulation and validation of this assertion we will apply and modify findings from [7,8]. Our resulting approximate variational principle will be formulated for a general subspace V of C 2 because later we want to use optimization techniques in a Hilbert space setting.
Let us introduce the following notation: Moreover, we use the norm In addition, let V be a Banach space with Then there exist for every γ ∈ (0, 1 2C ) a vector field v ∈ V and a λ ε ∈ Λ with Proof. Let us consider the ball B 1 for all λ, λ+k ∈ Λ and w, w+h ∈ B 1 2C withh = h•(id+w) −1 . The existence of a (λ ε , v) ∈ Λ×V satisfying the first three estimates (15), (16), and (17) (17). In (19) we set λ = λ ε , w = v and replace (k, h) by t(k, h), t > 0, which yields Dividing by (k,h) R I ×V and recalling the definition ofC finishes the proof.
That can be seen from applying the chain rule to h =h • (id + v).
We will need an estimate on the volume of the symmetric difference of a domain and its perturbed version to prove the main result of this section below.
Lemma 3.8. Let Γ ∈ G 0 be a domain with finite perimeter and h ∈ C 2 0 (Ω, R d ) a vector field with h C 2 sufficiently small. As usual, let Γ h denote the perturbed domain. Then the following estimates hold for the volume of the symmetric difference Proof. Let Γ be the (countable) union of the disjoint connected domains Γ n . For each n we consider the tube T n h with radius h ∞ around the boundary ∂Γ n . Obviously, Γ∆Γ h ⊂ n T n h and thus In [23] an upper bound for the volumes of tubes of type T n h is given: This inequality is sharp if no cross-sections overlap. The constantC Γ n is an invariant of Γ n and is calculated in [3, Corollary 7.5.5] to bẽ where γ n denotes the genus of Γ n . It can be bounded bỹ As the Γ n 's are disjoint, Per(Γ) = n Per(Γ n ). If d = 2 we finally observe that We will now use Lemma 3.1 as well as the previous two to show that we have nearly stationary C 2 -domains near the minimizing domain. Theorem 3.9 (Approximate Variational Principle). Let (λ * , G * ) be a minimizer of J α and λ * an inner point of Λ. In case d = 3, assume that each component of G * is a finite union of disjoint connected domains. Then for any ε > 0 sufficiently small we can find an intensity vector λ ε ∈ Λ and an I-tuple of C 2 -domains G ε satisfying In case d = 3, each componentG ε i is a finite union of disjoint connected domains for ε 1 sufficiently small. Let N be the maximal number of disjoint domains. Due to the continuity of the norm term in J α and due to the above inequalities we get for an ε 2 > 0 getting smaller with ε 1 . Applying Lemma 3.6 with γ = √ ε 2 =: ε 3 we obtain a λ ε ∈ Λ, a C 2 -function h, and the Using Lemma 3.8 and setting C 2 = 0 and C 3 = 8πIN/3 we observe Thus, The right-hand side of the last estimate converges to 0 for ε 1 → 0. Choosing now ε 1 sufficiently small shows the assertion.

Restriction to Star-shaped Domains
In this section we set the stage for the use of optimization methods to solve the minimization problem (7). All usual optimization methods require an underlying linear space which a set of domains does not provide directly. A standard and intuitive way to overcome this issue is working with parameterizations of the boundary. Here we will assume that the domains G i are star-shaped with respect to a known point m i . This presumption may be weakened by describing the boundaries ∂G i by closed curves, but then more effort is needed to prevent self-intersections of the boundary. For an idea of the latter approach see [13].
4.1. The Minimization Problem for Star-shaped Domains. Let the Ω i 's be star-shaped with respect to m i ∈ Ω i . Further, we consider only domains G i with the same property. In other words, we suppose that there exist functions r Ω i ∈ L ∞ (S d−1 ) and points m i ∈ Ω i such that r Ω i (θ)θ + m i , θ ∈ S d−1 , is a parameterization of the boundary ∂Ω i . Furthermore, we restrict our search for the support of the ith source to the set which can be identified with Again we will use the abbreviations For r ∈ R we understand J α (λ, r) to be the value of J α evaluated at (λ, G r ), where G r ∈ L is the tuple of domains represented by r. In the same way we understand expressions like F (λ, r) and Per(r).
With these definitions we are now able to state the minimization problem under consideration Remark 4.1. For ease of presentation and of coding we assume for the following analysis as well as the numerical experiments in Section 6 that all center points m i of the star-shaped domains Ω i are known. Indeed, one can argue to have some estimates of the m i 's from the measurements taken by CCD (charge-coupled device) image sensors, see [4]. In Section 6 we present one experiment where the center point is not known exactly ( Figure 5). However, considering the center points as unkowns is no problem in principle.

4.2.
Analysis of the Reformulated Minimization Problem. Now, as we have an underlying linear structure, we can address the question of convexity of the functional J α . Convexity of the minimization functional is an important property, since then every stationary point is a global minimizer. Unfortunately, J α is non-convex. Indeed, it is possible to construct a counterexample for the easy case that D, µ are constant, I = 1, and the support of the source is a ball. However, we can show that problem (20) possesses a solution relying on techniques used to prove existence for general domains. As in Section 2 we further obtain stability and the regularization property. Proof. Let {(λ n , r n )} n be a minimizing sequence that decays in J α . We denote by G n the tuple of domains parameterized by r n . As in the proof of Theorem 2.1 there exists a subsequence {G n k } converging to G * in the sense that For elements G n k i and G nm i we have the following relation Since the sequence {G n k i } k is convergent, it is especially a Cauchy sequence. Equality (21) reveals that (r n k i ) d is a Cauchy sequence in L 1 as well and therefore convergent. We denote its limit byr i ∈ L 1 and observer i ≥ 0 almost everywhere as {r n i } ⊂ R. The L 1 -convergence implies pointwise convergence almost everywhere, i.e.
i r n k i +r i | : d = 3, and 0 ≤ r n k i ≤ r Ω i , the dominated convergence theorem yields Let now G r * i be the domain parameterized by r * i =r since the limit is unique. Moreover, r * i ∈ R i holds as the set R i is also closed in L 1 . In the same manner as in the proof of Theorem 2.1 we see that (λ * , r * ) is indeed a solution of problem (20).
Combining the techniques of the last proof with the stability and regularization results of Section 2 we receive analogous results for star-shaped domains. Since the proofs are straightforward, we omit them.
Theorem 4.4 (Stability). Let g n → g in L 2 and denote by J n α the functional J α with g substituted by g n . Then the sequence of minimizers (λ n , r n ) of J n α over Λ × R possesses a subsequence converging to a minimizer of J α over Λ × R in R I × L 1 (S d−1 ) I .
Furthermore, every convergent subsequence of {(λ n , r n )} n converges in R I × L 1 (S d−1 ) I to a minimizer of J α .

Approximation by Smooth
Parameterizations. Similar to Section 3 we will develop an approximate variational principle for star-shaped domains. This result will be the justification to use optimization methods that converge to a critical point in the following sections.
To prove the main theorem below we need an approximation result for star-shaped domains analogous to Lemma 3.1.
e. such that the star-shaped domain Γ parameterized by ρ has finite perimeter. Then there exists a sequence {ρ n } n ⊂ C ∞ (S d−1 ) with ρ n − ρ L p → 0 and Per(ρ n ) → Per(ρ) as n → ∞.
Proof. We recall that the perimeter of Γ is given by, cf. [10], Using polar coordinates we observe Let Γ n be the domain parameterized by ρ n . By ν n we denote the unit outward normal of Γ n and ϕ n ∈ C 1 (R d , R d ) is an extension of ν n satisfying ϕ n ∞ ≤ 1. Applying first the dominated convergence theorem and then Gauss's theorem we deduce that Please note that the last equation holds true because Γ n is a smooth domain.
Similar to the proof of Theorem 4.3 we first see that and then conclude that Obviously, it follows that Choosing ε 1 arbitrarily small shows the assertion.

Numerical Schemes
In this section we develop descent methods to minimize J α for star-shaped domains. Since this functional is not differentiable with respect to general domains, we restrict ourselves to a dense subspace U ⊂ C 2 (S d−1 ) I . We assume U to be a Hilbert space. In view of Theorem 4.7 there exist smooth almost stationary points in any neighborhood of a minimizer and we therefore expect a descent method that approaches a stationary point in Λ×U also approaches a minimizer of J α .
Further, we have to implement the constraints λ ∈ Λ and 0 ≤ r i ≤ r Ω i in the optimization process and possibly a boundedness of r i away from zero. The latter property may be necessary to show convergence of the scheme. Therefore, we define the closed and convex subset C := Λ × R ad ⊂ Λ × U ∩ R and denote the convex projection onto C by P C . All schemes we consider to solve min need the gradient of J α as well as the projection operator P C . In a first step we provide these quantities.

Gradient and Projection.
The gradient of J α has to satisfy where J α is known from Theorem 3.5: Obviously, the gradient depends on the choice of the Hilbert space U . We start with calculating the L 2 -gradient. Later U will be chosen to be a periodic Sobolev space H s p where we can work with a Fourier expansion of the parameterization. In this context we will get the H s p -gradient by multiplying the jth Fourier coefficient of the L 2 -gradient by (1 + j 2 ) −s in case d = 2 and by (j + 1/2) −2s in case d = 3. 1 The components of the L 2 -gradient have to satisfy where Φ 1 is the parameterization of the unit ball and gr Φ ρ is the Gramian determinant of the derivative of the parameterization Φ ρ of Γ. In the two-dimensional case the last equality reduces to Herein the L 2 -adjoint operator of ∂ r i F (λ, r) is given by with the solution w of the adjoint boundary value problem −div(D∇w) + µw = 0 in Ω, i.e. of Ω D(∇w · ∇v + µwv) dx + 1 2 ∂Ω wv ds = 1 2 ∂Ω ψv ds for all v ∈ H 1 (Ω). This representation of ∂ r i F (λ, r) * can be seen from according to the weak formulation of the transmission boundary value problem (11). Finally, we derive the projection operator onto the set C. It is well-known that the projection in λ onto the interval Λ = [λ i , λ i ] is The projection in r onto R ad , depends again on the choice of U and cannot be expressed explicitly in general. Since in the numerical experiments the iterates stay in R ad in case of suitable initial values, the projection onto R ad is only of interest from a theoretical point of view.

A Gradient Method.
In [15] the projected gradient method specified in Algorithm 1 is presented for constrained optimization in Hilbert spaces. The step size σ k is chosen by the projected Armijo rule: The largest σ k ∈ { 1 2 n : n ∈ N 0 } is chosen satisfying with some constant γ ∈]0, 1[. Under a Hölder-continuity assumption on the gradient of the minimization functional a convergence result for the projected gradient method under the projected Armijo rule is established in [15]. Though we could only achieve a local Lipschitz-continuity of the gradient on Λ × R ad with R ad = {r ∈ U ∩ R | r i ≥ ε} for any ε > 0, we expect Algorithm 1 to converge also in our setting.

Split
Approach. Ramlau and Ring [18] propose a split approach where first the intensity is minimized while freezing the domain and then the domain is updated using the new intensity. Inspired by them, we split the kth iteration into the following two steps: The step size σ k is chosen as above (projected Armijo rule). This leads to Algorithm 2.
(S3) Set h k r = − grad J α (λ k+1 , r k ) r . (S4) Choose σ k by a projected step size rule such that Let us point out that the optimization problem in step (S2) possesses a solution, since J α is a quadratic function in λ and the set Λ is compact. Standard quadratic programming can be used to solve this problem [17]. However, the solution may not be unique unless the matrix K = Aχ G i , Aχ G j L 2 i,j is positive definite.
In the case I = 1 the optimization problem in (S2) is obviously uniquely solvable. In this situation, similar to the unconstrained case in [18], the split approach can be viewed as a descent method for the reduced functional J α (r) := J α (λ(r), r) with λ(r) := arg min λ∈Λ J α (λ, r), as − grad J α (λ(r), r) r is a descent direction forJ α (r) for every r in the interior of R ad .

Numerical Experiments
In this section we present some numerical experiments of the developed Mumford-Shah like approach for BLT, in order to see if this technique is feasible to reconstruct photon sources or not. For the sake of simplicity we restrict ourselves to the situation where the source term q consists of only one characteristic function: q = λχ G . The more general situation (4) poses no principal problems and corresponding numerical results shall be published elsewhere. Our discretization of r requires a matched discretization of the following quantities: 1. the source term q, i.e. the characteristic function χ G , 2. the L 2 -adjoint of ∂ r F (λ, r), see (25), and 3. the gradient of the perimeter, see (14). Recall that both latter objects appear in the second component of the L 2 -gradient ∂ r J α (λ, r) derived in (24).
In the following we describe in detail how we handle above quantities: 1. Let G M be the star-like domain parameterized by r M . Then we interpolate the characteristic function of G M in the finite element space to obtain the source function q h . Now the FEM solver of MATLAB can be straightforwardly applied to evaluate the forward operator A. the product of the mean curvature and the parameterization, by the trapezoidal rule, but this time with equidistant nodes. This is possible as H ∂G M r M is explicitly known over the interval [0, 2π]. We choose the number of nodes to be greater than max{2M +1, 1/h}. Thus, the error is at least of order h. As mentioned in the previous chapter, we do not implement the projection onto R ad , since for suitable initial values the iterates stay in this set. Only the projection of λ onto Λ is used.
The Hilbert space U is chosen to be H 3 p [0, 2π] ⊂ C 2 p [0, 2π] where the subscript p indicates periodic boundary conditions. So the developed theory is applicable. Hettlich [14] reports only a little difference between numerical simulations in the H s -and in the L 2 -setting. Therefore, we also perform some experiments using the L 2 -gradient. 6.2. Model with one Source. For our computations we use the phantom shown in Figure 1. The phantom has the shape of a circular disk with radius 10 and the origin as midpoint. It consists of four different types of tissue, namely bone (B), heart (H), lung (L), and muscle (M). According to [4], realistic optical parameters for these tissues are with µ being the reduced scattering coefficient. Having µ and µ we can derive the diffusion coefficient by the relation The source is placed around the midpoint (6, −3) and its boundary is parameterized by r(ϑ) = 2 − 0.5 cos ϑ + 0.25 sin ϑ − 0.1 sin(3ϑ) with intensity λ = 1. On a mesh with meshsize 0.2 we produce the synthetic data, whereas the inverse problem is solved on a coarser mesh with h = 0.5 in order to avoid the most obvious inverse crime 3 . By linear interpolation we transform the data from the finer grid to the coarser. The relative interpolation error of about 2.4% may be seen as a 'modeling' error.
The maximal degree M of the trigonometric polynomial is set to 8, since there is no big influence of overestimating the degree of the parameterization. However, the maximal degree M should not be chosen to small, as this can cause loss of details. We choose the regularization parameter α manually by visually inspecting the results. For the intensity λ we allow a variation of 30 % of the exact one, i.e. we set Λ = [0.7, 1.3].
In all experiments we start with initial values λ 0 = 1.1 and r 0 ≡ 2.5. Our termination criterion is taken from [16,Chapter 5.4.1]: In the notation of Algorithm 1 and 2, the gradient iteration is stopped if (h k λ , h k r ) R×U ≤ τ a + τ r (h 0 λ , h 0 r ) R×U and the split approach if h k r U ≤ τ a + τ r h 0 r U . The relative and absolute tolerances are chosen as τ r = τ a = 0.005 for both numerical schemes. Further, the parameter γ in the projected Armijo rule is set to 10 −4 and the step size σ is bounded from below by 2 −10 .
6.2.1. H 3 -setting. In Figure 2 two reconstructions by the gradient method are shown for slightly different regularization parameters. In all our experiments we observe a plateau behavior in the regularization parameter: the reconstruction of (λ, r M ) is pretty much stable over a whole range of α-values. However, at certain tipping points the character of the reconstruction changes dramatically. Such a tipping point behavior is demonstrated in Figure 2.
For a reconstruction using the split approach see Figure 3 (left).
6.2.2. L 2 -setting. Figures 4 (left) and 3 (right) display reconstructions by the gradient method and by the split approach, respectively. One observes that the L 2 -setting leads to a better approximation of the domain than the H 3 -regime. Due to the intrinsic smoothing property of the H 3 -gradient the reconstructed domains in the H 3 -setting resemble circular disks. 3 Still we commit a kind of inverse crime as we use the diffusion model for generating the data and for solving the inverse problem. In future work we plan to obtain the data via the radiative transport equation.  We also perform a numerical experiment where we corrupt the artificial data by 3% relative Gaussian noise with respect to a discrete L 2 (∂Ω)-norm. The reconstruction is shown in Figure 4 (right). The difference to the noise-free reconstruction, Figure 4 (left), is gradual because the regularizing effect of the low degree of r M (M = 8) dominates.
Finally, we come back to Remark 4.1. In our inverse solver we work with the midpoint (5, −2) which is different from the midpoint (6, −3) used for generating the data. The reconstruction in  Figure 5 exhibits the expected behavior: After 100 iterations (left) the size of the reconstructed source is comparable to the size of the original one, but it is located farer away from the boundary. After termination of the method, i.e. after 436 iterations (right) the reconstructed source support lies almost completely in the exact one, though it is smaller due to the penalization of the perimeter. In order to fit the photon flux over the surface, this leads in both cases to an over-estimation of the intensity. We emphasize that the reconstructed intensity coincides with the upper bound of the interval Λ = [0.7, 1.3] where we restrict the intensity to a priori. However, choosing the upper bound larger shows similar behavior. Large regularization parameters lead to small support of the sources with high intensities and small parameters cause lower intensities with larger source supports: this is the non-uniqueness of the BLT inverse source problem [21]. Nevertheless, incorporating precise a priori knowledge about the source, e.g. used marker and cell properties, via Λ and α into the reconstruction process will lead to useful results.

Outlook
As the diffusion approximation is only a simplified model of the propagation of photons in tissue, the natural next step is to extend the stated framework to the more realistic model based on the radiative transfer equation. Also in this setting the theory of Section 2 is directly applicable. By contrast, more work has to been done to obtain results as in Section 3. In particular the domain derivative for the forward operator based on the radiative transfer equation has to be developed, which is under investigations right now.