Geometric Reconstruction in Bioluminescence Tomography (BLT) - - PowerPoint PPT Presentation

Geometric Reconstruction in Bioluminescence Tomography (BLT)

Andreas Rieder

jointly with Tim Kreutzmann

KIT – University of the State of Baden-W¨ urttemberg and National Research Center of the Helmholtz Association

FAKULT ¨ AT F ¨ UR MATHEMATIK – INSTITUT F ¨ UR ANGEWANDTE UND NUMERISCHE MATHEMATIK

www.kit.edu (Wang et al. ’06)

Overview

2 / 32 c Andreas Rieder – Geometric Reconstruction in Bioluminescence Tomography CBM 2013, Rio de Janeiro

Mathematical model Inverse problem: formulation & uniqueness Inverse problem: reformulation & stabilization Gradient of the minimization functional Numerical experiments in 2D Summary

Mathematical model

⊲

Mathematical model Inverse problem: formulation & uniqueness Inverse problem: reformulation & stabilization Gradient of the mini- mization functional Numerical experi- ments in 2D Summary

3 / 32 c Andreas Rieder – Geometric Reconstruction in Bioluminescence Tomography CBM 2013, Rio de Janeiro

The (stationary) radiative transfer equation (RTE)

(Boltzmann transport equation)

4 / 32 c Andreas Rieder – Geometric Reconstruction in Bioluminescence Tomography CBM 2013, Rio de Janeiro

Let u(x, θ) be the photon flux (radiance) in direction θ ∈ S2 about x ∈ Ω ⊂

• R3. Then,

θ · ∇u(x, θ) + µ(x)u(x, θ) = µs(x)

• S2

η(θ · θ′)u(x, θ′)dθ′ + q(x, θ) u(x, θ) = g−(x, θ), x ∈ ∂Ω, n(x) · θ ≤ 0 g(x) = 1 4π

• S2

n(x) · θu(x, θ)dθ, x ∈ ∂Ω where µ = µs + µa and µs / µa scattering/absorption coefficients η scattering kernel (

• S2 η(θ · θ′)dθ′ = 1)

q source term µs = 0: RTE yields integral eqs. of transmission and emission tomography

(F. Natterer & F. W¨ ubbeling, Math. Methods in Image Reconstr., SIAM, ’01)

Diffusion approximation: setting

5 / 32 c Andreas Rieder – Geometric Reconstruction in Bioluminescence Tomography CBM 2013, Rio de Janeiro

Assume that u(x, θ) = u0(x) + 3θ · u1(x) where u0(x) = 1 4π

• S2

u(x, θ)dθ ∈ R and u1(x) = 1 4π

• S2

θu(x, θ)dθ ∈ R3. By the Funk-Hecke theorem,

• S2

θη(θ · θ′)dθ = η θ′ where η =

• S2

θ′ · θ η(θ · θ′)dθ is the scattering anisotropy.

Diffusion approximation: derivation

6 / 32 c Andreas Rieder – Geometric Reconstruction in Bioluminescence Tomography CBM 2013, Rio de Janeiro

Integrate RTE over S2, multiply RTE by θ, integrate again, and assume g−(x, θ) = g−(x). Then, −∇ · (D∇u0) + µau0 = q0 := 1 4π

• S2

q(·, θ)dθ, u0 + 2D∂nu0 = g− on ∂Ω, D∂nu0 = −g on ∂Ω, where D = 1 3(µ − ηµs) is the diffusion coefficient (reduced scattering coefficient).

Diffusion approximation: final equation

7 / 32 c Andreas Rieder – Geometric Reconstruction in Bioluminescence Tomography CBM 2013, Rio de Janeiro

Change of notation: u = u0, q = q0, µ = µa, and g = −g. The photon density u obeys the BVP −∇ · (D∇u) + µu = q in Ω, u + 2D∂nu = g− on ∂Ω. The measurements are given by D∂nu = g

• n ∂Ω.

Assume g− = 0 (no photons penetrate the object from outside).

Inverse problem: formulation & uniqueness

Mathematical model

⊲

