- Email: [email protected]
A new total variational regularization method for nonlinear inverse problems in fluorescence molecular tomography
Journal of Computational and Applied Mathematics 365 (2020) 112408
Contents lists available at ScienceDirect
Journal of Computational and Applied Mathematics journal homepage: www.elsevier.com/locate/cam
A new total variational regularization method for nonlinear inverse problems in fluorescence molecular tomography✩ ∗
Li Li a , , Chen Chen a , Bo Bi b a b
School of Mathematics, Harbin Institute of Technology, Harbin, 150001, PR China School of Mathematics and Statics, Northeast Petroleum University, Daqing, 163318, PR China
article
info
Article history: Received 15 November 2018 Received in revised form 31 May 2019 Keywords: Ill-posed problems Total variational regularization Bregman distance Homotopy
a b s t r a c t Fluorescence molecular tomography (FMT) is an imaging way of in vivo optical imaging, which is widely used in biomedical photon imaging for its higher sensitivity. For years, the propagation process of photons in tissue is simulated based on the equation of radiative transfer. In this work, we analyze an inverse problem in FMT, which uses the reading of body surface detector to reconstruct the optical parameter distribution. We propose a novel total variational regularization method to solve the absorption parameter inversion problem in FMT. The proof of the proposed method is also provided. Finally, numerical experiments demonstrate the feasibility of the proposed method. © 2019 Elsevier B.V. All rights reserved.
1. Introduction Since the classical medical tomography depends highly on the physical properties and physiological change of an organization, the lesion of tissue cannot be discovered timely, which delays the treatment of the patient. Biomedical photon imaging is a non-invasive imaging technique, which can overcome the weakness of traditional medical imaging technology and satisfy the demand of medical diagnosis for comfort and nondestructive nature through detecting minute changes in tissue. Diffusion optical tomography (DOT), Fluorescence molecule tomography (FMT), and Bioluminescence tomography (BLT) are some typical biomedical photon imaging technology. In these imaging techniques, FMT can be traced back to 1980s, and is important for its simplicity, low cost, and ability to monitor an organization for a long time. Moreover, this kind of imaging result can satisfy our need for quantitative analysis. FMT was formally produced in 1990s, which obtains the high resolution images by inversing the optical parameters of biological tissue, such as absorption coefficient and reflection coefficient [1]. Moreover, the feasibility of this technology has been proved in the experiments on animals [2]. FMT is used to reconstruct the distribution of fluorescence concentration according to the fluorescence characteristics of near-infrared fluorescent chromosomes. Normally the radiative transport equation (RTE), which is used to describe the propagation or radiation of particles inside a medium, can be used in reconstruction of FMT as a forward model. If the function spaces are chosen appropriately, the solvability of the RTE and the properties of the solutions are all intensively studied [3–6]. Analytical method and numerical method are two main methods for solving RTE equation. The analytical method is considered to solve the problems when the shape of the tissue is regular. Cong considered the tissue as an infinite medium, and the light source as a regular geometry, then proposed an algorithm to obtain the analytic ✩ This work is supported by a scholarship from the China Scholarship Council (No. 201806125006). ∗ Corresponding author. E-mail address: [email protected] (L. Li). https://doi.org/10.1016/j.cam.2019.112408 0377-0427/© 2019 Elsevier B.V. All rights reserved.
2
L. Li, C. Chen and B. Bi / Journal of Computational and Applied Mathematics 365 (2020) 112408
solution [7]. However, the numerical method is more preferred to the situation that the medium is irregular, such as FEM [8], SOR iteration [9], Gauss Seidel iteration, and random statistical method [10–14], et al. Since RTE-based inverse problems are difficult either in theoretical analysis or in numerical computation, more investigations prefer to the diffusion approximation of the RTE (DA) [15], which is first-order phase approximation of RTE and is used to simulate the light migration in photoacoustic imaging. Based on gradient, Klose et al. developed an iterative algorithm to reconstruct distribution of tissue parameters by using the gradient information of operator F about the optical parameter [16]. Similarly, Kui Ren used the following relation
∂p |(0,T )×∂ Ω = Λ(c , σ , g) ∂n to reconstruct (c , σ ) in one-step reconstruction algorithm [17]. Compared with the traditional Newton type algorithm, the method based on gradient does not need to calculate the Jacobi matrix many times, and avoids the repeated matrix inversion, which greatly simplifies the calculation process [9]. Generally speaking, the inverse problems of FMT are ill-posed due to its high degree of absorption and scattering of light through tissue. Effective regularization methods should be utilized to deal with the ill-posedness. Tikhonov regularization [18], Landweber iteration [19] and Levenberg–Marquardt iteration [20], et al. are well-known widespread regularization for solving ill-posed problems. Here, L2 norm is taken as penalty term, and is feasible for solving the problem of smoothness. However, properties of some non-smooth edges are neglected. Compared with L2 norm, total variational norm preserves local smoothness and piecewise constant in the reconstructed images. Therefore, the regularization with the penalty of total variational norm can achieve better reconstructions compared to some classical iteration methods with the penalty term of L2 norm. Gao and Zhao investigated the regularization with both l1 and total-variational norm for bioluminescence tomography based on radiative transfer equation [21]. In 2010, Hyde et al. proposed to solve FMT problems by using two step regularization methods [22]. Thereafter, Dutta et al. introduced the total variational regularization method with l1 penalty, and put forward a new algorithm to deal with the problem of sparse images reconstructed by FMT in the local smooth high-density region [23], which combined with sparsity constrained method and smoothness condition. Joyita et al. proposed an approach involving combination of L1 and total variation norm penalties, which could suppress spurious background signals, enforce sparsity, and preserve local smoothness and piecewise constancy in the reconstructed images [24]. In 2015, B. Bi et al. investigated the feasibility of the sparsity reconstruction approach for DOT based on RTE [25]. Two years later, combined with total variation penalty, B. Bi and J. Tang studied the L∞ norm fitting method for DOT based on RTE [26]. Based on the properties of preserving shapes of total variational regularization and the wide convergent range of homotopy method, this article constructs a new iterative algorithm for solving this problem. The rest of the paper is organized as follows. Firstly, we prove that the error of the new method decrease monotonically in Bregman distance. Secondly, the new method with the noisy data is proved to be a regularization method theoretically. Finally, in the section of numerical simulation, the propagation of photons in tissues is simulated by using RTE, and the new iteration is applied to inverse the absorption parameters in FMT . 2. New total variational regularization and its convergence analysis Considering nonlinear operator equation: F (u) = z
(2.1)
where F is a nonlinear operator, F : D ⊂ L (Ω ) → Z . F is continuous, weakly sequenced closed, and Fréchet differential. u and z are the identified parameter and the observed data, respectively. Usually, we cannot obtain the exact data z, but the noisy data z δ satisfies ∥z δ − z ∥ ≤ δ , where δ is the error level. In fact, the problem (2.1) with noisy data is 2
F (u) = z δ
(2.2)
In order to solve (2.2), we consider the following least-square functional: min J(u) = min
u∈D (F )
u∈D (F )
1 2
∥F (u) − z δ ∥2
(2.3)
Its corresponding Euler equation is F ′ (u)∗ (F (u) − z δ ) = 0
(2.4)
Introduce the parameter t, we construct the homotopy: H(u, t) = t [F ′ (u)∗ (F (u) − z δ )] + (1 − t)(u − u0 ) = 0
(2.5)
u0 is an initial value. As t = 0, the solution of Eq. (2.5) is u = u0 . As t = 1, the solution of Eq. (2.5) is the solution of Euler equation (2.4). If the initial value is known, we can find the final solution of Eq. (2.2) along the homotopy curve.
L. Li, C. Chen and B. Bi / Journal of Computational and Applied Mathematics 365 (2020) 112408
3
[0, 1] is divided into 0 = t0 < t1 < · · · < tl = 1. t = tk (0 ≤ k ≤ l) is considered. Assume the iteration has reached ujk , j
and the higher-order term is ignored. We linearize F (u) at the point uk : j
j
j
F (u) ≈ F (uk ) + F ′ (uk )(u − uk )
(2.6) j+1
Substituting (2.6) into (2.5), and taking the computed u as uk , j+1
uk
= ujk − [tk F ′ (ujk )∗ F ′ (ujk ) + (1 − tk )I ]∗ [tk F ′ (ujk )∗ (F (ujk ) − z δ ) + (1 − tk )(ujk − u0 )]
(2.7)
In the following, the total variational regularization is introduced. Based on Bregman distance, a new iterative algorithm for solving nonlinear ill-posed problems is constructed:
⎧ j+1 j j j j ⎪ uk = arg minu∈D(F ) {⟨[tk F ′ (uk )∗ F ′ (uk ) + (1 − tk )I ]∗ [tk F ′ (uk )∗ (F (uk ) − z δ ) ⎪ ⎪ ⎪ ⎪ ⎨ + (1 − tk )(ujk − u0 )], u − ujk ⟩ + αkj DJ j (u, ujk )} ξk
(2.8)
⎪ ⎪ ξkj+1 = ξkj − (αkj )−1 [tk F ′ (ujk )∗ F ′ (ujk ) + (1 − tk )I ]∗ [tk F ′ (ujk )∗ (F (ujk ) − z δ ) ⎪ ⎪ ⎪ ⎩ + (1 − tk )(ujk − u0 )] J
j
j
Here D j (u, uk ) is defined as the Bregman distance between u and uk : ξk J
j
j
j
j
D j (u, uk ) = J(u) − J(uk ) − ⟨ξk , u − uk ⟩ ξk j
j
J
j
J
j
j
where ξk ∈ ∂ J(uk ) is subgradient. Obviously, D j (u, uk ) ≥ 0, D j (u, uk ) = 0 if and only if u = uk . J(·) is defined as ξk
ξk
J(·) = J1 (·) + J2 (·) where J1 (·) =| · |BV (Ω ) +χD(F ) (·),
J 2 (· ) =
κ 2
∥ · ∥2L2 (Ω )
Suppose D(F ) is convex, then J(·) and J1 (·) are convex. Any ξ ∈ ∂ J(x) can be decomposed to ξ = κ x + p, p ∈ ∂ J1 (x). κ = 1 in this paper. (2.8) is a nest iteration. The inner iteration takes k to be fixed, 0 ≤ k ≤ l, and iterative on j till jk . Here jk is defined as jk = min{j : ∥tk F ′ (uk )∗ (F (uk ) − z δ ) + (1 − tk )(uk − u0 )∥ ≤ τ δ} j
j
j
(2.9)
In (2.9), τ > ρ/δ is parameter in Morozov discrepancy principle, ρ is the radius of the neighborhood of u0 . By the inner j j iteration, we get ukk , and take it as the initial value of the (k + 1)th step, i.e. u0k+1 = ukk . The iteration will stop till the j ∗
stopping rule ∥F (uk ) − z δ ∥ ≤ τ δ is satisfied, and ukk∗ is obtained. The stopping step k∗ = min{k : ∥F (uk ) − z δ ∥ ≤ τ δ}. The choice of the sequence {tk } is flexible. We can guarantee k∗ is finite by selecting some tk . If tk∗ has not reached 1, but the stopping rule is satisfied, we stop the iteration. In the following, we prove that the error of the iterative sequence in Bregman distance is decreasing. j
Theorem 1.
j
Let u∗ ∈ Bρ (u0 ) is the solution of (2.1). Assume the operator F ′ (·) is bounded, that is,
∥F ′ (u)∥ ≤ 1,
u ∈ Bρ (u0 )
(2.10)
and F satisfies the nonlinear condition
∥F (u) − F (u˜ ) − F ′ (u˜ )(u − u˜ )∥ ≤ η∥F (u) − F (u˜ )∥ 0<η<
1 2
,
u, u˜ ∈ Bρ (u0 )
(2.11)
Moreover, we assume the following inequality holds
∥F (ujk ) − z − F ′ (ujk )(ujk − u∗ )∥ ≤ ηk ∥tk F ′ (ujk )∗ (F (ujk ) − z) + (1 − tk )(ujk − u0 )∥
(2.12)
where ηk satisfies 0 ≤ ηk ≤
τ − tk − (1 − tk ) ρδ τ tk + (tk )2
(2.13)
4
L. Li, C. Chen and B. Bi / Journal of Computational and Applied Mathematics 365 (2020) 112408
αkj is selected appropriately, and satisfies ⎧ 2 j ⎪ k ⎪ ⎪ ⎨ (t − 1) ρ + 1 − t η − tk (ηk tk +1) ≤ αk ≤ 1, 1 ≤ j ≤ j k k k τδ τ 2 ⎪ ⎪ ≤ αkj ≤ 1, j = 0 ⎪ ⎩ t (η t +1) ρ (tk+1 − 1) τ δ + 1 − tk+1 ηk − k+1 kτk+1
(2.14)
Therefore, we get the following conclusion J
D
j+1
j+1
ξk
J
j
(u∗ , uk ) − D j (u∗ , uk ) ≤ 0,
0 ≤ j ≤ jk ,
ξk
0 ≤ k ≤ k∗
(2.15)
Proof. First, we consider the case j ̸ = 0. Let
ξkj = ujk + pjk ,
j
j
pk ∈ ∂ J1 (uk )
(2.16)
Then
ξkj+1 − ξkj = ujk+1 − ujk + pjk+1 − pjk j+1
Integrating with pk
(2.17)
− pjk on both sides,
⟨ξkj+1 − ξkj , pjk+1 − pjk ⟩ = ⟨ujk+1 − ujk , pjk+1 − pjk ⟩ + ∥pjk+1 − pjk ∥2
(2.18)
By (2.8),
ξkj+1 − ξkj = −(αkj )−1 [tk F ′ (ujk )∗ F ′ (ujk ) + (1 − tk )I ]∗ [tk F ′ (ujk )∗ (F (ujk ) − z δ ) + (1 − tk )(ujk − u0 )] Then
⟨−(αkj )−1 [tk F ′ (ujk )∗ F ′ (ujk ) + (1 − tk )I ]∗ [tk F ′ (ujk )∗ (F (ujk ) − z δ ) + (1 − tk )(ujk − u0 )], pjk+1 − pjk ⟩ = ⟨ujk+1 − ujk , pjk+1 − pjk ⟩ + ∥pjk+1 − pjk ∥2 j+1
Since J1 (·) is convex, ⟨uk
− ujk , pjk+1 − pjk ⟩ ≥ 0,
∥pjk+1 − pjk ∥ ≤ −(αkj )−1 ∥[tk F ′ (ujk )∗ F ′ (ujk ) + (1 − tk )I ]∗ [tk F ′ (ujk )∗ (F (ujk ) − z δ ) + (1 − tk )(ujk − u0 )]∥
(2.19)
Let u, u˜ , uˆ ∈ D(F ),
ξ˜ ∈ ∂ J(u˜ ),
ξˆ ∈ ∂ J(uˆ )
By Bregman distance, J
J
J
ξ
ξ
ξ
D ˜ (u, u˜ ) − D ˆ (u, uˆ ) + D ˆ (u˜ , uˆ ) = ⟨ξ˜ − ξˆ , u˜ − u⟩
(2.20)
Taking
ξ˜ = ξkj+1 ,
u = u∗ ,
j+1
ξˆ = ξkj ,
u˜ = uk ,
j
uˆ = uk
into (2.20), J
D
j+1
j+1
ξk
J
j
j+1
(u∗ , uk ) − D j (u∗ , uk ) ≤ ⟨ξk ξk
− ξkj , ujk+1 − u∗ ⟩
= ⟨ξkj+1 − ξkj , ujk+1 − ujk ⟩ + ⟨ξkj+1 − ξkj , ujk − u∗ ⟩
(2.21)
Combined with (2.19),
⟨ξkj+1 − ξkj , ujk+1 − ujk ⟩ = ⟨ξkj+1 − ξkj , ξkj+1 − ξkj − (pjk+1 − pjk )⟩ = ∥ξkj+1 − ξkj ∥2 − ⟨ξkj+1 − ξkj , pkj+1 − pjk ⟩ ≤
2 (α
j 2 k)
∥ [tk F ′ (ujk )∗ F ′ (ujk ) + (1 − tk )I ]∗ [tk F ′ (ujk )∗ (F (ujk ) − z δ )
+ (1 − tk )(ujk − u0 )] ∥2
(2.22)
L. Li, C. Chen and B. Bi / Journal of Computational and Applied Mathematics 365 (2020) 112408
5
Substituting the second equation in (2.8) and (2.22) into (2.21), J
j+1
D
j+1
ξk
J
j
(u∗ , uk ) − D j (u∗ , uk ) ξk
2
∥[tk F ′ (uk )∗ F ′ (uk ) + (1 − tk )I ]∗ [tk F ′ (uk )∗ (F (uk ) − z δ ) + (1 − tk )(uk − u0 )]∥2 j (αk )2 1 − j ⟨[tk F ′ (ujk )∗ F ′ (ujk ) + (1 − tk )I ]∗ [tk F ′ (ujk )∗ (F (ujk ) − z δ ) + (1 − tk )(ujk − u0 )], ujk − u∗ ⟩ ≤
j
j
j
j
j
αk
2
=
∥tk F ′ (ujk )∗ (F (ujk ) − z δ ) + (1 − tk )(ujk − u0 )∥2
j 2 k)
(α 2
−
j
(αk )2
⟨tk F ′ (ujk )∗ (F (ujk ) − z δ ) + (1 − tk )(ujk − u0 ),
{I − [tk F ′ (ujk )∗ F ′ (ujk ) + (1 − tk )I ][tk F ′ (ujk )∗ F ′ (ujk ) + (1 − tk )I ]∗ } · [tk F ′ (ujk )∗ (F (ujk ) − z δ ) + (1 − tk )(ujk − u0 )]⟩ 1
−
∥tk F
αkj 1
′
δ
j j (uk )∗ (F (uk )
− z ) + (1 −
j j tk F ′ (uk )∗ (F (uk ) z δ ) j k j j j j tk F ′ (uk )∗ F (uk ) z F ′ (uk )(uk
+
α
[
≤[ +
2 j 2 k)
(α 1
1
αkj
]∥tk F ′ (ujk )∗ (F (ujk ) − z δ ) + (1 − tk )(ujk − u0 )∥2
⟨[tk F ′ (ujk )∗ (F (ujk ) − z δ ) + (1 − tk )(ujk − u0 )], tk F ′ (ujk )∗ [F (ujk ) − z − F ′ (ujk )(ujk − u∗ )]⟩
αkj
+[
− u∗ ) + z − z δ ] + (1 − tk )(u∗ − u0 )⟩
− −
−
− u0 )∥
(2.23) 2
+ (1 − tk )(ujk − u0 )],
−
⟨[
j tk )(uk
tk δ
(1 − tk )ρ
+
j k
α
α
j k
]∥tk F ′ (ujk )∗ (F (ujk ) − z δ ) + (1 − tk )(ujk − u0 )∥
Substituting (2.12) into (2.23), J
D
j+1
j+1
ξk
J
j
(u∗ , uk ) − D j (u∗ , uk ) ξk
≤ ∥tk F ′ (ujk )∗ (F (ujk ) − z δ ) + (1 − tk )(ujk − u0 )∥ · [(
2 j 2 k)
(α 1
1
−
αkj
)∥tk F ′ (uk )∗ (F (uk ) − z δ ) + (1 − tk )(uk − u0 )∥ j
j j t t F ′ (uk )∗ (F (uk ) j k k k k j j tk F ′ (uk )∗ (F (uk ) z δ )
+
η∥
α
≤∥ [(
j
2 j 2 k)
(α tk
+
j k
α
−
1 j k
α
+
α
− z δ + z δ − z) + (1 − tk )(ujk − u0 )∥ +
tk δ j k
α
+
(1 − tk )ρ
αkj
]
(2.24)
+ (1 − tk )(ujk − u0 )∥·
− t k ηk j k
j
)∥tk F ′ (uk )∗ (F (uk ) − z δ ) + (1 − tk )(uk − u0 )∥ j
(ηk tk + 1)δ + (1 − tk )
ρ αkj
j
j
]
By the stopping rule, ∥tk F ′ (uk )∗ (F (uk ) − z δ ) + (1 − tk )(uk − u0 )∥ > τ δ . Then (2.24) can be simplified to j
D
J
j+1
j+1
ξk
J
j
j
j
(u∗ , uk ) − D j (u∗ , uk ) ξk
≤ ∥tk F ′ (ujk )∗ (F (ujk ) − z δ ) + (1 − tk )(ujk − u0 )∥ · 2 1 tk ηk tk j j j j [( j − j + j + (ηk tk + 1))∥tk F ′ (uk )∗ (F (uk ) − z δ ) + (1 − tk )(uk )(uk − u0 )∥ (αk )2 αk αk τ αkj ρ + (1 − tk ) j ] αk
(2.25)
6
L. Li, C. Chen and B. Bi / Journal of Computational and Applied Mathematics 365 (2020) 112408 j
Here we take βk :=
2
1
−
j
(αk )2
+
j
αk
tk η k j
αk
tk
+
j
τ αk
(ηk tk + 1). Since
2
αkj ≥
(tk −
ρ 1) τ δ
+ 1 − tk ηk −
tk (ηk tk +1)
τ
we get
βkj ≤ (tk − 1)
ρ
≤0
j k
α τδ
Therefore,
βkj ∥tk F ′ (ujk )∗ (F (ujk ) − z δ ) + (1 − tk )(ujk − u0 )∥ + (1 − tk ) ≤ βkj τ δ + (1 − tk )
ρ
ρ αkj
αkj
≤0 At last, it is proved that J
D
j+1
j+1
ξk
J
j
(u∗ , uk ) − D j (u∗ , uk ) ≤ 0
(2.26)
ξk
In the same way, (2.26) holds when j = 0. Finally, we prove that the error of the iterative sequence in the Bregman distance is decreasing. □ j
J
Theorem 2. Let u0 ∈ D(F ) ∩ BV (Ω ), ξ0 ∈ L2 (Ω ). u∗ ∈ Bρ (u0 ) is the solution of F (u) = z, and uk ∈ Bρ (u0 ), D 0 (u∗ , u0k ) ≤ ξk
j∗ ukk∗
ρ2 2
.
j
We can get according to the iteration sequence {uk } in (2.8) when the stopping rule reaches k∗ = k∗ (δ ). If F satisfies (2.10)–(2.11), (2.12)–(2.14) hold, and
∥F (ujk ) − F (ukk∗ )∥ ≤ 2∥tk F ′ (ujk )∗ (F (ujk ) − z δ ) + (1 − tk )(ujk − u0 )∥ j ∗
then the number of the iteration
(2.27) jk∗
∑k∗
is limited. Taking δn > 0, as δn → 0, there exists a subsequence of {uk∗n } convergent
k=0 jk
n
jk∗
in weak∗ topology, and the limit u ∈ S(z). If S(z) ∩ Bρ (u∗ ) = {u∗ }, then uk∗n ⇀ u∗ , where k∗n is the stopping index. n
Proof. First, we prove that j ∗
J
D
j ∗
ξkk∗
∑k∗
k=0 jk
is finite. As k ≥ 0, Eq. (2.25) holds. Therefore,
J
(u∗ , ukk∗ ) − Dξ0 (u∗ , u0 ) ∗
≤
k j=jk −1 ∑ ∑ k=0
∥tk F ′ (ujk )∗ (F (ujk ) − z δ ) + (1 − tk )(ujk − u0 )∥·
j=0
[βkj ∥tk F ′ (ujk )∗ (F (ujk ) − z δ ) + (1 − tk )(ujk − u0 )∥ + (1 − tk )
(2.28)
ρ αkj
]
Since ∥tk F ′ (uk )∗ (F (uk ) − z δ ) + (1 − tk )(uk − u0 )∥ ≥ τ δ , it follows that j
j
j
∗
J
Dξ0 (u∗ , u0 ) ≥ D j
j ∗
J j
∗ k ∗ (u , uk∗ ) +
ξkk∗
k j=jk −1 ∑ ∑
−τ δ[βkj τ δ + (1 − tk )
j=0
k=0
j
ρ αkj j
Combine αk and βk , we make sure that there is a upper bound to βk τ δ + (1 − tk )
γk = − max{τ δ[βkj τ δ + (1 − tk ) j
ρ
(2.29)
] ρ j
αk
on j. Let
]} > 0
(2.30)
jk γk
(2.31)
αkj
Eq. (2.29) can be simplified as ∗
J
J
Dξ0 (u , u0 ) ≥ D ∗
(u , ∗
j ∗
ξkk∗
j ∗ ukk∗ )
+
k ∑ k=0
Obviously, the total iteration number
∑k∗
k=0 jk
is finite.
L. Li, C. Chen and B. Bi / Journal of Computational and Applied Mathematics 365 (2020) 112408
7
In the following, We prove the iterative sequence is convergent. j ∗
⟨ξk∗ , ukk∗ − u∗ ⟩ j ∗
= ⟨ξ0 , ukk∗ − u∗ ⟩ ∗
k jk −1 ∑ ∑ j − ⟨(αk )−1 [tk F ′ (ujk )∗ F ′ (ujk ) + (1 − tk )I ]∗ [tk F ′ (ujk )∗ (F (ujk ) − z δ ) k=0 j=0 j ∗
+ (1 − tk )(ujk − u0 )], ukk∗ − u∗ ⟩
(2.32)
∗
j ∗
≤ ∥ξ0 ∥ · ∥ukk∗ − u∗ ∥ −
k jk −1 ∑ ∑
(αk )−1 ∥tk F ′ (uk )∗ (F (uk ) − z δ ) + (1 − tk )(uk − u0 )∥2 j
j
j
j
k=0 j=0 k∗ jk −1
(αk )−1 ⟨tk F ′ (uk )∗ (F (uk ) − z δ ) + (1 − tk )(uk − u0 ),
∑∑
+
j
j
j
j
k=0 j=0
tk F (uk )∗ [F (uk ) − z δ − F ′ (uk )(ukk∗ − u∗ )] + (1 − tk )(uk − u0 − ukk∗ + u∗ )⟩ ′
j
j
j ∗
j
j
j ∗
For simplicity, let Lk := tk F ′ (uk )∗ (F (uk ) − z δ ) + (1 − tk )(uk − u0 ). By (2.12), the underlined part in Eq. (2.32) is j
j
j
⟨Lk , tk F ′ (ujk )∗ [F (ujk ) − z − F ′ (ujk )(ujk − u∗ )]⟩ + ⟨Lk , tk F ′ (ujk )∗ (z − z δ )⟩ j ∗
j ∗
+ ⟨Lk , (1 − tk )(ujk − u0 − ukk∗ + u∗ )⟩ + ⟨Lk , tk F ′ (ujk )∗ [−F ′ (ujk )(ukk∗ − ujk )]⟩ ≤ tk ηk ∥Lk ∥2 + tk δ (tk ηk + 1)∥Lk ∥ j ∗
j ∗
+ ⟨Lk , (1 − tk )(ujk − u0 − ukk∗ + u∗ )⟩ + ⟨Lk , tk F ′ (ujk )∗ [−F ′ (ujk )(ukk∗ − ujk )]⟩ By the nonlinear condition (2.11), j ∗
j ∗
∥F ′ (ujk )(ukk∗ − ujk )∥ ≤ (1 + ηk )∥F (ukk∗ ) − F (ujk )∥ j
Since uk ∈ Bρ (u0 ), j ∗
j ∗
∥ujk − u0 − ukk∗ + u∗ ∥ = ∥(ujk − u∗ ) + (u∗ − ukk∗ ) + (u∗ − u0 )∥ ≤ 5ρ then
⟨Lk , tk F ′ (ujk )∗ [F (ujk ) − z − F ′ (ujk )(ujk − u∗ )]⟩ + ⟨Lk , tk F ′ (ujk )∗ (z − z δ )⟩ j ∗
j ∗
+ ⟨Lk , (1 − tk )(ujk − u0 − ukk∗ + u∗ )⟩ + ⟨Lk , tk F ′ (ujk )∗ [−F ′ (ujk )(ukk∗ − u∗ )]⟩ ≤ tk ηk ∥Lk ∥ + tk δ (tk ηk + 1)∥Lk ∥ + 5ρ (1 − tk )∥Lk ∥ 2
j ∗
+ tk (1 + ηk )∥Lk ∥ · ∥F (ukk∗ ) − F (ujk )∥ Moreover, by (2.27) ,
⟨Lk , tk F ′ (ujk )∗ [F (ujk ) − z − F ′ (ujk )(ujk − u∗ )]⟩ + ⟨Lk , tk F ′ (ujk )∗ (z − z δ )⟩ j ∗
j ∗
+ ⟨Lk , (1 − tk )(ujk − u0 − ukk∗ + u∗ )⟩ + ⟨Lk , tk F ′ (ujk )∗ [−F ′ (ujk )(ukk∗ − u∗ )]⟩ ≤ tk ηk ∥Lk ∥2 + [tk δ (tk ηk + 1) + 5ρ (1 − tk )]∥Lk ∥ + 2tk (1 + ηk )∥Lk ∥2 Substituting it into Eq. (2.32), j ∗
⟨ξk∗ , ukk∗ − u∗ ⟩ ∗
≤ ∥ξ0 ∥ · ∥
j ∗ ukk∗
∗
−u ∥−
k jk −1 ∑ ∑
j
(αk )−1 [1 − 2tk − 3tk ηk ]∥Lk ∥2
k=0 j=0 k∗ jk −1
+
∑∑
j
(αk )−1 {[tk2 ηk + tk ]δ + 5ρ (1 − tk )}∥Lk ∥
k=0 j=0
By (2.26), we developed J
D
j ∗ ξkk∗
j ∗
J
(u∗ , ukk∗ ) ≤ Dξ 0 (u∗ , u0k ) ≤ k
ρ2 2
8
L. Li, C. Chen and B. Bi / Journal of Computational and Applied Mathematics 365 (2020) 112408
and we can estimate the following inequality by the Bregman distance 1
j ∗
J
D
j ∗
ξkk∗
(u∗ , ukk∗ ) ≥
2
j ∗
∥ukk∗ − u∗ ∥2
then j ∗
∥ukk∗ − u∗ ∥ ≤ ρ Therefore, j ∗
⟨ξk∗ , ukk∗ − u∗ ⟩ ∗
≤ ∥ξ0 ∥ρ −
k jk −1 ∑ ∑
j
(αk )−1 [1 − 2tk − 3tk ηk ]∥Lk ∥2
k=0 j=0 k∗ jk −1
+
∑∑
j
(αk )−1 {[tk2 ηk + tk ]δ + 5ρ (1 − tk )}∥Lk ∥
k=0 j=0
Moreover, Bregman distance D
j ∗
J j ∗
ξkk∗
(u∗ , ukk∗ ) ≥ 0, it follows that j ∗
j ∗
j ∗
J(ukk∗ ) ≤ J(u∗ ) − ⟨ξkk∗ , u∗ − ukk∗ ⟩ k∗
≤ J(u∗ ) + ∥ξ0 ∥ρ −
jk − 1 ∑∑
j
(αk )−1 [1 − 2tk − 3tk ηk ]∥Lk ∥2
k=0 j=0 ∗
+
k jk −1 ∑ ∑
j
(αk )−1 {[tk2 ηk + tk ]δ + 5ρ (1 − tk )}∥Lk ∥
(2.33)
k=0 j=0 jk∗
jk∗
n
n
We can prove ∥F (uk∗n ) − z δn ∥ → 0 as δn → 0, then ∥F (uk∗n ) − z ∥ → 0. Therefore, we can derive ∥Lk ∥ < ∞ and conclude j ∗
jk∗
that J(ukk∗ ) is uniformly bounded on small δ , and there is a subsequence of {uk∗n } is weak∗ topology convergent in BV(Ω ). n
jk∗
Moreover, F is continuous, weakly sequenced closed, then the limit u of the convergent subsequence of {uk∗n } belongs to S(z). If S(z)
⋂
jk∗
Bρ (u∗ ) = {u∗ }, then uk∗n ⇀ u∗ . n
n
□
3. Numerical analysis The principle of fluorescence molecule tomography(FMT) is that the fluorescence reporter group(such as GFP, RFP) is injected into the organism, and is excited by a specific wavelength of excited light. When the group absorbs the excited light, the energy level transition occurs to the state of high energy, and gradually attenuates back to the ground state, and the fluorescence with longer wavelength than the incident light is emitted. After excitation and fluorescence light interacts with tissue, they are captured by the surface detector. According to the data received by the probe, the mathematical model was built, and the optical parameters distribution of the tissue was reconstructed [27,28], and the high resolution image was obtained. The absorption coefficient and the scattering coefficient represent the probability that a photon is absorbed and scattered at the unit length, respectively. The unit is mm−1 or cm−1 . In this section, we will consider to reconstruct the absorption coefficient of the tissue. There are two kinds of light in fluorescent molecular tomography, one is an excitation light from the external light source(the wavelength is assumed as m), and the other is the fluorescent produced by the fluorescent markers that is injected into the target tissue under the radiation of excitation light (the wavelength is assumed as n). The chromophores and fluorophore in the tissue can absorb both kinds of light. We denote the absorption coefficients that the chromophores, the fluorophore and the tissue to the excitation light and fluorescent as µai,m , µai,n , µaf ,m , µaf ,n , µs,m , µs,n , respectively. Usually, the other optical parameters are known, and are used to reconstruct µai,m , µaf ,m . Now we consider the operator equation F (µai,m , µaf ,m ) = z δ z δ denotes the reading of the detector on the surface. Assume that there are d detectors, z δ ∈ Rd . F : D → Rd ,
(µai,m , µaf ,m ) ↦ → z δ
(3.1)
where µai,m , µaf ,m denote the absorption coefficient of fluorescence by chromoplast and fluorophore. D = {(µai,m , µaf ,m ) ∈ L∞ (X ) × L∞ (X )}.
L. Li, C. Chen and B. Bi / Journal of Computational and Applied Mathematics 365 (2020) 112408
9
Fig. 1. Light enters the columnar tissue.
The direct problem is to solve the reading of the detector on surface z δ by the absorption coefficients µai,m , µaf ,m of chromophore and fluorophore for fluorescence. On the contrary, the inverse problem is to solve µai,m , µaf ,m by z δ . Later, we will use Eq. (2.8) to solve our problem. In Sections 3.2 and 3.3, we will give the models for solving the direct problem and inverse problem of FMT, respectively. 3.1. Radiation Transfer Equation for FMT Radiation Transfer Equation (RTE) can effectively simulate the propagation of photons in tissues and is widely used in biomedical imaging. Radiation transfer equation is used to simulate the scattering, reflection and absorption of light in tissue, and uses radiation rate and phase function p(sˆ, sˆ′ ) to describe the intensity of light and the probability distribution of photon scattering, respectively. Radiation rate φ (r , sˆ) represents the flux density of average power in a given direction sˆ and in unit solid angle, and the unit is W m−2 sr−1 Hz−1 . The phase function p(sˆ, sˆ′ ) represents the probability of being incident on a scattering particle from the direction sˆ′ , and then being scattered into the direction sˆ. p(sˆ, sˆ′ )dΩ ′ describes the probability that a photon incident from the direction sˆ′ is scattered into some solid angle in the direction sˆ,
∫ 4π
p(sˆ, sˆ′ )dΩ ′ = 1
It is generally assumed that the phase function is only related to the angle between the incident and exit directions, p(sˆ, sˆ′ ) = p(sˆ′ , sˆ) = p(cosθ ) The Henyey–Greenstein function is chosen as the phase function of the model, i.e., p(cos θ ) =
1 − g2
1
4π (1 + g 2 − 2gcosθ ) 32
where g is an anisotropic coefficient, and is defined as
∫ g = ⟨cosθ⟩ =
4π
p(sˆ, sˆ′ )(sˆ · sˆ′ )dΩ ′
∫
4π
p(sˆ, sˆ′ )dΩ ′
∫ = 4π
p(sˆ, sˆ′ )(sˆ, sˆ′ )dΩ ′
Examine a columnar structure with a length of ds and a unit cross section as shown in Fig. 1. Light having an emissivity of φ (r , sˆ) is incident on the cylindrical voxel. The reason for the change in emissivity can be reduced to the following aspects: (1) The emissivity will disappear due to absorption and scattering by the tissue. The reduction in emissivity can be expressed as: dφ (r , sˆ) = −(µa + µs )dsφ (r , sˆ) (2) The energy obtained by external scattering. This part of the increment can be expressed as dφ (r , sˆ) = µs ds
∫ 4π
p(sˆ, sˆ′ )φ (r , sˆ′ )dsˆ′
(3) If the light source emits photon in the body, the emissivity will also be increased. Let the power radiated from the internal light source along the direction sˆ in unit angle to be q(r , sˆ). The increment is dφ (r , sˆ) = q(r , sˆ)ds
10
L. Li, C. Chen and B. Bi / Journal of Computational and Applied Mathematics 365 (2020) 112408
Fig. 2. Angular dispersion diagram.
These three factors are considered, the radiation transfer equation can be derived by the law of conservation of energy as dφ (r , sˆ) ds
= −(µa + µs )φ (r , sˆ) + µs
∫ 4π
p(sˆ, sˆ′ )φ (r , sˆ′ )dsˆ′ + q(r , sˆ)
which is also written as Hamilton operator equation: (s · ∇ + µa (r) + µs (r))φ (r , sˆ) = µs (r)
∫ 4π
p(sˆ, sˆ′ )φ (r , sˆ′ )dsˆ′ + q(r , sˆ)
RTE can reflect the energy balance relation of every points in the media. Fluorescence Molecule Tomography (FMT) includes two radiation transfer equations, which describe propagation processes that the fluorescence produced from exciting light and fluorescent groups in the tissue. The propagation equation of exciting light is (s · ∇ + µai,m (r) + µaf ,m (r) + µs,m (r))φm (r , sˆ) = µs,m (r)
∫ 4π
p(sˆ, sˆ′ )φm (r , sˆ′ )dsˆ′
The propagation equation of fluorescent is (s · ∇ + µai,n (r) + µaf ,n (r) + µs,n (r))φn (r , sˆ)
= µs,n (r)
∫ 4π
p(sˆ, sˆ′ )φn (r , sˆ′ )dsˆ′ + ηµaf ,m (r)
∫ 4π
φm (r , sˆ′ )dsˆ′
The boundary condition is
φm (r , sˆ) = gm (r , sˆ),
φn (r , sˆ) = 0 on Γ−
Here φm (r , sˆ′ ) and φn (r , sˆ′ ) denote the emissivity of exciting light and fluorescent, and η is determined by the property of fluorescent group in the tissue. Γ− is incident boundary. 3.2. Direct model of FMT The direct problem of FMT is to solve the detector readings z δ on the body by using the absorption coefficients µai,m and µaf ,m of fluorescence by chromoplast and fluorophore. Aiming at 2D system, we use the finite element method to discrete the space to get the matrix equation, and then to solve this equation by SOR iteration method. At last, we construct the complete direct model. First of all, the data is spatially discrete and the region is triangulated. The region Ω is decomposed into N triangular areas, which are discretized into I × J points r = (xi , yj ), i = 0, . . . , I, j = 0, . . . , J. Then the angle sˆ is discretized, and [0, 2π ) is evenly divided into K directions. The length of each cell is △θ , which refers to Fig. 2. Here sˆ = θk , k = 1, . . . , K . The radiation rate can be approximated in sections:
φ (xi , yj , θ ) ≈
K ∑ k=1
φ k,i,j Lk (θ )
L. Li, C. Chen and B. Bi / Journal of Computational and Applied Mathematics 365 (2020) 112408
11
φ k,i,j , k = 1, . . . , K is the radiation rate along the kth direction at (xi , yj ), Lk (θ ) is piecewise linear function of θ , { 1, k = k′ Lk (θk′ ) = 0, others After dispersion, the RTE equation of fluorescence molecular tomography can be reduced to the following equations:
∂φmk,i,j ∂φmk,i,j + sin θk + (µai,m + µaf ,m + µs,m )φmk,i,j ∂x ∂y ∫ k,i,j ′ (θ )dθ ′ = µs,m p(θk , θ ′ )φm
cos θk
2π
∂φnk,i,j ∂φnk,i,j cos θk + sin θk + (µai,n + µaf ,n + µs,n )φnk,i,j ∂x ∂y ∫ ∫ = µs,n p(θk , θ ′ )φnk,i,j (θ ′ )dθ ′ + ηµaf ,m φmk,i,j (θ ′ )dθ ′ 2π
φmk,i,j = gmk,i,j ,
2π
φnk,i,j = 0 on Γ−
The integral on the right-hand side of the equation is expressed as
∫ 2π
p(θk , θ ′ )φ k,i,j (θ ′ )dθ ′ =
K ∑
αk′ pk,k′ φ k ,i,j and ′
k′ =1
∫ 2π
φ k,i,j (θ ′ )dθ ′ =
K ∑
βk′ φ k ,i,j ′
k′ =1
αk′ and βk′ are selected weight. Then we use the upwind difference scheme to approximate the partial derivative of the radiation rate φ at (xi , yj ). Then 2D RTE is obtained after the angle is discretized: φmk,i,j − φmk,i−1,j φmk,i,j − φmk,i,j−1 + sin θk + (µai,m + µaf ,m + µs,m )φmk,i,j △x △y K ∑ ′ = µs,m αk′ pk,k′ φmk ,i,j
cos θk
k′ =1 k,i,j n
φnk,i,j − φnk,i,j−1 − φnk,i−1,j + sin θk + (µai,n + µaf ,n + µs,n )φnk,i,j △x △y K K ∑ ∑ ′ ′ = µs,n αk′ pk,k′ φnk ,i,j + ηµaf ,m βk′ φmk ,i,j , cos θk
φ
k′ =1
k′ =1
φmk,i,j = gmk,i,j ,
φnk,i,j = 0 on Γ−
(3.2)
Obviously, Eq. (3.2) can be written as matrix equation: Aφ = b. Use SOR method, the following iteration algorithm is obtained:
[φmk,i,j ]l+1 = ρ
µs,m
∑K
k′ =1
αk′ pk,k′ [φmk ,i,j ]l + [φmk,i−1,j ]l cos θk /△x + [φmk,i,j−1 ]l sin θk /△y cos θk /△x + sin θk /△y + (µai,m + µaf ,m + µs,m ) ′
+ (1 − ρ )[φmk,i,j ]l [φnk,i,j ]l+1 =ρ
µs,n
(3.3)
∑K
k′ =1
αk′ pk,k′ [φ
k′ ,i,j l n
∑K
k′ ,i,j l m
k,i−1,j l n
k,i,j−1 l n
] + ηµaf ,m k′ =1 βk′ [φ ] + [φ ] cos θk /△x + [φ cos θk /△x + sin θk /△y + (µai,n + µaf ,n + µs,n )
] sin θk /△y
+ (1 − ρ )[φnk,i,j ]l where ρ is super relaxation coefficient. The detector reading at the boundary lattice point can be calculated according to the following formula
Φ (xi , yj ) =
∫
φ (xi , yj , sˆ)dsˆ ≈ AP
k2 ∑
ak φk,i,j = Φij
k=k1
The weight ak can be obtained by generalized Trapezoidal rule. Therefore, we will get the corresponding detector readings F (µai,m , µaf ,m ) when the optical parameters µai,m and µaf ,m are given arbitrary.
12
L. Li, C. Chen and B. Bi / Journal of Computational and Applied Mathematics 365 (2020) 112408
Fig. 3. Subdivision graph of tissue.
Fig. 4. True absorption coefficient and scattering coefficient distribution.
3.3. Inverse model of FMT The inverse problem of FMT is to use the reading z δ of known detectors to inverse the absorption coefficients µai,m and µaf ,m of chromosome and fluorophore to fluorescence. Set u = (µai,m , µaf ,m ). The nonlinear equation of FMT imaging can be written as F (u) = z δ . There are d detectors, z δ ∈ Rd . Applying the iteration algorithm (2.8) to inverse the absorption coefficients µai,m and the scattering coefficients µaf ,m of the chromophores to the excited light. 3.4. Numerical results Consider the 2D tissue (Fig. 3). The graph is discretized into 276 triangulars, 32 directional angles and 156 vertices. Place 16 excitation light sources and 16 receivers at regular intervals along the region boundary. Two abnormal bodies are placed in the tissue (refer to Fig. 4(a) and (b)). The absorption coefficients of the anomalous bodies are 0.4, and the absorption coefficient of background of the tissue is 0.2. The scattering coefficients of the anomalous bodies are 4, and the scattering coefficient of the background of the tissue is 2. We only focus on the inversion of absorption coefficient of one abnormal body (Fig. 4(a)).
L. Li, C. Chen and B. Bi / Journal of Computational and Applied Mathematics 365 (2020) 112408
13
Fig. 5. (a)–(d). The reconstructions of four methods for the absorption coefficient of the tissue.
Table 1 Comparison of four regularization methods.
Iteration steps Running time (s) Relative error
Tikhonov regularizatioin
l1 regularization
Total variational regularization
New regularization
24 578.34 0.2864
9 209.4288 0.2378
15 349.79 0.2541
12 330.62 0.2338
The simulation results of the absorption coefficient are obtained by using Tikhonov regularization, regularization with l1 norm, total variational regularization, and the new regularization method proposed in this paper, respectively. As shown in Fig. 5(d), the proposed method can exactly inverse the position of abnormal body, and depict clearly the boundary information of abnormal body. Compared with Tikhonov regularization and total variational regularization, the proposed new method runs faster, and the inversion of the abnormal body is closer to the true value. Furthermore, the inversion effect of the proposed new method on the background is good. From Table 1 and Fig. 5, l1 regularization method is the fastest, and the inversion effect on the abnormal body and background are better in these four methods. However, the proposed new method has the smallest relative error. Tikhonov regularization method is relatively worse in both computation time and inversion effect.
14
L. Li, C. Chen and B. Bi / Journal of Computational and Applied Mathematics 365 (2020) 112408
4. Conclusion Total variational regularization can effectively inverse the piecewise constant function, and can preserve the boundary of large object, thus are widely used to the field of image processing. Moreover, homotopy method can enlarge the range of iterative initial value. Therefore, a new total variational regularization method is proposed based on total variational regularization and fix-point homotopy. Moreover, the convergence of the new method is proved in theory. In the stage of numerical simulation, the absorption coefficient inversion problem of fluorescence molecular tomography (FMT) was solved, and the radiation transfer equation (RTE) was used to simulate the propagation of photons in tissues. Experimental results show that the new algorithm can accurately inverse the position of anomalous body and clearly depict the edge information of anomalous body. Although four regularization methods can calculate the position of abnormal body, Tikhonov regularization method and total variational regularization are not very good at restoring the image background. In contrast, l1 method and the new approach better preserve the background information, their whole simulation results are closer to the true value, and their accuracy are also slightly higher. However, the relative error of the proposed new method is smallest. Therefore, the proposed new method is a stable and effective regularization method. References [1] V. Ntziachristos, C. Bremer, E.E. Graves, et al., In vivo tomographic imaging of near-infrared fluorescent probes, Mol. Imaging 1 (2) (2002) 82–88. [2] V. Ntziachristos, C.H. Tung, C. Bremer, et al., Fluorescence molecular tomography resolves protease activity in vivo, Nature Med. 8 (7) (2002) 757–760. [3] S. Chandrasekhar, Radiative Transfer, Clarendon, London, 1950. [4] K.M. Case, P.F. Zweifel, Linear Transport Theory, Addison-Wesley, Reading, MA, 1967. [5] V. Agoshkov, Boundary Value Problems for Transport Equations, Birkhauser, Boston, MA, 1998. [6] G. Choppin, J.-O. Liljenzin, J. Rydberg, Radiochemistry and Nuclear Chemistry, third ed., Butterworth-Heinemann, Woburn, MA, 2002. [7] W.X. Cong, L.V. Wang, G. Wang, Formulation of photon diffusion from spherical bioluminescent sources in an infinite homogeneous medium, Biomed. Eng. Online 3 (1) (2004) 1–6. [8] C.H. Contag, B.D. Ross, It’s not just about anatomy: in vivo bioluminescence imaging as an eyepiece into biology, J. Mag. Reson. Imaging 16 (4) (2002) 378–387. [9] A.D. Klose, U. Netz, J. Beuthan, et al., Optical tomography using the time-independent equation of radiative transfer−Part 1: forward model, J. Quant. Spectrosc. Radiat. Transfer 72 (5) (2002) 691–713. [10] G. Hao, K.Z. Hong, A fast forward solver of radiative transfer equation, Transp. Theory Stat. Phys. 38 (3) (2009) 149–192. [11] B.C. Wilson, G. Adam, A Monte Carlo model for the absorption and flux distributions of light in tissue, Med. Phys. 10 (6) (1983) 824. [12] L. Wilson, S.L. Jacques, L. Zheng, MCML-Monte Carlo Modeling of light transport in multi-layered tissues, Comput. Methods Programs Biomed. 47 (2) (1995) 131–146. [13] D. Boas, J. Culver, J. Stott, et al., Three dimensional Monte Carlo code for photon migration through complex heterogeneous media including the adult human head, Opt. Express 10 (3) (2002) 159–170. [14] A. Kienle, M.S. Patterson, Improved solutions of the steady-state and the time-resolved diffusion equations for equations for reflectance from a semi-infinite turbid medium, J. Opt. Soc. Amer. A 14 (1) (1997) 246–254. [15] H. Egger, M. Schlottbom, Analysis and regularization of problems in diffuse optical tomography, SIAM J. Math. Anal. 42 (2010) 1934–1948. [16] A.D. Klose, A.H. Hielscher, Iterative reconstruction scheme for optical tomography based on the equation of radiative transfer, Med. Phys. 26 (8) (1999) 1698–1707. [17] K. Ren, T. Ding, S. Vallelian, A one-step reconstruction algorithm for quantitative photoaccoustic imaging, Inverse Problems 31 (2015) 095005. [18] A.N. Tikhonov, V.Y. Arsenin, Solutions of Ill-Posed Problems, Wiley, Washington, DC, 1977. [19] M. Hanke, A. Neubauer, O. Scherzer, A convergence analysis of the Landweber iteration for nonlinear ill-posed problems, Numer. Math. 72 (1995) 21–37. [20] Parimah Kazemi, Robert J. Renka, A Levenberg–Marquardt method based on Sobolev gradients, Nonlinear Anal. TMA 75 (16) (2012) 6170–6179. [21] H. Gao, H.K. Zhao, Multilevel bioluminescence tomography based on radiative transfer equation, Part 2: total variation and l1 data fidelity, Opt. Express 18 (3) (2010) 2894. [22] D. Hyde, E.L. Miller, D.H. Brooks, et al., Data specific spatially varying regularization for multimodal fluorescence molecular tomography, IEEE Trans. Med. Imaging 29 (2) (2010) 365–374. [23] J. Dutta, S. Ahn, C. Li, et al., Joint L1 and total variation regularization for fluorescence molecular tomography, Phys. Med. Biol. 57 (6) (2012) 1459–1476. [24] Joyita Dutta, Snagtae Ahn, Changqing Li, Simon R. Cherry, Richard M. Leahy, Joint L1 and total variation regularization for fluorescnece molecular tomography, Phys. Med. Biol. 57 (6) (2012) 1459–1476. [25] B. Bi, B. Han, W. Han, J. Tang, L. Li, Image reconstruction for diffuse optical tomography based on radiative transfer equation, Comput. Math. Methods Med. 2015 (2015). [26] B. Bi, J. Tang, Diffuse optical tomography based on radiative transfer equation with L∞ data fidelity and total variation penalty regularization, J. Med. Imaging Health Inf. 7 (6) (2017) 1212–1218. [27] I. Georgakoudi, B.C. Jacobson, D.J. Van, et al., Fluorescence, reflectance, and light-scattering spectroscopy for evaluating dysplasia in patients with Barrett’s esophagus, Gastroenterology 120 (7) (2001) 1620–1629. [28] Vasilis Ntziachristos, Jorge Ripoll, Lihong V. Wang, Ralph Weissleder, Looking and listening to light: the evolution of whole-body photonic imaging, Nature Biotechnol. 23 (2005) 313–320.
Contents lists available at ScienceDirect
Journal of Computational and Applied Mathematics journal homepage: www.elsevier.com/locate/cam
A new total variational regularization method for nonlinear inverse problems in fluorescence molecular tomography✩ ∗
Li Li a , , Chen Chen a , Bo Bi b a b
School of Mathematics, Harbin Institute of Technology, Harbin, 150001, PR China School of Mathematics and Statics, Northeast Petroleum University, Daqing, 163318, PR China
article
info
Article history: Received 15 November 2018 Received in revised form 31 May 2019 Keywords: Ill-posed problems Total variational regularization Bregman distance Homotopy
a b s t r a c t Fluorescence molecular tomography (FMT) is an imaging way of in vivo optical imaging, which is widely used in biomedical photon imaging for its higher sensitivity. For years, the propagation process of photons in tissue is simulated based on the equation of radiative transfer. In this work, we analyze an inverse problem in FMT, which uses the reading of body surface detector to reconstruct the optical parameter distribution. We propose a novel total variational regularization method to solve the absorption parameter inversion problem in FMT. The proof of the proposed method is also provided. Finally, numerical experiments demonstrate the feasibility of the proposed method. © 2019 Elsevier B.V. All rights reserved.
1. Introduction Since the classical medical tomography depends highly on the physical properties and physiological change of an organization, the lesion of tissue cannot be discovered timely, which delays the treatment of the patient. Biomedical photon imaging is a non-invasive imaging technique, which can overcome the weakness of traditional medical imaging technology and satisfy the demand of medical diagnosis for comfort and nondestructive nature through detecting minute changes in tissue. Diffusion optical tomography (DOT), Fluorescence molecule tomography (FMT), and Bioluminescence tomography (BLT) are some typical biomedical photon imaging technology. In these imaging techniques, FMT can be traced back to 1980s, and is important for its simplicity, low cost, and ability to monitor an organization for a long time. Moreover, this kind of imaging result can satisfy our need for quantitative analysis. FMT was formally produced in 1990s, which obtains the high resolution images by inversing the optical parameters of biological tissue, such as absorption coefficient and reflection coefficient [1]. Moreover, the feasibility of this technology has been proved in the experiments on animals [2]. FMT is used to reconstruct the distribution of fluorescence concentration according to the fluorescence characteristics of near-infrared fluorescent chromosomes. Normally the radiative transport equation (RTE), which is used to describe the propagation or radiation of particles inside a medium, can be used in reconstruction of FMT as a forward model. If the function spaces are chosen appropriately, the solvability of the RTE and the properties of the solutions are all intensively studied [3–6]. Analytical method and numerical method are two main methods for solving RTE equation. The analytical method is considered to solve the problems when the shape of the tissue is regular. Cong considered the tissue as an infinite medium, and the light source as a regular geometry, then proposed an algorithm to obtain the analytic ✩ This work is supported by a scholarship from the China Scholarship Council (No. 201806125006). ∗ Corresponding author. E-mail address: [email protected] (L. Li). https://doi.org/10.1016/j.cam.2019.112408 0377-0427/© 2019 Elsevier B.V. All rights reserved.
2
L. Li, C. Chen and B. Bi / Journal of Computational and Applied Mathematics 365 (2020) 112408
solution [7]. However, the numerical method is more preferred to the situation that the medium is irregular, such as FEM [8], SOR iteration [9], Gauss Seidel iteration, and random statistical method [10–14], et al. Since RTE-based inverse problems are difficult either in theoretical analysis or in numerical computation, more investigations prefer to the diffusion approximation of the RTE (DA) [15], which is first-order phase approximation of RTE and is used to simulate the light migration in photoacoustic imaging. Based on gradient, Klose et al. developed an iterative algorithm to reconstruct distribution of tissue parameters by using the gradient information of operator F about the optical parameter [16]. Similarly, Kui Ren used the following relation
∂p |(0,T )×∂ Ω = Λ(c , σ , g) ∂n to reconstruct (c , σ ) in one-step reconstruction algorithm [17]. Compared with the traditional Newton type algorithm, the method based on gradient does not need to calculate the Jacobi matrix many times, and avoids the repeated matrix inversion, which greatly simplifies the calculation process [9]. Generally speaking, the inverse problems of FMT are ill-posed due to its high degree of absorption and scattering of light through tissue. Effective regularization methods should be utilized to deal with the ill-posedness. Tikhonov regularization [18], Landweber iteration [19] and Levenberg–Marquardt iteration [20], et al. are well-known widespread regularization for solving ill-posed problems. Here, L2 norm is taken as penalty term, and is feasible for solving the problem of smoothness. However, properties of some non-smooth edges are neglected. Compared with L2 norm, total variational norm preserves local smoothness and piecewise constant in the reconstructed images. Therefore, the regularization with the penalty of total variational norm can achieve better reconstructions compared to some classical iteration methods with the penalty term of L2 norm. Gao and Zhao investigated the regularization with both l1 and total-variational norm for bioluminescence tomography based on radiative transfer equation [21]. In 2010, Hyde et al. proposed to solve FMT problems by using two step regularization methods [22]. Thereafter, Dutta et al. introduced the total variational regularization method with l1 penalty, and put forward a new algorithm to deal with the problem of sparse images reconstructed by FMT in the local smooth high-density region [23], which combined with sparsity constrained method and smoothness condition. Joyita et al. proposed an approach involving combination of L1 and total variation norm penalties, which could suppress spurious background signals, enforce sparsity, and preserve local smoothness and piecewise constancy in the reconstructed images [24]. In 2015, B. Bi et al. investigated the feasibility of the sparsity reconstruction approach for DOT based on RTE [25]. Two years later, combined with total variation penalty, B. Bi and J. Tang studied the L∞ norm fitting method for DOT based on RTE [26]. Based on the properties of preserving shapes of total variational regularization and the wide convergent range of homotopy method, this article constructs a new iterative algorithm for solving this problem. The rest of the paper is organized as follows. Firstly, we prove that the error of the new method decrease monotonically in Bregman distance. Secondly, the new method with the noisy data is proved to be a regularization method theoretically. Finally, in the section of numerical simulation, the propagation of photons in tissues is simulated by using RTE, and the new iteration is applied to inverse the absorption parameters in FMT . 2. New total variational regularization and its convergence analysis Considering nonlinear operator equation: F (u) = z
(2.1)
where F is a nonlinear operator, F : D ⊂ L (Ω ) → Z . F is continuous, weakly sequenced closed, and Fréchet differential. u and z are the identified parameter and the observed data, respectively. Usually, we cannot obtain the exact data z, but the noisy data z δ satisfies ∥z δ − z ∥ ≤ δ , where δ is the error level. In fact, the problem (2.1) with noisy data is 2
F (u) = z δ
(2.2)
In order to solve (2.2), we consider the following least-square functional: min J(u) = min
u∈D (F )
u∈D (F )
1 2
∥F (u) − z δ ∥2
(2.3)
Its corresponding Euler equation is F ′ (u)∗ (F (u) − z δ ) = 0
(2.4)
Introduce the parameter t, we construct the homotopy: H(u, t) = t [F ′ (u)∗ (F (u) − z δ )] + (1 − t)(u − u0 ) = 0
(2.5)
u0 is an initial value. As t = 0, the solution of Eq. (2.5) is u = u0 . As t = 1, the solution of Eq. (2.5) is the solution of Euler equation (2.4). If the initial value is known, we can find the final solution of Eq. (2.2) along the homotopy curve.
L. Li, C. Chen and B. Bi / Journal of Computational and Applied Mathematics 365 (2020) 112408
3
[0, 1] is divided into 0 = t0 < t1 < · · · < tl = 1. t = tk (0 ≤ k ≤ l) is considered. Assume the iteration has reached ujk , j
and the higher-order term is ignored. We linearize F (u) at the point uk : j
j
j
F (u) ≈ F (uk ) + F ′ (uk )(u − uk )
(2.6) j+1
Substituting (2.6) into (2.5), and taking the computed u as uk , j+1
uk
= ujk − [tk F ′ (ujk )∗ F ′ (ujk ) + (1 − tk )I ]∗ [tk F ′ (ujk )∗ (F (ujk ) − z δ ) + (1 − tk )(ujk − u0 )]
(2.7)
In the following, the total variational regularization is introduced. Based on Bregman distance, a new iterative algorithm for solving nonlinear ill-posed problems is constructed:
⎧ j+1 j j j j ⎪ uk = arg minu∈D(F ) {⟨[tk F ′ (uk )∗ F ′ (uk ) + (1 − tk )I ]∗ [tk F ′ (uk )∗ (F (uk ) − z δ ) ⎪ ⎪ ⎪ ⎪ ⎨ + (1 − tk )(ujk − u0 )], u − ujk ⟩ + αkj DJ j (u, ujk )} ξk
(2.8)
⎪ ⎪ ξkj+1 = ξkj − (αkj )−1 [tk F ′ (ujk )∗ F ′ (ujk ) + (1 − tk )I ]∗ [tk F ′ (ujk )∗ (F (ujk ) − z δ ) ⎪ ⎪ ⎪ ⎩ + (1 − tk )(ujk − u0 )] J
j
j
Here D j (u, uk ) is defined as the Bregman distance between u and uk : ξk J
j
j
j
j
D j (u, uk ) = J(u) − J(uk ) − ⟨ξk , u − uk ⟩ ξk j
j
J
j
J
j
j
where ξk ∈ ∂ J(uk ) is subgradient. Obviously, D j (u, uk ) ≥ 0, D j (u, uk ) = 0 if and only if u = uk . J(·) is defined as ξk
ξk
J(·) = J1 (·) + J2 (·) where J1 (·) =| · |BV (Ω ) +χD(F ) (·),
J 2 (· ) =
κ 2
∥ · ∥2L2 (Ω )
Suppose D(F ) is convex, then J(·) and J1 (·) are convex. Any ξ ∈ ∂ J(x) can be decomposed to ξ = κ x + p, p ∈ ∂ J1 (x). κ = 1 in this paper. (2.8) is a nest iteration. The inner iteration takes k to be fixed, 0 ≤ k ≤ l, and iterative on j till jk . Here jk is defined as jk = min{j : ∥tk F ′ (uk )∗ (F (uk ) − z δ ) + (1 − tk )(uk − u0 )∥ ≤ τ δ} j
j
j
(2.9)
In (2.9), τ > ρ/δ is parameter in Morozov discrepancy principle, ρ is the radius of the neighborhood of u0 . By the inner j j iteration, we get ukk , and take it as the initial value of the (k + 1)th step, i.e. u0k+1 = ukk . The iteration will stop till the j ∗
stopping rule ∥F (uk ) − z δ ∥ ≤ τ δ is satisfied, and ukk∗ is obtained. The stopping step k∗ = min{k : ∥F (uk ) − z δ ∥ ≤ τ δ}. The choice of the sequence {tk } is flexible. We can guarantee k∗ is finite by selecting some tk . If tk∗ has not reached 1, but the stopping rule is satisfied, we stop the iteration. In the following, we prove that the error of the iterative sequence in Bregman distance is decreasing. j
Theorem 1.
j
Let u∗ ∈ Bρ (u0 ) is the solution of (2.1). Assume the operator F ′ (·) is bounded, that is,
∥F ′ (u)∥ ≤ 1,
u ∈ Bρ (u0 )
(2.10)
and F satisfies the nonlinear condition
∥F (u) − F (u˜ ) − F ′ (u˜ )(u − u˜ )∥ ≤ η∥F (u) − F (u˜ )∥ 0<η<
1 2
,
u, u˜ ∈ Bρ (u0 )
(2.11)
Moreover, we assume the following inequality holds
∥F (ujk ) − z − F ′ (ujk )(ujk − u∗ )∥ ≤ ηk ∥tk F ′ (ujk )∗ (F (ujk ) − z) + (1 − tk )(ujk − u0 )∥
(2.12)
where ηk satisfies 0 ≤ ηk ≤
τ − tk − (1 − tk ) ρδ τ tk + (tk )2
(2.13)
4
L. Li, C. Chen and B. Bi / Journal of Computational and Applied Mathematics 365 (2020) 112408
αkj is selected appropriately, and satisfies ⎧ 2 j ⎪ k ⎪ ⎪ ⎨ (t − 1) ρ + 1 − t η − tk (ηk tk +1) ≤ αk ≤ 1, 1 ≤ j ≤ j k k k τδ τ 2 ⎪ ⎪ ≤ αkj ≤ 1, j = 0 ⎪ ⎩ t (η t +1) ρ (tk+1 − 1) τ δ + 1 − tk+1 ηk − k+1 kτk+1
(2.14)
Therefore, we get the following conclusion J
D
j+1
j+1
ξk
J
j
(u∗ , uk ) − D j (u∗ , uk ) ≤ 0,
0 ≤ j ≤ jk ,
ξk
0 ≤ k ≤ k∗
(2.15)
Proof. First, we consider the case j ̸ = 0. Let
ξkj = ujk + pjk ,
j
j
pk ∈ ∂ J1 (uk )
(2.16)
Then
ξkj+1 − ξkj = ujk+1 − ujk + pjk+1 − pjk j+1
Integrating with pk
(2.17)
− pjk on both sides,
⟨ξkj+1 − ξkj , pjk+1 − pjk ⟩ = ⟨ujk+1 − ujk , pjk+1 − pjk ⟩ + ∥pjk+1 − pjk ∥2
(2.18)
By (2.8),
ξkj+1 − ξkj = −(αkj )−1 [tk F ′ (ujk )∗ F ′ (ujk ) + (1 − tk )I ]∗ [tk F ′ (ujk )∗ (F (ujk ) − z δ ) + (1 − tk )(ujk − u0 )] Then
⟨−(αkj )−1 [tk F ′ (ujk )∗ F ′ (ujk ) + (1 − tk )I ]∗ [tk F ′ (ujk )∗ (F (ujk ) − z δ ) + (1 − tk )(ujk − u0 )], pjk+1 − pjk ⟩ = ⟨ujk+1 − ujk , pjk+1 − pjk ⟩ + ∥pjk+1 − pjk ∥2 j+1
Since J1 (·) is convex, ⟨uk
− ujk , pjk+1 − pjk ⟩ ≥ 0,
∥pjk+1 − pjk ∥ ≤ −(αkj )−1 ∥[tk F ′ (ujk )∗ F ′ (ujk ) + (1 − tk )I ]∗ [tk F ′ (ujk )∗ (F (ujk ) − z δ ) + (1 − tk )(ujk − u0 )]∥
(2.19)
Let u, u˜ , uˆ ∈ D(F ),
ξ˜ ∈ ∂ J(u˜ ),
ξˆ ∈ ∂ J(uˆ )
By Bregman distance, J
J
J
ξ
ξ
ξ
D ˜ (u, u˜ ) − D ˆ (u, uˆ ) + D ˆ (u˜ , uˆ ) = ⟨ξ˜ − ξˆ , u˜ − u⟩
(2.20)
Taking
ξ˜ = ξkj+1 ,
u = u∗ ,
j+1
ξˆ = ξkj ,
u˜ = uk ,
j
uˆ = uk
into (2.20), J
D
j+1
j+1
ξk
J
j
j+1
(u∗ , uk ) − D j (u∗ , uk ) ≤ ⟨ξk ξk
− ξkj , ujk+1 − u∗ ⟩
= ⟨ξkj+1 − ξkj , ujk+1 − ujk ⟩ + ⟨ξkj+1 − ξkj , ujk − u∗ ⟩
(2.21)
Combined with (2.19),
⟨ξkj+1 − ξkj , ujk+1 − ujk ⟩ = ⟨ξkj+1 − ξkj , ξkj+1 − ξkj − (pjk+1 − pjk )⟩ = ∥ξkj+1 − ξkj ∥2 − ⟨ξkj+1 − ξkj , pkj+1 − pjk ⟩ ≤
2 (α
j 2 k)
∥ [tk F ′ (ujk )∗ F ′ (ujk ) + (1 − tk )I ]∗ [tk F ′ (ujk )∗ (F (ujk ) − z δ )
+ (1 − tk )(ujk − u0 )] ∥2
(2.22)
L. Li, C. Chen and B. Bi / Journal of Computational and Applied Mathematics 365 (2020) 112408
5
Substituting the second equation in (2.8) and (2.22) into (2.21), J
j+1
D
j+1
ξk
J
j
(u∗ , uk ) − D j (u∗ , uk ) ξk
2
∥[tk F ′ (uk )∗ F ′ (uk ) + (1 − tk )I ]∗ [tk F ′ (uk )∗ (F (uk ) − z δ ) + (1 − tk )(uk − u0 )]∥2 j (αk )2 1 − j ⟨[tk F ′ (ujk )∗ F ′ (ujk ) + (1 − tk )I ]∗ [tk F ′ (ujk )∗ (F (ujk ) − z δ ) + (1 − tk )(ujk − u0 )], ujk − u∗ ⟩ ≤
j
j
j
j
j
αk
2
=
∥tk F ′ (ujk )∗ (F (ujk ) − z δ ) + (1 − tk )(ujk − u0 )∥2
j 2 k)
(α 2
−
j
(αk )2
⟨tk F ′ (ujk )∗ (F (ujk ) − z δ ) + (1 − tk )(ujk − u0 ),
{I − [tk F ′ (ujk )∗ F ′ (ujk ) + (1 − tk )I ][tk F ′ (ujk )∗ F ′ (ujk ) + (1 − tk )I ]∗ } · [tk F ′ (ujk )∗ (F (ujk ) − z δ ) + (1 − tk )(ujk − u0 )]⟩ 1
−
∥tk F
αkj 1
′
δ
j j (uk )∗ (F (uk )
− z ) + (1 −
j j tk F ′ (uk )∗ (F (uk ) z δ ) j k j j j j tk F ′ (uk )∗ F (uk ) z F ′ (uk )(uk
+
α
[
≤[ +
2 j 2 k)
(α 1
1
αkj
]∥tk F ′ (ujk )∗ (F (ujk ) − z δ ) + (1 − tk )(ujk − u0 )∥2
⟨[tk F ′ (ujk )∗ (F (ujk ) − z δ ) + (1 − tk )(ujk − u0 )], tk F ′ (ujk )∗ [F (ujk ) − z − F ′ (ujk )(ujk − u∗ )]⟩
αkj
+[
− u∗ ) + z − z δ ] + (1 − tk )(u∗ − u0 )⟩
− −
−
− u0 )∥
(2.23) 2
+ (1 − tk )(ujk − u0 )],
−
⟨[
j tk )(uk
tk δ
(1 − tk )ρ
+
j k
α
α
j k
]∥tk F ′ (ujk )∗ (F (ujk ) − z δ ) + (1 − tk )(ujk − u0 )∥
Substituting (2.12) into (2.23), J
D
j+1
j+1
ξk
J
j
(u∗ , uk ) − D j (u∗ , uk ) ξk
≤ ∥tk F ′ (ujk )∗ (F (ujk ) − z δ ) + (1 − tk )(ujk − u0 )∥ · [(
2 j 2 k)
(α 1
1
−
αkj
)∥tk F ′ (uk )∗ (F (uk ) − z δ ) + (1 − tk )(uk − u0 )∥ j
j j t t F ′ (uk )∗ (F (uk ) j k k k k j j tk F ′ (uk )∗ (F (uk ) z δ )
+
η∥
α
≤∥ [(
j
2 j 2 k)
(α tk
+
j k
α
−
1 j k
α
+
α
− z δ + z δ − z) + (1 − tk )(ujk − u0 )∥ +
tk δ j k
α
+
(1 − tk )ρ
αkj
]
(2.24)
+ (1 − tk )(ujk − u0 )∥·
− t k ηk j k
j
)∥tk F ′ (uk )∗ (F (uk ) − z δ ) + (1 − tk )(uk − u0 )∥ j
(ηk tk + 1)δ + (1 − tk )
ρ αkj
j
j
]
By the stopping rule, ∥tk F ′ (uk )∗ (F (uk ) − z δ ) + (1 − tk )(uk − u0 )∥ > τ δ . Then (2.24) can be simplified to j
D
J
j+1
j+1
ξk
J
j
j
j
(u∗ , uk ) − D j (u∗ , uk ) ξk
≤ ∥tk F ′ (ujk )∗ (F (ujk ) − z δ ) + (1 − tk )(ujk − u0 )∥ · 2 1 tk ηk tk j j j j [( j − j + j + (ηk tk + 1))∥tk F ′ (uk )∗ (F (uk ) − z δ ) + (1 − tk )(uk )(uk − u0 )∥ (αk )2 αk αk τ αkj ρ + (1 − tk ) j ] αk
(2.25)
6
L. Li, C. Chen and B. Bi / Journal of Computational and Applied Mathematics 365 (2020) 112408 j
Here we take βk :=
2
1
−
j
(αk )2
+
j
αk
tk η k j
αk
tk
+
j
τ αk
(ηk tk + 1). Since
2
αkj ≥
(tk −
ρ 1) τ δ
+ 1 − tk ηk −
tk (ηk tk +1)
τ
we get
βkj ≤ (tk − 1)
ρ
≤0
j k
α τδ
Therefore,
βkj ∥tk F ′ (ujk )∗ (F (ujk ) − z δ ) + (1 − tk )(ujk − u0 )∥ + (1 − tk ) ≤ βkj τ δ + (1 − tk )
ρ
ρ αkj
αkj
≤0 At last, it is proved that J
D
j+1
j+1
ξk
J
j
(u∗ , uk ) − D j (u∗ , uk ) ≤ 0
(2.26)
ξk
In the same way, (2.26) holds when j = 0. Finally, we prove that the error of the iterative sequence in the Bregman distance is decreasing. □ j
J
Theorem 2. Let u0 ∈ D(F ) ∩ BV (Ω ), ξ0 ∈ L2 (Ω ). u∗ ∈ Bρ (u0 ) is the solution of F (u) = z, and uk ∈ Bρ (u0 ), D 0 (u∗ , u0k ) ≤ ξk
j∗ ukk∗
ρ2 2
.
j
We can get according to the iteration sequence {uk } in (2.8) when the stopping rule reaches k∗ = k∗ (δ ). If F satisfies (2.10)–(2.11), (2.12)–(2.14) hold, and
∥F (ujk ) − F (ukk∗ )∥ ≤ 2∥tk F ′ (ujk )∗ (F (ujk ) − z δ ) + (1 − tk )(ujk − u0 )∥ j ∗
then the number of the iteration
(2.27) jk∗
∑k∗
is limited. Taking δn > 0, as δn → 0, there exists a subsequence of {uk∗n } convergent
k=0 jk
n
jk∗
in weak∗ topology, and the limit u ∈ S(z). If S(z) ∩ Bρ (u∗ ) = {u∗ }, then uk∗n ⇀ u∗ , where k∗n is the stopping index. n
Proof. First, we prove that j ∗
J
D
j ∗
ξkk∗
∑k∗
k=0 jk
is finite. As k ≥ 0, Eq. (2.25) holds. Therefore,
J
(u∗ , ukk∗ ) − Dξ0 (u∗ , u0 ) ∗
≤
k j=jk −1 ∑ ∑ k=0
∥tk F ′ (ujk )∗ (F (ujk ) − z δ ) + (1 − tk )(ujk − u0 )∥·
j=0
[βkj ∥tk F ′ (ujk )∗ (F (ujk ) − z δ ) + (1 − tk )(ujk − u0 )∥ + (1 − tk )
(2.28)
ρ αkj
]
Since ∥tk F ′ (uk )∗ (F (uk ) − z δ ) + (1 − tk )(uk − u0 )∥ ≥ τ δ , it follows that j
j
j
∗
J
Dξ0 (u∗ , u0 ) ≥ D j
j ∗
J j
∗ k ∗ (u , uk∗ ) +
ξkk∗
k j=jk −1 ∑ ∑
−τ δ[βkj τ δ + (1 − tk )
j=0
k=0
j
ρ αkj j
Combine αk and βk , we make sure that there is a upper bound to βk τ δ + (1 − tk )
γk = − max{τ δ[βkj τ δ + (1 − tk ) j
ρ
(2.29)
] ρ j
αk
on j. Let
]} > 0
(2.30)
jk γk
(2.31)
αkj
Eq. (2.29) can be simplified as ∗
J
J
Dξ0 (u , u0 ) ≥ D ∗
(u , ∗
j ∗
ξkk∗
j ∗ ukk∗ )
+
k ∑ k=0
Obviously, the total iteration number
∑k∗
k=0 jk
is finite.
L. Li, C. Chen and B. Bi / Journal of Computational and Applied Mathematics 365 (2020) 112408
7
In the following, We prove the iterative sequence is convergent. j ∗
⟨ξk∗ , ukk∗ − u∗ ⟩ j ∗
= ⟨ξ0 , ukk∗ − u∗ ⟩ ∗
k jk −1 ∑ ∑ j − ⟨(αk )−1 [tk F ′ (ujk )∗ F ′ (ujk ) + (1 − tk )I ]∗ [tk F ′ (ujk )∗ (F (ujk ) − z δ ) k=0 j=0 j ∗
+ (1 − tk )(ujk − u0 )], ukk∗ − u∗ ⟩
(2.32)
∗
j ∗
≤ ∥ξ0 ∥ · ∥ukk∗ − u∗ ∥ −
k jk −1 ∑ ∑
(αk )−1 ∥tk F ′ (uk )∗ (F (uk ) − z δ ) + (1 − tk )(uk − u0 )∥2 j
j
j
j
k=0 j=0 k∗ jk −1
(αk )−1 ⟨tk F ′ (uk )∗ (F (uk ) − z δ ) + (1 − tk )(uk − u0 ),
∑∑
+
j
j
j
j
k=0 j=0
tk F (uk )∗ [F (uk ) − z δ − F ′ (uk )(ukk∗ − u∗ )] + (1 − tk )(uk − u0 − ukk∗ + u∗ )⟩ ′
j
j
j ∗
j
j
j ∗
For simplicity, let Lk := tk F ′ (uk )∗ (F (uk ) − z δ ) + (1 − tk )(uk − u0 ). By (2.12), the underlined part in Eq. (2.32) is j
j
j
⟨Lk , tk F ′ (ujk )∗ [F (ujk ) − z − F ′ (ujk )(ujk − u∗ )]⟩ + ⟨Lk , tk F ′ (ujk )∗ (z − z δ )⟩ j ∗
j ∗
+ ⟨Lk , (1 − tk )(ujk − u0 − ukk∗ + u∗ )⟩ + ⟨Lk , tk F ′ (ujk )∗ [−F ′ (ujk )(ukk∗ − ujk )]⟩ ≤ tk ηk ∥Lk ∥2 + tk δ (tk ηk + 1)∥Lk ∥ j ∗
j ∗
+ ⟨Lk , (1 − tk )(ujk − u0 − ukk∗ + u∗ )⟩ + ⟨Lk , tk F ′ (ujk )∗ [−F ′ (ujk )(ukk∗ − ujk )]⟩ By the nonlinear condition (2.11), j ∗
j ∗
∥F ′ (ujk )(ukk∗ − ujk )∥ ≤ (1 + ηk )∥F (ukk∗ ) − F (ujk )∥ j
Since uk ∈ Bρ (u0 ), j ∗
j ∗
∥ujk − u0 − ukk∗ + u∗ ∥ = ∥(ujk − u∗ ) + (u∗ − ukk∗ ) + (u∗ − u0 )∥ ≤ 5ρ then
⟨Lk , tk F ′ (ujk )∗ [F (ujk ) − z − F ′ (ujk )(ujk − u∗ )]⟩ + ⟨Lk , tk F ′ (ujk )∗ (z − z δ )⟩ j ∗
j ∗
+ ⟨Lk , (1 − tk )(ujk − u0 − ukk∗ + u∗ )⟩ + ⟨Lk , tk F ′ (ujk )∗ [−F ′ (ujk )(ukk∗ − u∗ )]⟩ ≤ tk ηk ∥Lk ∥ + tk δ (tk ηk + 1)∥Lk ∥ + 5ρ (1 − tk )∥Lk ∥ 2
j ∗
+ tk (1 + ηk )∥Lk ∥ · ∥F (ukk∗ ) − F (ujk )∥ Moreover, by (2.27) ,
⟨Lk , tk F ′ (ujk )∗ [F (ujk ) − z − F ′ (ujk )(ujk − u∗ )]⟩ + ⟨Lk , tk F ′ (ujk )∗ (z − z δ )⟩ j ∗
j ∗
+ ⟨Lk , (1 − tk )(ujk − u0 − ukk∗ + u∗ )⟩ + ⟨Lk , tk F ′ (ujk )∗ [−F ′ (ujk )(ukk∗ − u∗ )]⟩ ≤ tk ηk ∥Lk ∥2 + [tk δ (tk ηk + 1) + 5ρ (1 − tk )]∥Lk ∥ + 2tk (1 + ηk )∥Lk ∥2 Substituting it into Eq. (2.32), j ∗
⟨ξk∗ , ukk∗ − u∗ ⟩ ∗
≤ ∥ξ0 ∥ · ∥
j ∗ ukk∗
∗
−u ∥−
k jk −1 ∑ ∑
j
(αk )−1 [1 − 2tk − 3tk ηk ]∥Lk ∥2
k=0 j=0 k∗ jk −1
+
∑∑
j
(αk )−1 {[tk2 ηk + tk ]δ + 5ρ (1 − tk )}∥Lk ∥
k=0 j=0
By (2.26), we developed J
D
j ∗ ξkk∗
j ∗
J
(u∗ , ukk∗ ) ≤ Dξ 0 (u∗ , u0k ) ≤ k
ρ2 2
8
L. Li, C. Chen and B. Bi / Journal of Computational and Applied Mathematics 365 (2020) 112408
and we can estimate the following inequality by the Bregman distance 1
j ∗
J
D
j ∗
ξkk∗
(u∗ , ukk∗ ) ≥
2
j ∗
∥ukk∗ − u∗ ∥2
then j ∗
∥ukk∗ − u∗ ∥ ≤ ρ Therefore, j ∗
⟨ξk∗ , ukk∗ − u∗ ⟩ ∗
≤ ∥ξ0 ∥ρ −
k jk −1 ∑ ∑
j
(αk )−1 [1 − 2tk − 3tk ηk ]∥Lk ∥2
k=0 j=0 k∗ jk −1
+
∑∑
j
(αk )−1 {[tk2 ηk + tk ]δ + 5ρ (1 − tk )}∥Lk ∥
k=0 j=0
Moreover, Bregman distance D
j ∗
J j ∗
ξkk∗
(u∗ , ukk∗ ) ≥ 0, it follows that j ∗
j ∗
j ∗
J(ukk∗ ) ≤ J(u∗ ) − ⟨ξkk∗ , u∗ − ukk∗ ⟩ k∗
≤ J(u∗ ) + ∥ξ0 ∥ρ −
jk − 1 ∑∑
j
(αk )−1 [1 − 2tk − 3tk ηk ]∥Lk ∥2
k=0 j=0 ∗
+
k jk −1 ∑ ∑
j
(αk )−1 {[tk2 ηk + tk ]δ + 5ρ (1 − tk )}∥Lk ∥
(2.33)
k=0 j=0 jk∗
jk∗
n
n
We can prove ∥F (uk∗n ) − z δn ∥ → 0 as δn → 0, then ∥F (uk∗n ) − z ∥ → 0. Therefore, we can derive ∥Lk ∥ < ∞ and conclude j ∗
jk∗
that J(ukk∗ ) is uniformly bounded on small δ , and there is a subsequence of {uk∗n } is weak∗ topology convergent in BV(Ω ). n
jk∗
Moreover, F is continuous, weakly sequenced closed, then the limit u of the convergent subsequence of {uk∗n } belongs to S(z). If S(z)
⋂
jk∗
Bρ (u∗ ) = {u∗ }, then uk∗n ⇀ u∗ . n
n
□
3. Numerical analysis The principle of fluorescence molecule tomography(FMT) is that the fluorescence reporter group(such as GFP, RFP) is injected into the organism, and is excited by a specific wavelength of excited light. When the group absorbs the excited light, the energy level transition occurs to the state of high energy, and gradually attenuates back to the ground state, and the fluorescence with longer wavelength than the incident light is emitted. After excitation and fluorescence light interacts with tissue, they are captured by the surface detector. According to the data received by the probe, the mathematical model was built, and the optical parameters distribution of the tissue was reconstructed [27,28], and the high resolution image was obtained. The absorption coefficient and the scattering coefficient represent the probability that a photon is absorbed and scattered at the unit length, respectively. The unit is mm−1 or cm−1 . In this section, we will consider to reconstruct the absorption coefficient of the tissue. There are two kinds of light in fluorescent molecular tomography, one is an excitation light from the external light source(the wavelength is assumed as m), and the other is the fluorescent produced by the fluorescent markers that is injected into the target tissue under the radiation of excitation light (the wavelength is assumed as n). The chromophores and fluorophore in the tissue can absorb both kinds of light. We denote the absorption coefficients that the chromophores, the fluorophore and the tissue to the excitation light and fluorescent as µai,m , µai,n , µaf ,m , µaf ,n , µs,m , µs,n , respectively. Usually, the other optical parameters are known, and are used to reconstruct µai,m , µaf ,m . Now we consider the operator equation F (µai,m , µaf ,m ) = z δ z δ denotes the reading of the detector on the surface. Assume that there are d detectors, z δ ∈ Rd . F : D → Rd ,
(µai,m , µaf ,m ) ↦ → z δ
(3.1)
where µai,m , µaf ,m denote the absorption coefficient of fluorescence by chromoplast and fluorophore. D = {(µai,m , µaf ,m ) ∈ L∞ (X ) × L∞ (X )}.
L. Li, C. Chen and B. Bi / Journal of Computational and Applied Mathematics 365 (2020) 112408
9
Fig. 1. Light enters the columnar tissue.
The direct problem is to solve the reading of the detector on surface z δ by the absorption coefficients µai,m , µaf ,m of chromophore and fluorophore for fluorescence. On the contrary, the inverse problem is to solve µai,m , µaf ,m by z δ . Later, we will use Eq. (2.8) to solve our problem. In Sections 3.2 and 3.3, we will give the models for solving the direct problem and inverse problem of FMT, respectively. 3.1. Radiation Transfer Equation for FMT Radiation Transfer Equation (RTE) can effectively simulate the propagation of photons in tissues and is widely used in biomedical imaging. Radiation transfer equation is used to simulate the scattering, reflection and absorption of light in tissue, and uses radiation rate and phase function p(sˆ, sˆ′ ) to describe the intensity of light and the probability distribution of photon scattering, respectively. Radiation rate φ (r , sˆ) represents the flux density of average power in a given direction sˆ and in unit solid angle, and the unit is W m−2 sr−1 Hz−1 . The phase function p(sˆ, sˆ′ ) represents the probability of being incident on a scattering particle from the direction sˆ′ , and then being scattered into the direction sˆ. p(sˆ, sˆ′ )dΩ ′ describes the probability that a photon incident from the direction sˆ′ is scattered into some solid angle in the direction sˆ,
∫ 4π
p(sˆ, sˆ′ )dΩ ′ = 1
It is generally assumed that the phase function is only related to the angle between the incident and exit directions, p(sˆ, sˆ′ ) = p(sˆ′ , sˆ) = p(cosθ ) The Henyey–Greenstein function is chosen as the phase function of the model, i.e., p(cos θ ) =
1 − g2
1
4π (1 + g 2 − 2gcosθ ) 32
where g is an anisotropic coefficient, and is defined as
∫ g = ⟨cosθ⟩ =
4π
p(sˆ, sˆ′ )(sˆ · sˆ′ )dΩ ′
∫
4π
p(sˆ, sˆ′ )dΩ ′
∫ = 4π
p(sˆ, sˆ′ )(sˆ, sˆ′ )dΩ ′
Examine a columnar structure with a length of ds and a unit cross section as shown in Fig. 1. Light having an emissivity of φ (r , sˆ) is incident on the cylindrical voxel. The reason for the change in emissivity can be reduced to the following aspects: (1) The emissivity will disappear due to absorption and scattering by the tissue. The reduction in emissivity can be expressed as: dφ (r , sˆ) = −(µa + µs )dsφ (r , sˆ) (2) The energy obtained by external scattering. This part of the increment can be expressed as dφ (r , sˆ) = µs ds
∫ 4π
p(sˆ, sˆ′ )φ (r , sˆ′ )dsˆ′
(3) If the light source emits photon in the body, the emissivity will also be increased. Let the power radiated from the internal light source along the direction sˆ in unit angle to be q(r , sˆ). The increment is dφ (r , sˆ) = q(r , sˆ)ds
10
L. Li, C. Chen and B. Bi / Journal of Computational and Applied Mathematics 365 (2020) 112408
Fig. 2. Angular dispersion diagram.
These three factors are considered, the radiation transfer equation can be derived by the law of conservation of energy as dφ (r , sˆ) ds
= −(µa + µs )φ (r , sˆ) + µs
∫ 4π
p(sˆ, sˆ′ )φ (r , sˆ′ )dsˆ′ + q(r , sˆ)
which is also written as Hamilton operator equation: (s · ∇ + µa (r) + µs (r))φ (r , sˆ) = µs (r)
∫ 4π
p(sˆ, sˆ′ )φ (r , sˆ′ )dsˆ′ + q(r , sˆ)
RTE can reflect the energy balance relation of every points in the media. Fluorescence Molecule Tomography (FMT) includes two radiation transfer equations, which describe propagation processes that the fluorescence produced from exciting light and fluorescent groups in the tissue. The propagation equation of exciting light is (s · ∇ + µai,m (r) + µaf ,m (r) + µs,m (r))φm (r , sˆ) = µs,m (r)
∫ 4π
p(sˆ, sˆ′ )φm (r , sˆ′ )dsˆ′
The propagation equation of fluorescent is (s · ∇ + µai,n (r) + µaf ,n (r) + µs,n (r))φn (r , sˆ)
= µs,n (r)
∫ 4π
p(sˆ, sˆ′ )φn (r , sˆ′ )dsˆ′ + ηµaf ,m (r)
∫ 4π
φm (r , sˆ′ )dsˆ′
The boundary condition is
φm (r , sˆ) = gm (r , sˆ),
φn (r , sˆ) = 0 on Γ−
Here φm (r , sˆ′ ) and φn (r , sˆ′ ) denote the emissivity of exciting light and fluorescent, and η is determined by the property of fluorescent group in the tissue. Γ− is incident boundary. 3.2. Direct model of FMT The direct problem of FMT is to solve the detector readings z δ on the body by using the absorption coefficients µai,m and µaf ,m of fluorescence by chromoplast and fluorophore. Aiming at 2D system, we use the finite element method to discrete the space to get the matrix equation, and then to solve this equation by SOR iteration method. At last, we construct the complete direct model. First of all, the data is spatially discrete and the region is triangulated. The region Ω is decomposed into N triangular areas, which are discretized into I × J points r = (xi , yj ), i = 0, . . . , I, j = 0, . . . , J. Then the angle sˆ is discretized, and [0, 2π ) is evenly divided into K directions. The length of each cell is △θ , which refers to Fig. 2. Here sˆ = θk , k = 1, . . . , K . The radiation rate can be approximated in sections:
φ (xi , yj , θ ) ≈
K ∑ k=1
φ k,i,j Lk (θ )
L. Li, C. Chen and B. Bi / Journal of Computational and Applied Mathematics 365 (2020) 112408
11
φ k,i,j , k = 1, . . . , K is the radiation rate along the kth direction at (xi , yj ), Lk (θ ) is piecewise linear function of θ , { 1, k = k′ Lk (θk′ ) = 0, others After dispersion, the RTE equation of fluorescence molecular tomography can be reduced to the following equations:
∂φmk,i,j ∂φmk,i,j + sin θk + (µai,m + µaf ,m + µs,m )φmk,i,j ∂x ∂y ∫ k,i,j ′ (θ )dθ ′ = µs,m p(θk , θ ′ )φm
cos θk
2π
∂φnk,i,j ∂φnk,i,j cos θk + sin θk + (µai,n + µaf ,n + µs,n )φnk,i,j ∂x ∂y ∫ ∫ = µs,n p(θk , θ ′ )φnk,i,j (θ ′ )dθ ′ + ηµaf ,m φmk,i,j (θ ′ )dθ ′ 2π
φmk,i,j = gmk,i,j ,
2π
φnk,i,j = 0 on Γ−
The integral on the right-hand side of the equation is expressed as
∫ 2π
p(θk , θ ′ )φ k,i,j (θ ′ )dθ ′ =
K ∑
αk′ pk,k′ φ k ,i,j and ′
k′ =1
∫ 2π
φ k,i,j (θ ′ )dθ ′ =
K ∑
βk′ φ k ,i,j ′
k′ =1
αk′ and βk′ are selected weight. Then we use the upwind difference scheme to approximate the partial derivative of the radiation rate φ at (xi , yj ). Then 2D RTE is obtained after the angle is discretized: φmk,i,j − φmk,i−1,j φmk,i,j − φmk,i,j−1 + sin θk + (µai,m + µaf ,m + µs,m )φmk,i,j △x △y K ∑ ′ = µs,m αk′ pk,k′ φmk ,i,j
cos θk
k′ =1 k,i,j n
φnk,i,j − φnk,i,j−1 − φnk,i−1,j + sin θk + (µai,n + µaf ,n + µs,n )φnk,i,j △x △y K K ∑ ∑ ′ ′ = µs,n αk′ pk,k′ φnk ,i,j + ηµaf ,m βk′ φmk ,i,j , cos θk
φ
k′ =1
k′ =1
φmk,i,j = gmk,i,j ,
φnk,i,j = 0 on Γ−
(3.2)
Obviously, Eq. (3.2) can be written as matrix equation: Aφ = b. Use SOR method, the following iteration algorithm is obtained:
[φmk,i,j ]l+1 = ρ
µs,m
∑K
k′ =1
αk′ pk,k′ [φmk ,i,j ]l + [φmk,i−1,j ]l cos θk /△x + [φmk,i,j−1 ]l sin θk /△y cos θk /△x + sin θk /△y + (µai,m + µaf ,m + µs,m ) ′
+ (1 − ρ )[φmk,i,j ]l [φnk,i,j ]l+1 =ρ
µs,n
(3.3)
∑K
k′ =1
αk′ pk,k′ [φ
k′ ,i,j l n
∑K
k′ ,i,j l m
k,i−1,j l n
k,i,j−1 l n
] + ηµaf ,m k′ =1 βk′ [φ ] + [φ ] cos θk /△x + [φ cos θk /△x + sin θk /△y + (µai,n + µaf ,n + µs,n )
] sin θk /△y
+ (1 − ρ )[φnk,i,j ]l where ρ is super relaxation coefficient. The detector reading at the boundary lattice point can be calculated according to the following formula
Φ (xi , yj ) =
∫
φ (xi , yj , sˆ)dsˆ ≈ AP
k2 ∑
ak φk,i,j = Φij
k=k1
The weight ak can be obtained by generalized Trapezoidal rule. Therefore, we will get the corresponding detector readings F (µai,m , µaf ,m ) when the optical parameters µai,m and µaf ,m are given arbitrary.
12
L. Li, C. Chen and B. Bi / Journal of Computational and Applied Mathematics 365 (2020) 112408
Fig. 3. Subdivision graph of tissue.
Fig. 4. True absorption coefficient and scattering coefficient distribution.
3.3. Inverse model of FMT The inverse problem of FMT is to use the reading z δ of known detectors to inverse the absorption coefficients µai,m and µaf ,m of chromosome and fluorophore to fluorescence. Set u = (µai,m , µaf ,m ). The nonlinear equation of FMT imaging can be written as F (u) = z δ . There are d detectors, z δ ∈ Rd . Applying the iteration algorithm (2.8) to inverse the absorption coefficients µai,m and the scattering coefficients µaf ,m of the chromophores to the excited light. 3.4. Numerical results Consider the 2D tissue (Fig. 3). The graph is discretized into 276 triangulars, 32 directional angles and 156 vertices. Place 16 excitation light sources and 16 receivers at regular intervals along the region boundary. Two abnormal bodies are placed in the tissue (refer to Fig. 4(a) and (b)). The absorption coefficients of the anomalous bodies are 0.4, and the absorption coefficient of background of the tissue is 0.2. The scattering coefficients of the anomalous bodies are 4, and the scattering coefficient of the background of the tissue is 2. We only focus on the inversion of absorption coefficient of one abnormal body (Fig. 4(a)).
L. Li, C. Chen and B. Bi / Journal of Computational and Applied Mathematics 365 (2020) 112408
13
Fig. 5. (a)–(d). The reconstructions of four methods for the absorption coefficient of the tissue.
Table 1 Comparison of four regularization methods.
Iteration steps Running time (s) Relative error
Tikhonov regularizatioin
l1 regularization
Total variational regularization
New regularization
24 578.34 0.2864
9 209.4288 0.2378
15 349.79 0.2541
12 330.62 0.2338
The simulation results of the absorption coefficient are obtained by using Tikhonov regularization, regularization with l1 norm, total variational regularization, and the new regularization method proposed in this paper, respectively. As shown in Fig. 5(d), the proposed method can exactly inverse the position of abnormal body, and depict clearly the boundary information of abnormal body. Compared with Tikhonov regularization and total variational regularization, the proposed new method runs faster, and the inversion of the abnormal body is closer to the true value. Furthermore, the inversion effect of the proposed new method on the background is good. From Table 1 and Fig. 5, l1 regularization method is the fastest, and the inversion effect on the abnormal body and background are better in these four methods. However, the proposed new method has the smallest relative error. Tikhonov regularization method is relatively worse in both computation time and inversion effect.
14
L. Li, C. Chen and B. Bi / Journal of Computational and Applied Mathematics 365 (2020) 112408
4. Conclusion Total variational regularization can effectively inverse the piecewise constant function, and can preserve the boundary of large object, thus are widely used to the field of image processing. Moreover, homotopy method can enlarge the range of iterative initial value. Therefore, a new total variational regularization method is proposed based on total variational regularization and fix-point homotopy. Moreover, the convergence of the new method is proved in theory. In the stage of numerical simulation, the absorption coefficient inversion problem of fluorescence molecular tomography (FMT) was solved, and the radiation transfer equation (RTE) was used to simulate the propagation of photons in tissues. Experimental results show that the new algorithm can accurately inverse the position of anomalous body and clearly depict the edge information of anomalous body. Although four regularization methods can calculate the position of abnormal body, Tikhonov regularization method and total variational regularization are not very good at restoring the image background. In contrast, l1 method and the new approach better preserve the background information, their whole simulation results are closer to the true value, and their accuracy are also slightly higher. However, the relative error of the proposed new method is smallest. Therefore, the proposed new method is a stable and effective regularization method. References [1] V. Ntziachristos, C. Bremer, E.E. Graves, et al., In vivo tomographic imaging of near-infrared fluorescent probes, Mol. Imaging 1 (2) (2002) 82–88. [2] V. Ntziachristos, C.H. Tung, C. Bremer, et al., Fluorescence molecular tomography resolves protease activity in vivo, Nature Med. 8 (7) (2002) 757–760. [3] S. Chandrasekhar, Radiative Transfer, Clarendon, London, 1950. [4] K.M. Case, P.F. Zweifel, Linear Transport Theory, Addison-Wesley, Reading, MA, 1967. [5] V. Agoshkov, Boundary Value Problems for Transport Equations, Birkhauser, Boston, MA, 1998. [6] G. Choppin, J.-O. Liljenzin, J. Rydberg, Radiochemistry and Nuclear Chemistry, third ed., Butterworth-Heinemann, Woburn, MA, 2002. [7] W.X. Cong, L.V. Wang, G. Wang, Formulation of photon diffusion from spherical bioluminescent sources in an infinite homogeneous medium, Biomed. Eng. Online 3 (1) (2004) 1–6. [8] C.H. Contag, B.D. Ross, It’s not just about anatomy: in vivo bioluminescence imaging as an eyepiece into biology, J. Mag. Reson. Imaging 16 (4) (2002) 378–387. [9] A.D. Klose, U. Netz, J. Beuthan, et al., Optical tomography using the time-independent equation of radiative transfer−Part 1: forward model, J. Quant. Spectrosc. Radiat. Transfer 72 (5) (2002) 691–713. [10] G. Hao, K.Z. Hong, A fast forward solver of radiative transfer equation, Transp. Theory Stat. Phys. 38 (3) (2009) 149–192. [11] B.C. Wilson, G. Adam, A Monte Carlo model for the absorption and flux distributions of light in tissue, Med. Phys. 10 (6) (1983) 824. [12] L. Wilson, S.L. Jacques, L. Zheng, MCML-Monte Carlo Modeling of light transport in multi-layered tissues, Comput. Methods Programs Biomed. 47 (2) (1995) 131–146. [13] D. Boas, J. Culver, J. Stott, et al., Three dimensional Monte Carlo code for photon migration through complex heterogeneous media including the adult human head, Opt. Express 10 (3) (2002) 159–170. [14] A. Kienle, M.S. Patterson, Improved solutions of the steady-state and the time-resolved diffusion equations for equations for reflectance from a semi-infinite turbid medium, J. Opt. Soc. Amer. A 14 (1) (1997) 246–254. [15] H. Egger, M. Schlottbom, Analysis and regularization of problems in diffuse optical tomography, SIAM J. Math. Anal. 42 (2010) 1934–1948. [16] A.D. Klose, A.H. Hielscher, Iterative reconstruction scheme for optical tomography based on the equation of radiative transfer, Med. Phys. 26 (8) (1999) 1698–1707. [17] K. Ren, T. Ding, S. Vallelian, A one-step reconstruction algorithm for quantitative photoaccoustic imaging, Inverse Problems 31 (2015) 095005. [18] A.N. Tikhonov, V.Y. Arsenin, Solutions of Ill-Posed Problems, Wiley, Washington, DC, 1977. [19] M. Hanke, A. Neubauer, O. Scherzer, A convergence analysis of the Landweber iteration for nonlinear ill-posed problems, Numer. Math. 72 (1995) 21–37. [20] Parimah Kazemi, Robert J. Renka, A Levenberg–Marquardt method based on Sobolev gradients, Nonlinear Anal. TMA 75 (16) (2012) 6170–6179. [21] H. Gao, H.K. Zhao, Multilevel bioluminescence tomography based on radiative transfer equation, Part 2: total variation and l1 data fidelity, Opt. Express 18 (3) (2010) 2894. [22] D. Hyde, E.L. Miller, D.H. Brooks, et al., Data specific spatially varying regularization for multimodal fluorescence molecular tomography, IEEE Trans. Med. Imaging 29 (2) (2010) 365–374. [23] J. Dutta, S. Ahn, C. Li, et al., Joint L1 and total variation regularization for fluorescence molecular tomography, Phys. Med. Biol. 57 (6) (2012) 1459–1476. [24] Joyita Dutta, Snagtae Ahn, Changqing Li, Simon R. Cherry, Richard M. Leahy, Joint L1 and total variation regularization for fluorescnece molecular tomography, Phys. Med. Biol. 57 (6) (2012) 1459–1476. [25] B. Bi, B. Han, W. Han, J. Tang, L. Li, Image reconstruction for diffuse optical tomography based on radiative transfer equation, Comput. Math. Methods Med. 2015 (2015). [26] B. Bi, J. Tang, Diffuse optical tomography based on radiative transfer equation with L∞ data fidelity and total variation penalty regularization, J. Med. Imaging Health Inf. 7 (6) (2017) 1212–1218. [27] I. Georgakoudi, B.C. Jacobson, D.J. Van, et al., Fluorescence, reflectance, and light-scattering spectroscopy for evaluating dysplasia in patients with Barrett’s esophagus, Gastroenterology 120 (7) (2001) 1620–1629. [28] Vasilis Ntziachristos, Jorge Ripoll, Lihong V. Wang, Ralph Weissleder, Looking and listening to light: the evolution of whole-body photonic imaging, Nature Biotechnol. 23 (2005) 313–320.