Inverse problem: formulation & uniqueness Inverse problem: reformulation & stabilization Gradient of the mini- mization functional Numerical experi- ments in 2D Summary

8 / 32 c Andreas Rieder – Geometric Reconstruction in Bioluminescence Tomography CBM 2013, Rio de Janeiro

Inverse problem of BLT (in the diffusive regime)

9 / 32 c Andreas Rieder – Geometric Reconstruction in Bioluminescence Tomography CBM 2013, Rio de Janeiro

Define the (linear) forward operator A: L2(Ω) → H− 1

2 (∂Ω),

q → D∂nu , where u solves the BVP with g− = 0: −∇ · (D∇u) + µu = q in Ω, u + 2D∂nu = 0

• n ∂Ω.

BLT Problem: Given g ∈ R(A), find a source q ∈ L2(Ω) satisfying Aq = g.

Null Space of A

10 / 32 c Andreas Rieder – Geometric Reconstruction in Bioluminescence Tomography CBM 2013, Rio de Janeiro

Lemma (Wang, Li & Jiang ’04, Kreutzmann ’13): There is an isomorphism Φ: H1(Ω) → H1(Ω)′ such that N(A) = Φ

• H1

0(Ω)

• ∩ L2(Ω).

If D ∈ W 1,∞ then N(A) = Φ

• H1

0(Ω) ∩ H2(Ω)

• .

Proof: Define Φ: H1(Ω) → H1(Ω)′, u → (Φu)(v) = a(u, v) where a(u, v) =

• Ω
• D∇u · ∇v + µuv
• dx + 1

2

• ∂Ω

uv ds.

Singular Functions of A: L2(Ω0) → L2(∂Ω)

11 / 32 c Andreas Rieder – Geometric Reconstruction in Bioluminescence Tomography CBM 2013, Rio de Janeiro

H L L M B

−8 −6 −4 −5 5 σ1 = 0.34221 −0.4 −0.3 −0.2 −0.1 −8 −6 −4 −5 5 σ2 = 0.2932 −0.2 0.2 −8 −6 −4 −5 5 σ3 = 0.23248 −0.2 0.2 −8 −6 −4 −5 5 σ4 = 0.17673 −0.4 −0.2 0.2 −8 −6 −4 −5 5 σ5 = 0.1296 −0.2 0.2 0.4 −8 −6 −4 −5 5 σ6 = 0.096659 −0.4 −0.2 0.2 −8 −6 −4 −5 5 σ7 = 0.070445 −0.2 0.2 −8 −6 −4 −5 5 σ8 = 0.04619 −0.2 0.2 −8 −6 −4 −5 5 σ9 = 0.035356 −0.2 0.2 0.4 −8 −6 −4 −5 5 σ10 = 0.021841 −0.2 0.2 0.4 −8 −6 −4 −5 5 σ11 = 0.013265 −0.4 −0.2 0.2 0.4 −8 −6 −4 −5 5 σ12 = 0.0089158 −0.4 −0.2 0.2

Can we restore uniqueness by a priori information?

12 / 32 c Andreas Rieder – Geometric Reconstruction in Bioluminescence Tomography CBM 2013, Rio de Janeiro

Consider, for instance, q = λχS where λ ≥ 0 is a constant and S ⊂ Ω.

Can we restore uniqueness by a priori information?

12 / 32 c Andreas Rieder – Geometric Reconstruction in Bioluminescence Tomography CBM 2013, Rio de Janeiro

Consider, for instance, q = λχS where λ ≥ 0 is a constant and S ⊂ Ω. Lemma (Wang, Li & Jiang ’04): There exist z ∈ Ω, λ1 = λ2 and r1 = r2 such that A(λ1χB1) = A(λ2χB2) with Bk = Brk(z).

Inverse problem: reformulation & stabilization

Mathematical model Inverse problem: formulation & uniqueness

⊲

Inverse problem: reformulation & stabilization Reformulation Tikhonov-like regularization Existence of a minimizer & stability Regularization property Gradient of the mini- mization functional Numerical experi- ments in 2D Summary

13 / 32 c Andreas Rieder – Geometric Reconstruction in Bioluminescence Tomography CBM 2013, Rio de Janeiro

Reformulation

14 / 32 c Andreas Rieder – Geometric Reconstruction in Bioluminescence Tomography CBM 2013, Rio de Janeiro

Ansatz: q =

I

• i=1

λiχSi where Si ⊂ Ω, λi ∈ [λi, λi] = Λi, and I ∈ N. For the ease of presentation: I = 1. Define the nonlinear operator F : Λ × L → L2(∂Ω), (λ, S) → D∂nu|∂Ω where L is the set of all measurable subsets of Ω. Note: F(λ, S) = λAχS BLT Problem: Given measurements g, find an intensity λ ∈ Λ and a domain S ∈ L such that F(λ, S) = g.

Tikhonov-like regularization

15 / 32 c Andreas Rieder – Geometric Reconstruction in Bioluminescence Tomography CBM 2013, Rio de Janeiro

Minimize Jα(λ, S) = 1 2F(λ, S) − g2

L2 + αPer(S) over Λ × L

where α > 0 is the regularization parameter and Per(S) is the perimeter

• f S:

Per(S) = |D(χS)|, with |D(·)| denoting the BV-semi-norm (Ramlau & Ring ’07, ’10).

AG Sahin, Univ. Mainz

Existence of a minimizer & stability

16 / 32 c Andreas Rieder – Geometric Reconstruction in Bioluminescence Tomography CBM 2013, Rio de Janeiro

Theorem: For all α > 0 and g ∈ L2(∂Ω) there exists a minimizer (λ∗, S∗) ∈ Λ × L, that is, Jα(λ∗, S∗) ≤ Jα(λ, S) for all (λ, S) ∈ Λ × L. Theorem: Let gn → g in L2 as n → ∞ and let (λn, Sn) minimize Jn

α(λ, S) = 1 2F(λ, S) − gn2 L2 + αPer(S) over Λ × L.

Then there exists a subsequence {(λnk, Snk)}k converging to a minimizer (λ∗, S∗) ∈ Λ × L of Jα in the sense that λnkχSnk − λ∗χS∗L2 → 0 as k → ∞. Furthermore, every convergent subsequence of {(λn, Sn)}n converges to a minimizer of Jα.

Regularization property

17 / 32 c Andreas Rieder – Geometric Reconstruction in Bioluminescence Tomography CBM 2013, Rio de Janeiro

Theorem: Let g be in range(F) and let δ → α(δ) where α(δ) → 0 and δ2 α(δ) → 0 as δ → 0. In addition, let {δn}n be a positive null sequence and {gn}n such that gn − gL2 ≤ δn. Then, the sequence {(λn, Sn)} of minimizers of Jn

α(δn) possesses a sub-

sequence converging to a solution (λ+, S+) where S+ = arg min{Per(S) : S ∈ L s.t. ∃λ ∈ Λ with F(λ, S) = g}. Furthermore, every convergent subsequence of {(λn, Sn)}n converges to a pair (λ†, S†) with above property.

Gradient of the minimization functional

Mathematical model Inverse problem: formulation & uniqueness Inverse problem: reformulation & stabilization

⊲

Gradient of the minimization functional Domain derivative: general definition Domain derivative

• f F (λ, ·): S →

L2(∂Ω) Domain derivative of Per: S → R Derivative of Jα : Λ × S → R Approximate vari- ational principle (Ekeland 1974) Numerical experi- ments in 2D Summary

18 / 32 c Andreas Rieder – Geometric Reconstruction in Bioluminescence Tomography CBM 2013, Rio de Janeiro

Domain derivative: general definition

19 / 32 c Andreas Rieder – Geometric Reconstruction in Bioluminescence Tomography CBM 2013, Rio de Janeiro

Let Γ ∈ S = { Γ ⊂ Ω : ∂ Γ ∈ C2} and let h ∈ C2

0(Ω, Rd). Define

Γh = {x + h(x) : x ∈ Γ}. If h is small enough, say if hC2 < 1/2, then Γh ∈ S. By the domain derivative of Φ: S → Y about Γ we understand Φ′(Γ) ∈ L(C2, Y ) satisfying Φ(Γh) − Φ(Γ) − Φ′(Γ)hY = o(hC2) where Y is a normed space.

Domain derivative of F(λ, ·): S → L2(∂Ω)

20 / 32 c Andreas Rieder – Geometric Reconstruction in Bioluminescence Tomography CBM 2013, Rio de Janeiro

Reminder: F(λ, S) = λAχS Lemma: We have that ∂SF(λ, S)h = u′|∂Ω where u′ ∈ H1(Ω\∂S) solves the transmission bvp −∇ · (D∇u′) + µu′ = 0 in Ω\∂S, 2D∂nu′ + u′ = 0

• n ∂Ω,

[u′]± = 0

• n ∂S,
• D∂nu′

± = −λh · n

• n ∂S.

Proof: similar to Hettlich’s habilitation thesis 1999.

Domain derivative of Per: S → R

21 / 32 c Andreas Rieder – Geometric Reconstruction in Bioluminescence Tomography CBM 2013, Rio de Janeiro

Lemma (Simon 1980): We have that ∂SPer(S)h =

• ∂S

H∂S(h · n) ds where H∂S denotes the mean curvature of ∂S.

Derivative of Jα : Λ × S → R

22 / 32 c Andreas Rieder – Geometric Reconstruction in Bioluminescence Tomography CBM 2013, Rio de Janeiro

Jα(λ, S) = 1 2F(λ, S) − g2

L2 + αPer(S)

∂λF(λ, S)k = kAχS = F(k, S) Theorem: We have that J′

α(λ, S)(k, h) =

• F(λ, S) − g, F(k, S) + u′

L2(∂Ω) + α

• ∂S

H∂S(h · n) ds for k ∈ R, h ∈ C2

0(Ω, R3).

Proof: J′

α(λ, S)(k, h) = ∂λJα(λ, S)k + ∂SJα(λ, S)h

Approximate variational principle (Ekeland 1974)

23 / 32 c Andreas Rieder – Geometric Reconstruction in Bioluminescence Tomography CBM 2013, Rio de Janeiro

There exist smooth almost critical points of Jα near to any of its minimizers.

Approximate variational principle (Ekeland 1974)

23 / 32 c Andreas Rieder – Geometric Reconstruction in Bioluminescence Tomography CBM 2013, Rio de Janeiro

There exist smooth almost critical points of Jα near to any of its minimizers. Theorem: Let (λ∗, S∗) be a minimizer of Jα where λ∗ is an inner point of Λ. Then, for any ε > 0 sufficiently small there is a (λε, Sε) ∈ Λ × S with Jα(λε, Sε) − Jα(λ∗, S∗) ≤ ε, λεχSε − λ∗χS∗L1 ≤ ε, J′

α(λε, Sε)R×C2→R ≤ ε.

Approximate variational principle (Ekeland 1974)

23 / 32 c Andreas Rieder – Geometric Reconstruction in Bioluminescence Tomography CBM 2013, Rio de Janeiro

There exist smooth almost critical points of Jα near to any of its minimizers. Theorem: Let (λ∗, S∗) be a minimizer of Jα where λ∗ is an inner point of Λ. Then, for any ε > 0 sufficiently small there is a (λε, Sε) ∈ Λ × S with Jα(λε, Sε) − Jα(λ∗, S∗) ≤ ε, λεχSε − λ∗χS∗L1 ≤ ε, J′

α(λε, Sε)R×C2→R ≤ ε.

Proof: Key ingredient is To any bounded measurable Γ ⊂ Rd with finite perimeter exists a se- quence {Γn}n of C∞-domains such that

• Rd |χΓn − χΓ|dx → 0

and Per(Γn) → Per(Γ) as n → ∞.

Numerical experiments in 2D

Mathematical model Inverse problem: formulation & uniqueness Inverse problem: reformulation & stabilization Gradient of the mini- mization functional

⊲

Numerical experiments in 2D Star-shaped do- mains Algorithm: Projected Gradient Method The model H3-Reconstructions L2-Reconstructions Summary

24 / 32 c Andreas Rieder – Geometric Reconstruction in Bioluminescence Tomography CBM 2013, Rio de Janeiro

Star-shaped domains

25 / 32 c Andreas Rieder – Geometric Reconstruction in Bioluminescence Tomography CBM 2013, Rio de Janeiro

For the numerical experiments we consider a star-shaped domain only: S = {x ∈ R2 : x = m + t θ(ϑ)r(ϑ), 0 ≤ t ≤ 1, 0 ≤ ϑ ≤ 2π} where m is the center (assumed to be known) and r: [0, 2π] → [0, ∞[ parameterizes the boundary of S. All previous results hold in this setting as well if we work in a space of smooth parameterizations, say, r ∈ H3

p(0, 2π) ⊂ C2 p(0, 2π).

(λ, S) (λ, r) ∈ Λ × Rad where Rad =

• r ∈ H3

p(0, 2π) : r ≥ 0

• .

Gradient equation:

• gradJα(λ, r), (k, h)
• R×H3 = J′

α(λ, r)(k, h).

We have implemented star-shaped domains using trigonometric poly- nomials.

Algorithm: Projected Gradient Method

26 / 32 c Andreas Rieder – Geometric Reconstruction in Bioluminescence Tomography CBM 2013, Rio de Janeiro

(S0) Choose (λ0, r0) ∈ C := Λ × Rad, k := 0 (S1) Iterate (S2)-(S4) until (S2) Set ∇k := −gradJα(λk, rk). (S3) Choose σk by a projected step size rule such that Jα

• PC
• (λk, rk) + σk∇k
• < Jα(λk, rk).

(S4) Set (λk+1, rk+1) := PC

• (λk, rk) + σk∇k
• .

The model

27 / 32 c Andreas Rieder – Geometric Reconstruction in Bioluminescence Tomography CBM 2013, Rio de Janeiro

lung lung heart bone muscle source

H3-Reconstructions

28 / 32 c Andreas Rieder – Geometric Reconstruction in Bioluminescence Tomography CBM 2013, Rio de Janeiro

0.5 1 1.5 2 2.5 3 30 210 60 240 90 270 120 300 150 330 180 λ = 0.97843 for α = 0.00763 0.5 1 1.5 2 2.5 3 30 210 60 240 90 270 120 300 150 330 180 λ = 0.702 for α = 0.00762

Reconstructions (blue) and source (red). Left: 69 iterations, right: k = 17 iterations

L2-Reconstructions

29 / 32 c Andreas Rieder – Geometric Reconstruction in Bioluminescence Tomography CBM 2013, Rio de Janeiro

0.5 1 1.5 2 2.5 3 30 210 60 240 90 270 120 300 150 330 180 λ = 0.84033 for α = 0.0079 0.5 1 1.5 2 2.5 3 30 210 60 240 90 270 120 300 150 330 180 λ = 0.77183 for α = 0.008

Reconstructions (blue) and source (red). Left: 37 iterations, right: noisy data (rel. 3%), 24 iterations

L2-Reconstruction with variable center

30 / 32 c Andreas Rieder – Geometric Reconstruction in Bioluminescence Tomography CBM 2013, Rio de Janeiro

Summary

Mathematical model Inverse problem: formulation & uniqueness Inverse problem: reformulation & stabilization Gradient of the mini- mization functional Numerical experi- ments in 2D

⊲ Summary

What to remember from this talk

31 / 32 c Andreas Rieder – Geometric Reconstruction in Bioluminescence Tomography CBM 2013, Rio de Janeiro

What to remember from this talk

32 / 32 c Andreas Rieder – Geometric Reconstruction in Bioluminescence Tomography CBM 2013, Rio de Janeiro

Bioluminescence tomography images cells in vivo. From a mathemati- cal point of view it is an inverse source problem which suffers from non- uniqueness (diffusion approximation) and ill-posedness. To overcome these difficulties the sources are modeled as ”hot spots” leading to a nonlinear problem which is stabilized by a Tikhonov-like regularization penalizing the perimeter of the hot spots. The approximate variational principle justifies the restriction to hot spots with smooth boundaries. For star-shaped domains in 2D a projected steepest decent solver has been implemented and tested.