# Mixed Total Variation andL1 Regularization Method for Optical Tomography Based on Radiative Transfer Equation

Authors: Jinping Tang; Bo Han; Weimin Han; Bo Bi; Li Li
Journal: Computational and Mathematical Methods in Medicine
 (2017)

---

## Abstract

Optical tomography is an emerging and important molecular imaging modality. The aim of optical tomography is to reconstruct optical properties of human tissues. In this paper, we focus on reconstructing the absorption coefficient based on the radiative transfer equation (RTE). It is an ill-posed parameter identification problem. Regularization methods have been broadly applied to reconstruct the optical coefficients, such as the total variation (TV) regularization and theL1 regularization. In order to better reconstruct the piecewise constant and sparse coefficient distributions, TV and L1 norms are combined as the regularization. The forward problem is discretized with the discontinuous Galerkin method on the spatial space and the finite element method on the angular space. The minimization problem is solved by a Jacobian-based Levenberg-Marquardt type method which is equipped with a split Bregman algorithms for the L1 regularization. We use the adjoint method to compute the Jacobian matrix which dramatically improves the computation efficiency. By comparing with the other imaging reconstruction methods based on TV and L1 regularizations, the simulation results show the validity and efficiency of the proposed method.


---

## Body

## 1. Introduction

Optical tomography with near-infrared light is a promising technique for noninvasive studying of the functional characters of human tissues. One can find many applications of this technique, for example, the early detection of breast cancer, cervical cancer screening, monitoring of infant brain tissue oxygenation level and functional brain activation, and the study of dosimetry in photon dynamic therapy; see [1–4]. Unlike the X-ray tomography, optical tomography uses the near-infrared (NIR) light as the probing radiation. Different human tissues have different absorption and scattering properties which affect the transmission of the NIR light. Therefore, we can identify the location and quantity of abnormal tissues by the emerging light measured on the boundary.Optical tomography is a high contrast imaging modality. However, currently only low resolution reconstructions are possible, especially when using an unmodulated continuous wave source [5]. Another barrier of optical tomography is that reconstructed images have a poor quality particularly when abnormal targets are located deep in tissues. On the other hand, due to the limited measurement data, as well as the high scattering and high absorption properties, the identification problem of optical tomography usually is ill-posed and underdetermined. The ill-posedness makes the coefficient reconstruction sensitive to small perturbation from measurements, such as noise and computational error. Therefore, various reconstruction algorithms based on regularization are developed to obtain reasonable and stable reconstructions [6–9]. Recently, the regularization method with L1 norm has been used in optical tomography [10, 11]. In [10], Levenberg-Marquardt strategy is applied for solving the L2 regularization step of split Bregman algorithm. The contrast tests show the superiority of L1 regularization over L2 regularization. The numerical experiments therein also show that the L1 regularization has a better utility when the independent measurements are much more limited. In [11], linearized Bregman iteration based on the Bregman distance is exploited to minimize the sparse regularization. The experimental results in numerical simulation of anin vivo mouse demonstrate the effectiveness of the algorithm. However, the sparsity regularization may oversparsify the distribution of the coefficients.Among many types of regularizations, TV regularization is often used due to its abilities of well reconstructing the discontinuous and piecewise constant distributions. TV regularization was first introduced in the field of image processing and image reconstruction; see [12–14] and the references therein. To date, TV regularization has been the preferred regularization strategy of considerable recent works, like photoacoustic tomography (PAT) [15–17], bioluminescence tomography (BLT) [18, 19], fluorescence tomography (FT) [20], optical coherence tomography (OCT) [21], and so on. In [18], the split Bregman algorithm for TV regularization (SBRTV) is applied to the reconstruction of source distribution in BLT. The algorithm is evaluated with 2D and 3D simulations and 3Din vivo experiments. In the reported results, the source distribution can be reconstructed with better accuracy for the location with SBRTV regularization than with L2 or L1 regularizations. In [19], the synergistic combination of Bregman method and TV regularization is utilized to quantitatively improve the reconstruction of absorption and scattering coefficients for both Jacobian-based and gradient-based methods in quantitative PAT. The feasibility of the algorithm is checked with 3D simulations.For optical tomography, the coefficients are usually taken as piecewise constant. Since TV regularization is effective for piecewise constant reconstruction, it is a natural choice for reconstructing coefficients of piecewise constant distributions in optical tomography. However, due to the nonsmoothness and nondifferentiability, TV regularization is difficult to realize computationally. TV regularization for optical tomography based on the radiative transport equation has been studied in [22], where inexact Gauss-Newton method is used to solve the TV regularization problem. In [23], W1,2 norm is chosen as the penalty in the regularization strategy for optical tomography based on the frequency radiative transport equation.In this paper, we use regularization as a combination of TV andL1 norm. To reduce the computational complexity, we apply the split Bregman method to solve the joint regularization problem. It is reasonable to combine L1 regularization and TV regularization to improve the reconstruction quality in optical tomography. Moreover, based on the split Bregman method, the iterative algorithms can be computationally efficient. Various experiments in 2D are performed to evaluate the performance of the algorithm.The rest of this paper is organized as follows. In Section2, we briefly describe the forward and inverse problems for optical tomography. In Sections 3 and 4, we show the implementation details about computing the Fréchet derivative of the forward operator and the iterative procedure. In Section 5, numerical results are presented.


## 2. Optical Tomography

We consider two mathematical problems: the forward problem and the inverse problem. For the forward problem, based on the physical model of light propagation in tissues, for a given set of optical properties, we model the measurements on the boundary. For the inverse problem, the optical properties can be reconstructed by matching the predictions calculated from the forward problem and the measurements from the detectors.


### 2.1. Forward Problem

The light propagation in tissues is described by the radiative transport equation(1)ω·∇ux,ω+μax+μsxux,ω=μsx∫Ωkω·ω^ux,ω^dσω^inX×Ω.Here, X⊂Rd, d=2 or 3, denotes a bounded convex domain with a C1 boundary ∂X and Ω≔Sd-1 denotes the unit sphere of Rd. The variables x and ω denote the spatial position and the angular direction. u(x,ω) describes the density of photons. The expression ω·∇u(x,ω) denotes the directional derivative at position x along the direction ω. The nonnegative normalized phase function k(ω·ω^) is the probability that photons traveling in the direction ω^ are scattered into the direction ω. In optical tomography, the phase function usually is taken as the Henyey-Greenstein phase function (cf. [24]): in two dimensions, it is of the form (2)kω·ω^=1-g22π1+g2-2gω·ω^,where the parameter g∈(-1,1) is the anisotropy factor of the scattering medium. The absorption coefficient is denoted by μa(x), and the scattering coefficient is denoted by μs(x). For the optical parameters μa(x) and μs(x), the following conditions are assumed to hold throughout this presentation.Assumptions (A1)
The functionμa(x) is uniformly positive and bounded; that is, there exist two positive constants μa1 and μa2 such that 0<μa1≤μa≤μa2<∞ a.e. in X. (A2)
The functionμs(x) is uniformly positive and bounded; that is, there exist two positive constants μs1 and μs2 such that 0<μs1≤μs≤μs2<∞ a.e. in X.Equation (1) is supplemented by boundary conditions. Similar to the X-ray CT, optical tomography experiments acquire the current distribution of detectors on the boundary under multi-incidents. Let ζi, 1≤i≤s, be disjoint, connected subsets of ∂X. Corresponding to s incident sources uin,i on ζi,1≤i≤s, define ui(x,ω) by(3)ω·∇uix,ω+μax+μsxuix,ω=μsx∫Ωkω·ω^uix,ω^dσω^inX×Ω,uix,ω=uin,ix,ω,x∈ζi,ω·νx<0,0,otherwise,on∂X×Ω,where ν(x) is the unit outward normal vector at x∈∂X.Corresponding to each independent incidentuin,i(x,ω), the measurable quantity is the outgoing light Mi(x) on the boundary of domain, which can be written as(4)Mix=∫ω·νx>0ω·νxuix,ωdσω,x∈∂X.If we assume there are d detectors located at d different positions ξj,1≤j≤d, then the optical tomography experiment consists of exciting the domain X with a sequence of incident source uin,i and recording the corresponding measurements data Mj,i=Mixx=ξj. Here, j and i represent the row index and column index in matrix Mj,i, respectively.With above notations, a mathematical description of such an experiment is provided by the following nonlinear forward operator:(5)Fi:D⟶L2∂X,μa,μs⟼Mix,1≤i≤s,which maps prescribed optical parameters to the corresponding measurements data. Here, Fi denotes the ith forward operator corresponding to the ith incident source and the resulting measurement data on d detectors. The forward operator Fi is well defined for μa and μs in the set(6)D=μa,μssatisfying assumptions  (A1)-(A2)(cf. [25]).There are many references on the discretization of RTE; see, for instance, [26–30]. In this paper, we use continuous linear elements and discontinuous Galerkin method with piecewise linear functions to discretize the angular variables and the spatial variables, respectively. An upwind numerical flux is used to approximate the incoming flux through the surface of the control element and inflow boundary. After assembling the full discretization formulation and forming a system of tebra equation, we solve the linear system by Gauss-Seidel method. We use piecewise constants to approximate the absorption and scattering coefficients. We can express (μa,μs) as(7)μax≈∑k=1Nμa,kχkx,μsx≈∑k=1Nμs,kχkx,where N denotes the number of the elements, χk(x) denotes the character function corresponding to the kth element, and μa,k and μs,k denote the values of absorption and scattering coefficients on the kth element.


### 2.2. Inverse Problem

The inverse problem of optical tomography is to determine the unknown coefficientsμa and μs from the boundary detector readings. In this paper, we only reconstruct the absorption coefficient assuming that the scattering coefficient is known; then the forward operator Fi in fact acts on the unknown μa only. Thus, the inverse problem of optical tomography is to determine μa from the following system of nonlinear equations:(8)Fiμa=Mi,1≤i≤s.Then the optical tomography can be formulated as minimizing the difference between the measurement data and the model predictions(9)μa=arg⁡minD⁡12∑i=1sFiμa-MiL2∂X2,known as data fidelity. The inverse problem (9) is ill-posed in the sense that the amount of the measurements is quite limited compared with the number of the unknowns and that the measurements contain noises. To overcome the ill-posedness, the data fidelity should be combined with appropriate regularization. In the regularization strategy, we minimize the following objective functional:(10)μa=arg⁡minD⁡12∑i=1sFiμa-MiL2∂X2+αRμa.Here, R(μa) is the regularization penalty functional that enforcesa prior information on μa and α>0 is the regularization parameter that trades off the weight between the discrepancy term and the penalty functional.According to different ways of treating the derivative of data fidelity, there are two sorts of approaches for solving (10). One is to linearize the forward operator Fi near the nth iteration of μa, which is denoted as μan, as follows:(11)Fiμa≈Fiμan+Fi′μanμa-μan,where Fi′(μan) is the derivative of the forward operator Fi with respect to μan. The linearized formulation provides a good approximation when μan is close to the true value. As a result, the minimized problem (10) is treated by an iterative procedure as follows:(12)μan+1=arg⁡minD12∑i=1sFi′μanμa-μan+Fiμan-MiL2∂X2+αRμan.This Jacobian-based minimization problem can be solved with many iteratively optimization techniques, such as the Levenberg-Marquardt type method (LM). LM type method is the special case of Gauss-Newton type method. The standard LM method (or Gauss-Newton method) uses μa-μanL2(X)2 as the penalty. The LM type method defines the update μan+1 in a region around μan, while the Gauss-Newton method always defines μan+1 in a neighbourhood of the initial guess μa0 for each n≥0. Hence, from the optimization point of view, LM type method is more favourable in nature. Based on the Jacobian-based method, many optimization techniques such as split Bregman method can be used to solve (12). On the other hand, when μan is close to the true value, the Jacobian-based method shares the typical quadratic convergence from Newton method.Alternatively, (10) can be solved directly by minimizing the nonlinear functional with some gradient-based methods, like Quasi-Newton method, nonlinear conjugate gradient method, limited-memory BFGS method, and so on. Thus, the computation of the linearized Jacobian matrix is avoided and only the gradient of the nonlinear minimizing functional needs to be computed. The gradient-based method in general has superlinear convergence and is economic in memory storage which is suitable for large scale problems, such as 3D problems. Since our numerical experiments are all done in two dimensions and the problem scale is not so large as that in three dimensions, we use the Jacobian-based LM type method to solve (10).


## 2.1. Forward Problem

The light propagation in tissues is described by the radiative transport equation(1)ω·∇ux,ω+μax+μsxux,ω=μsx∫Ωkω·ω^ux,ω^dσω^inX×Ω.Here, X⊂Rd, d=2 or 3, denotes a bounded convex domain with a C1 boundary ∂X and Ω≔Sd-1 denotes the unit sphere of Rd. The variables x and ω denote the spatial position and the angular direction. u(x,ω) describes the density of photons. The expression ω·∇u(x,ω) denotes the directional derivative at position x along the direction ω. The nonnegative normalized phase function k(ω·ω^) is the probability that photons traveling in the direction ω^ are scattered into the direction ω. In optical tomography, the phase function usually is taken as the Henyey-Greenstein phase function (cf. [24]): in two dimensions, it is of the form (2)kω·ω^=1-g22π1+g2-2gω·ω^,where the parameter g∈(-1,1) is the anisotropy factor of the scattering medium. The absorption coefficient is denoted by μa(x), and the scattering coefficient is denoted by μs(x). For the optical parameters μa(x) and μs(x), the following conditions are assumed to hold throughout this presentation.Assumptions (A1)
The functionμa(x) is uniformly positive and bounded; that is, there exist two positive constants μa1 and μa2 such that 0<μa1≤μa≤μa2<∞ a.e. in X. (A2)
The functionμs(x) is uniformly positive and bounded; that is, there exist two positive constants μs1 and μs2 such that 0<μs1≤μs≤μs2<∞ a.e. in X.Equation (1) is supplemented by boundary conditions. Similar to the X-ray CT, optical tomography experiments acquire the current distribution of detectors on the boundary under multi-incidents. Let ζi, 1≤i≤s, be disjoint, connected subsets of ∂X. Corresponding to s incident sources uin,i on ζi,1≤i≤s, define ui(x,ω) by(3)ω·∇uix,ω+μax+μsxuix,ω=μsx∫Ωkω·ω^uix,ω^dσω^inX×Ω,uix,ω=uin,ix,ω,x∈ζi,ω·νx<0,0,otherwise,on∂X×Ω,where ν(x) is the unit outward normal vector at x∈∂X.Corresponding to each independent incidentuin,i(x,ω), the measurable quantity is the outgoing light Mi(x) on the boundary of domain, which can be written as(4)Mix=∫ω·νx>0ω·νxuix,ωdσω,x∈∂X.If we assume there are d detectors located at d different positions ξj,1≤j≤d, then the optical tomography experiment consists of exciting the domain X with a sequence of incident source uin,i and recording the corresponding measurements data Mj,i=Mixx=ξj. Here, j and i represent the row index and column index in matrix Mj,i, respectively.With above notations, a mathematical description of such an experiment is provided by the following nonlinear forward operator:(5)Fi:D⟶L2∂X,μa,μs⟼Mix,1≤i≤s,which maps prescribed optical parameters to the corresponding measurements data. Here, Fi denotes the ith forward operator corresponding to the ith incident source and the resulting measurement data on d detectors. The forward operator Fi is well defined for μa and μs in the set(6)D=μa,μssatisfying assumptions  (A1)-(A2)(cf. [25]).There are many references on the discretization of RTE; see, for instance, [26–30]. In this paper, we use continuous linear elements and discontinuous Galerkin method with piecewise linear functions to discretize the angular variables and the spatial variables, respectively. An upwind numerical flux is used to approximate the incoming flux through the surface of the control element and inflow boundary. After assembling the full discretization formulation and forming a system of tebra equation, we solve the linear system by Gauss-Seidel method. We use piecewise constants to approximate the absorption and scattering coefficients. We can express (μa,μs) as(7)μax≈∑k=1Nμa,kχkx,μsx≈∑k=1Nμs,kχkx,where N denotes the number of the elements, χk(x) denotes the character function corresponding to the kth element, and μa,k and μs,k denote the values of absorption and scattering coefficients on the kth element.


## 2.2. Inverse Problem

The inverse problem of optical tomography is to determine the unknown coefficientsμa and μs from the boundary detector readings. In this paper, we only reconstruct the absorption coefficient assuming that the scattering coefficient is known; then the forward operator Fi in fact acts on the unknown μa only. Thus, the inverse problem of optical tomography is to determine μa from the following system of nonlinear equations:(8)Fiμa=Mi,1≤i≤s.Then the optical tomography can be formulated as minimizing the difference between the measurement data and the model predictions(9)μa=arg⁡minD⁡12∑i=1sFiμa-MiL2∂X2,known as data fidelity. The inverse problem (9) is ill-posed in the sense that the amount of the measurements is quite limited compared with the number of the unknowns and that the measurements contain noises. To overcome the ill-posedness, the data fidelity should be combined with appropriate regularization. In the regularization strategy, we minimize the following objective functional:(10)μa=arg⁡minD⁡12∑i=1sFiμa-MiL2∂X2+αRμa.Here, R(μa) is the regularization penalty functional that enforcesa prior information on μa and α>0 is the regularization parameter that trades off the weight between the discrepancy term and the penalty functional.According to different ways of treating the derivative of data fidelity, there are two sorts of approaches for solving (10). One is to linearize the forward operator Fi near the nth iteration of μa, which is denoted as μan, as follows:(11)Fiμa≈Fiμan+Fi′μanμa-μan,where Fi′(μan) is the derivative of the forward operator Fi with respect to μan. The linearized formulation provides a good approximation when μan is close to the true value. As a result, the minimized problem (10) is treated by an iterative procedure as follows:(12)μan+1=arg⁡minD12∑i=1sFi′μanμa-μan+Fiμan-MiL2∂X2+αRμan.This Jacobian-based minimization problem can be solved with many iteratively optimization techniques, such as the Levenberg-Marquardt type method (LM). LM type method is the special case of Gauss-Newton type method. The standard LM method (or Gauss-Newton method) uses μa-μanL2(X)2 as the penalty. The LM type method defines the update μan+1 in a region around μan, while the Gauss-Newton method always defines μan+1 in a neighbourhood of the initial guess μa0 for each n≥0. Hence, from the optimization point of view, LM type method is more favourable in nature. Based on the Jacobian-based method, many optimization techniques such as split Bregman method can be used to solve (12). On the other hand, when μan is close to the true value, the Jacobian-based method shares the typical quadratic convergence from Newton method.Alternatively, (10) can be solved directly by minimizing the nonlinear functional with some gradient-based methods, like Quasi-Newton method, nonlinear conjugate gradient method, limited-memory BFGS method, and so on. Thus, the computation of the linearized Jacobian matrix is avoided and only the gradient of the nonlinear minimizing functional needs to be computed. The gradient-based method in general has superlinear convergence and is economic in memory storage which is suitable for large scale problems, such as 3D problems. Since our numerical experiments are all done in two dimensions and the problem scale is not so large as that in three dimensions, we use the Jacobian-based LM type method to solve (10).


## 3. Adjoint Problem

Since we adopt the Jacobian-based minimization method to solve (10) in this paper, the derivative of the forward operator becomes a matrix which is called Jacobi matrix. For convenience of representing the element of Jacobi matrix hereinafter, we use [Jμai] to represent the Jacobi matrix of Fi, with respect to μa=(μa,1,…,μa,N), 1≤i≤s; here μa,i denotes the ith element of the discretized μa. For each i, [Jμai]∈Rd×N is defined as follows:(13)
J
μ
a
i
=
∂
M
1
,
i
∂
μ
a
,
1
∂
M
1
,
i
∂
μ
a
,
2
⋯
∂
M
1
,
i
∂
μ
a
,
N
∂
M
2
,
i
∂
μ
a
,
1
∂
M
2
,
i
∂
μ
a
,
2
⋯
∂
M
2
,
i
∂
μ
a
,
N
⋯
⋯
⋮
⋯
∂
M
d
,
i
∂
μ
a
,
1
∂
M
d
,
i
∂
μ
a
,
2
⋯
∂
M
d
,
i
∂
μ
a
,
N
,
1
≤
i
≤
s
.[Jμai] can be computed by direct method which differentiates Mζi(ξj) with respect to the perturbation of μa on each element. In order to compute [Jμai], for each i, we need to solve forward problems d×N times. For a total of s sources, we need to compute s×d×N times, so the direct method is very time consuming. Therefore, we use the adjoint formulation to compute [Jμai] instead of direct method.We first consider the analytic Fréchet derivative of the forward operator for the boundary value problem of RTE given by [25].(14)Fi′μa∗Mζix=-∫Ωuix,ωφx,ωdσω,where Fi′(μa)∗ denotes the adjoint operator of Fi(μa) and φ(x,ω) is the solution of the adjoint RTE(15)-ω·∇φx,ω+μa+μsφx,ω=μs∫Ωkω,ω^φx,ω^dσω^inX×Ωwith boundary condition(16)φx,ω=ω·νxMix,x∈∂X,ω·νx>0,0,otherwise,on∂X×Ω.From the adjoint RTE (15) and the boundary condition (16), it seems that the adjoint problem should be solved in the reverse direction of propagation with a completely new computation solver. However, by the standard reciprocity theorem for the Boltzmann equation given in [31],(17)Gx,ω;x0,ω0=Gx0,-ω0;x,-ω,which states that the angular density at x in direction ω due to a source at x0 in direction ω0 is the same as the angular density at x0 in direction -ω0 due to a source at x in direction -ω. Here, G(·;·) presents the Green function of (1) for isotropic point source on the boundary. Then the adjoint problems (15) and (16) can be transformed into the same form as the forward problem, simply replacing direction ω with -ω. Then we just need to solve the following equation for the radiance φ^(x,-ω)=φ(x,ω) by the same forward solver:(18)ω·∇φ^x,ω+μa+μsφ^x,ω=μs∫Ωkω,ω^φ^x,ω^dσω^,with adjoint boundary condition(19)φ^x,ω=ω·νxMix,x∈∂X,ω·νx<0,0,otherwise.This means we need to solve (18) and (19) with the forward solver firstly and then reverse all directions on solution.If we consider one source positionζi and one detector position ξj, then the kth column of the Jacobi matrix [Jμai] can be computed as follow:(20)Jμai:,k=∫Xχkx∫Ωuix,ωφjx,ωdσωdx,where φj is the solution of the following boundary value problem:(21)-ω·∇φjx,ω+μax+μsxφjx,ω=μsx∫Ωkω·ω^φjx,ω^dσω^,φjx,ω=ω·νxMixx=ξj,ω·νx>0,0,otherwise.Thus for the adjoint method to compute [Jμai], for each 1≤i≤s, we only need to solve forward problems d times. For s sources, we need to solve forward problems s×d times, which dramatically reduces the computation burden and improves the efficiency of the algorithm.


## 4. Iterative Procedure

The unknownμa can be estimated through the regularization method. The quality of reconstructed image strongly depends on the choice of the penalty term R(μa). If we choose the L2 norm as penalty, that is, R(μa)=μa22, the reconstructed image usually blurs with a low resolution. If we choose L1 norm as the penalty, that is, R(μa)=μa1, the reconstructed image tends to find a sparse solution [10, 18]. To treat the discontinuity and the edges of different distribution regions, total variation is usually chosen as the penalty functional, that is,(22)Rμa=∫X∇μa.The symbol ∫X|∇μa| denotes the total variation seminorm [32] of μa∈L1(X) as follows:(23)∫X∇μa≔sup∫Xμadivφdx∣φ∈C0∞X;Rd,φL∞X;Rd≤1.To improve the reconstruction quality, we consider the total variation mixed with the L1 norm as the penalty. The functional to be minimized is of the form(24)Jμa=12∑i=1sFiμa-MiL2∂X2+α∫X∇μa+βμal1,where α,β are the regularization parameters. We will apply the split Bregman method to solve (24) [33]. Instead of (24), we consider the following constrained optimization problem:(25)infμa,D⁡12∑i=1sFiμa-MiL2∂X2+α∫X∇μa+βDl1such  thatD=μa.Solution of the above minimization problem can be obtained by solving the unconstrained optimization problem(26)μa,D=arg⁡min⁡12∑i=1sFiμa-MiL2∂X2+α∫X∇μa+βDl1+η2D-μal22,where η>0 is the split parameter. Now let us iteratively solve the following subproblems [34]:(27)μan+1,Dn+1=arg⁡min⁡12∑i=1sFiμa-MiL2∂X2+α∫X∇μa+βDl1+η2D-μa-bdn22,with the following update for bd: (28)bdn+1=bdn+μan+1-Dn+1.The minimization of subproblems in (1) can be iteratively solved by splitting it into the minimizations of μa and D separately. This suggests the following steps.Step 1.
(29) μ a n + 1 = arg ⁡ min μ a ⁡ 1 2 ∑ i = 1 s F i μ a - M i L 2 ∂ X 2 + α ∫ X ∇ μ a + η 2 D n - μ a - b d n 2 2 .Step 2.
(30) D n + 1 = arg ⁡ min D ⁡ β D l 1 + η 2 D - μ a n + 1 - b d n 2 2 .Step 3.
(31) b d n + 1 = b d n + μ a n + 1 - D n + 1 .For the solution of Step1, we use the Levenberg-Marquardt method, which has a high convergence rate. Thus we solve a minimization problem as follows.Step 1 ∗ (32) μ a n + 1 = arg ⁡ min μ a ⁡ 1 2 ∑ i = 1 s F i μ a n + J μ a n i μ a - μ a n - M i 2 2 + α ∫ X ∇ μ a + η 2 D n - μ a - b d n 2 2 ,where [Jμani] denotes the Jacobi matrix of the forward operator Fi with respect to μan. For solving Step  1∗, in order to avoid the nondifferentiability of the total variation term ∫X|∇μa| at a zero point, we approximate it with a smooth functional μaTV,ε defined as(33)μaTV,ε=∫X∇μa2+ϵ2dx,ϵ>0.The parameter ϵ is a constant, and it cannot be too large or too small. The Euler-Lagrange equation for Step1∗ is(34)∑i=1sJμaniTJμaniμa-μan-Riμan+αLμanμa+ηDn-μa-bdn=0,where ·T denotes the transpose, Ri(μan) is the short notion of Fi(μan)-Mi, the term L(μan)μa is defined as(35)Lμanμa=-∇·∇μa∇μan2+ϵ2.Then we can get the iterative update by solving the above equation with Newton method as follows:(36)∑i=1sJμaniTJμani+αLμan+ηIμa-μan=∑i=1sJμaniTRiμan+αLμanμan+ηDn-μan-bdn.Step 2 is an L1 norm regularization problem and it can be solved efficiently through the shrinkage operator, that is,(37)Dn+1=shrinkμan+1+bdn,βη,where the shrinkage operator(38)shrinkx,t=sin⁡xmax⁡x-t,0=x-t,x≥t,0,x<t,x+t,x≤-t.Implementation of the split Bregman formulation for our TV-L1 regularization is described in Algorithm 1.Algorithm 1:Reconstruction algorithm based on the split Bregman method for TV-L1 regularization.
Input: Given the initial guess μa0; the regularization parameter α>0, β>0, η>0; D0=bd0=0, and the tolerance ε. whileμan+1-μan>ε do For1≤i≤s, compute μan+1 with (36); ComputeDn+1=shrinkμan+1+bdn,β/η; Computebdn+1=bdn+μan+1-Dn+1. end while Output: Output an approximation μa=μan+1.


## 5. Reconstruction Results and Discussion

In this section, theTV-L1 regularization with the split Bregman formulation is applied to 2D test problems. We assume the distribution of the scattering coefficient is known. Reconstructions of spatially dependent distributions of absorption coefficient inside the medium are performed and discussed. In our simulation, a circular domain X which contains different inclusions is investigated. The radius of the circle is 10 mm. 12 sources and 12 detectors are located on the boundary of the domain with equal space. This yields totally 144 source-detector pairs to be used in the inversion.Noise-free synthetic data are generated by solving the forward problem on triangular meshes with the method we mentioned in Section2.1. Note that the meshes for the inverse problem are coarser than the meshes for the forward problem in order to avoid the “inverse crime” [35] and for the purpose of testing the robustness of the proposed algorithm. We will describe the meshes in each case. In all the cases, the angular space is discretized into 32 directions which equally divide the interval [0,2π].The background medium has scattering coefficient of 10 mm-1 and absorption coefficient of 0.01 mm-1, respectively, and these values keep the same throughout this paper. The initial guess of reconstruction is set to be identical to the properties of the background medium. That is, the iterative procedure started with the background value of the absorption coefficient.The incident impulse on the inflow boundary is settled as(39)uin,ix,ω=Bi,ωi,x,ω∈∂X×Ω,ω·νx<0,1≤i≤s,where Bi,ωi is a piecewise linear function whose spatial support is Si and achieves the value 1 at the center node of Si, where Si is the element through which theith incident impulse passes. The direction of uin,i points approximately from the center of Si to the center of X.Based on above settings, we use various simulations to validate the proposed split Bregman algorithm for ourTV-L1 regularization. Our purpose is to show the following results.First, by comparing the reconstruction with different anisotropic factorsg, the proposed algorithm works better when g is bigger.Second, for small sparse inclusions, theTV-L1 regularization can reconstruct the location and the quality of the coefficient more accurately than the TV regularization and the L1 regularization.Third, the proposed split algorithm forTV-L1 needs less computation time than the Levenberg-Marquardt algorithm for TV regularization. Even though our proposed algorithm needs more computation time than the split Bregman algorithm for L1 regularization [10], the reconstruction quality with the proposed algorithm is much better than that with the split Bregman algorithm for L1 regularization.


### 5.1. Simulation 1: Reconstruction with Different Anisotropic Factors

In this simulation, we reconstruct one circle inclusion with the radius 0.5 mm centered at (7.0 mm, 0.0 mm). In this circle, the absorption coefficientμa=0.02mm-1 and μs=20mm-1. Our goal is to reconstruct μa. For solving the inverse problem, we use a mesh of 772 nodes and 1484 triangular elements. To avoid the inverse crime, we use a mesh of 1296 nodes and 2440 triangular elements for the synthetic data. The simulated absorption distributions and the meshes can be seen in Figure 1. We solve the minimizer of J(μa) by using the split Bregman method for TV-L1 method with exact data, that is, δ=0, in this example. Then we compare the reconstruction results for different anisotropic factors g. For the comparison purpose, we use the same parameters α, β, η, and ε for different g. Noticing that the so-called exact data in fact contains noise, so the regularization parameter α,β should not be too small. Since the inclusion is very small in this example, we enhance more weight of L1 penalty than the weight of TV penalty. Here we take α=10-4 and β=10-3. The other two parameters are taken as η=10-5 and ϵ=10-6, respectively.Figure 1
The forward mesh, inverse mesh, and the true distribution ofμa in simulation 1.
(a)
Forward mesh (b)
Inverse mesh (c)
TrueμaThe reconstruction results can be seen in Figure2. The first row of Figure 2 is the reconstruction results for g=0.1 and g=0.4 from the left to the right. The second row of Figure 2 is the reconstruction results for g=0.7 and g=0.9 from the left to the right. As it can be seen, the reconstructions for the bigger g are more clear than that for smaller g. There are some apparent blurred parts in the left side of the inclusion in Figures 2(a) and 2(b). The blurred parts are much smaller in Figure 2(c) than in Figure 2(b) and almost disappeared in Figure 2(d).Figure 2
Reconstruction results with and without scaling strategy. (a) Reconstruction ofμa for g=0.1; (b) reconstruction of μa for g=0.4; (c) reconstruction of μa for g=0.7; (d) reconstruction of μa for g=0.9. The reconstructions are done by using the split Bregman algorithm for TV-L1 regularization with noise-free synthetic data.
(a) (b) (c) (d)From this example, we can find that the proposed algorithm works better for the bigger anisotropic factorg. The following simulations in this paper are all performed with g=0.9.


### 5.2. Simulation 2: Reconstruction with TV-L1 Regularization and the Standard TV Regularization

There are two groups of experiments in this section. In the first group of experiment, we design three types of inclusions to compare the performance of theTV-L1 regularization and the TV regularization with noise-free synthetic data. In the first case, a small inclusion with radius 0.5 mm centered at (7.0 mm, 0.0 mm) is designed for μa=0.02mm-1 and μs=20.0mm-1, respectively. In the second case, a middle inclusion with radius 2 mm centered at (5.0 mm, 0.0 mm) is designed for μa=0.02mm-1 and μs=20mm-1, respectively. In the third case, a middle inclusion with radius 4 mm centered at (3.0 mm, 0.0 mm) is designed for μa=0.02mm-1 and μs=20mm-1, respectively. The simulated true absorption coefficient can be seen on the first column of Figure 3. For the small inclusion case, we use a mesh of 772 nodes and 1484 triangular elements for the inverse problem and a mesh of 1267 nodes and 2440 triangular elements for the forward problem, the simulated true absorption coefficient and reconstruction results with TV-L1 regularization and TV regularization are shown in Figure 3. For the middle inclusion case, we use a mesh of 575 nodes and 1080 triangular elements for the inverse problem and a mesh of 813 nodes and 2642 triangular elements for the forward problem, the simulated true absorption coefficient and reconstruction results with TV-L1 regularization and TV regularization are shown in Figure 4. For the big inclusion case, we use a mesh of 583 nodes and 1096 triangular elements for the inverse problem and a mesh of 871 nodes and 3124 triangular elements for the forward problem, the simulated true absorption coefficient and reconstruction results with TV-L1 regularization and TV regularization are shown in Figure 5. The values of the parameters are shown in Table 1. Since we use the same ε=10-6 as the last section, we will not present it repeatedly in Table 1.Table 1
Parameters corresponding to various regularization.
Cases Small inclusion Middle inclusion Big inclusion ( α , β , η ) for TV-L1 ( 0.0003,0.003,0.0003 ) ( 0.003,0.003,0.0003 ) ( 0.003,0.0003,0.0003 ) ( α , β , η ) for TV ( 3 × 1 0 - 4 , 0,0 ) ( 3 × 1 0 - 4 , 0,0 ) ( 3 × 1 0 - 4 , 0,0 )Figure 3
Reconstruction results with TV-L1 regularization and TV regularization for one small inclusion. (a) True μa; (b) reconstruction of μa with TV-L1 penalty; (c) reconstruction of μa with only TV penalty.
(a) (b) (c)Figure 4
Reconstruction results with TV-L1 regularization and TV regularization for one middle inclusion. (a) True μa; (b) reconstruction of μa with TV-L1 penalty; (c) reconstruction of μa with only TV penalty.
(a) (b) (c)Figure 5
Reconstruction results with TV-L1 regularization and TV regularization for one big inclusion. (a) True μa; (b) reconstruction of μa with TV-L1 penalty; (c) reconstruction of μa with only TV penalty.
(a) (b) (c)Observing from the results, we can see that, for small and middle inclusions, theTV-L1 regularization performs better than the TV regularization in both localizing the location and quantifying the values; see Figures 3(b) and 4(b) for TV-L1 regularization and Figures 3(c) and 4(c) for TV regularization. In fact, these results are reasonable and can be interpreted. The TV-L1 regularization imposes both TV penalty and L1 penalty on μa. The TV penalty tends to find edges and the L1 penalty tends to find the sparse details of the inclusion. As can be seen from Table 1, for the big inclusion, we enhance the weight of TV penalty. When the inclusion gets smaller, we enhance the weight of L1 penalty.For the big inclusion,TV-L1 regularization and TV regularization perform no big differences, but we still can see that the blurred parts in Figure 5(b) with TV-L1 regularization are smaller than that in Figure 5(c) with TV regularization. We can see from Figures 5(b) and 5(c) that there are large area of blurs in the big inclusion case. It is reasonable in the sense that large area of absorption inclusion means more photons are absorbed in propagation process. Hence, we can alleviate this phenomenon by increasing the source-detector pairs. In other words, by increasing the measurements, one can alleviate the effect of absorption to some extent.In the second group of experiments of this section, three types of complicated and multi-inclusions are designed for comparing the reconstruction results byTV-L1 regularization, TV regularization, and L1 regularization.In the first case, we reconstruct two small circle inclusions centered at (7 mm, 0 mm) and (0 mm, 0 mm) with the same radius 0.5 mm. In both of the two inclusions,μa=0.02mm-1 and μs=20mm-1. We use a mesh of 1277 nodes and 2488 triangular elements for the inverse problem and a mesh of 1861 nodes and 3616 triangular elements for the forward problem. The true distributions of μa and reconstruction results with various regularizations can be seen in Figure 6.Figure 6
Reconstruction results with various regularization for two small inclusions. (a) Trueμa; (b) reconstruction of μa with TV-L1 penalty; (c) reconstruction of μa with TV penalty; (d) reconstruction of μa with L1 penalty.
(a) (b) (c) (d)In the second case, we reconstruct three small circles centered at (7 mm, 0 mm), (0 mm, −7 mm), and (0 mm, 0 mm) with the same radius, 0.5 mm. In all the three inclusions,μa=0.02mm-1 and μs=20mm-1. We use a mesh of 1671 nodes and 3264 triangular elements for the inverse problem and a mesh of 2141 nodes and 4176 triangulations for the forward problem. The true distributions of μa and reconstruction results with various regularizations can be seen in Figure 7.Figure 7
Reconstruction results with various regularization for three small inclusions. (a) Trueμa; (b) reconstruction of μa with TV-L1 penalty; (c) reconstruction of μa with TV penalty; (d) reconstruction of μa with L1 penalty.
(a) (b) (c) (d)In the third case, we design three inclusions. One big circle centered at (−5 mm, 2 mm) with radius 2 mm. Two small circles centered at (0 mm, −7 mm) and (7 mm, 0 mm) with the same radius 0.5 mm. In all the three cases, we take the same coefficient valuesμa=0.02mm-1 and μs=20mm-1. We use a mesh of 1323 nodes and 2568 triangulations for the inverse problem and a mesh of 1809 nodes and 3512 triangulations for the forward problem. The true distributions of μa and reconstruction results with various regularizations can be seen in Figure 8.Figure 8
Reconstruction results with various regularization for one big and two small inclusions. (a) Trueμa; (b) reconstruction of μa with TV-L1 penalty; (c) reconstruction of μa with TV penalty; (d) reconstruction of μa with L1 penalty.
(a) (b) (c) (d)In Table2, for all the three cases, we present the parameter settings for the iterative procedure (α, β, η), the steps of the iterations Iters, the relative residual error Eresi, and the relative solution error Eμa in the sense of L2 norm corresponding to each regularization method. Eresi and Eμa are defined as(40)Eresi=Fμan-FμatrueFμatrue,Eμa=μan-μatrueμatrue.Table 2
Regularization parameters, error estimates, and iteration number for different cases in the second example in Section5.2.
Cases Parameters TV-L1 TV L 1 Two small inclusions ( α , β , η ) ( 1 0 - 4 , 1 0 - 3 , 1 0 - 5 ) ( 1 0 - 4 , 0,0 ) ( 0,1 0 - 3 , 1 0 - 5 ) Iters 11 39 5 E r e s i 7.3 × 1 0 - 5 6.4 × 1 0 - 5 3.3 × 1 0 - 4 E μ a 13.29% 14.78% 15.76%  Three small inclusions ( α , β , η ) ( 1 0 - 4 , 1 0 - 3 , 1 0 - 5 ) ( 1 0 - 3 , 0,0 ) ( 0,1 0 - 3 , 1 0 - 5 ) Iters 8 79 6 E r e s i 8.6 × 1 0 - 5 7.0 × 1 0 - 5 6.3 × 1 0 - 4 E μ a 15.02% 15.80% 15.88%  One big and two small inclusions ( α , β , η ) ( 1 0 - 3 , 1 0 - 3 , 1 0 - 5 ) ( 1 0 - 3 , 0,0 ) ( 0,1 0 - 4 , 1 0 - 6 ) Iters 32 53 6 E r e s i 9.5 × 1 0 - 4 6.8 × 1 0 - 4 0.0019 E μ a 15.75% 15.47% 15.17%We say that in the first and the second case, compared with the inclusions near the boundary, the inclusion deep in the domain is difficult to identify. In the third case, it is more difficult to reconstruct the small inclusion than the big inclusion. From the reconstruction results, we can see that theTV-L1 regularization can preferably identify the location and the value of the inclusions. In other words, the advantages with TV-L1 regularization are indeed apparent over TV regularization and L1 regularization; see Figures 6(b), 7(b), and 8(b) for the reconstructions results with TV-L1 regularization. The LM type algorithm for the TV regularization converges slowly especially when the inclusions are very small; see Figures 6(c) and 7(c). When big inclusion is included, the LM for TV regularization can well identify it, but it still cannot perfectly reconstruct the value of the small inclusions; see Figure 8(c). The split Bregman algorithm for L1 regularization is an efficient algorithm with high convergent rate; we can see this from the number of iterations in Table 2. But, in our simulations, although the algorithm converges quickly, it loses validity when the inclusions are very small; see Figures 6(d) and 7(d). In Figure 8(d), the big inclusion is identified, but we only can see two blurred circles at where the two small inclusions are located. Moreover, we can find that the reconstruction may be oversparsified by using the L1 regularization.


### 5.3. Simulation 3: High Absorbing and Low Scattering Inclusions with TV-L1

As the last simulation, our purpose is to show the validity of our algorithm in the situation that high absorbing inclusion and low scattering inclusion are contained in the domainX. Moreover, we investigated the reconstruction under different noise level to show the robustness of the proposed algorithm. The noises added to the exact data Mitrue are by the following rule: Mi=Mitrue(1+δ∗random), where δ is the signal-to-noise ratio; random is a Gaussian random variable with zero mean and unity variation. The background values of the absorption and the scattering are 0.05mm-1 and 5mm-1, respectively. The absorbing inclusion for which we set the absorption coefficient is μa=0.1mm-1 which is centered at (3.5 mm, 3.5 mm) with radius 2 mm. The scattering inclusion for which we set the scattering coefficient is μs=10mm-1 which is centered at (0 mm, −5 mm) with radius 2 mm. For the discretization, we use a mesh of 463 nodes and 856 triangulations for the inverse problem and a mesh of 681 nodes and 1288 triangulations for the forward problem. See the first row of the Figure 9(a) for the true distribution of μa. In Figure 9(b), the reconstruction of μa with exact data is presented. The parameter value in this case is chosen as α=0.0005, β=0.0005, and η=10-6. Figure 9(c) presents the reconstruction with noise level 0.1%, in which the parameter is taken as α=0.0005, β=0.0005, and η=10-6. In Figure 9(d), the reconstruction with noise level 1% is shown, in which the parameter value is taken as α=0.005, β=0.0005, and η=10-6. From the four reconstructions under four different noise levels, we can find that as the noise level increased, the quality of the reconstruction gets worse. Here, we take the noise level under 1%, since the iteration will not converge if the noise level is bigger than 1% which reflects the severe ill-posedness of the problem. We also find that in order to obtain the reconstruction results with the relative error in L2 norm no more than 20% the noise level on the synthetic data should be less than 1%.Figure 9
Reconstruction ofμa with Algorithm 1 for high absorption and low scattering situation. (a) True μa; (b) reconstruction of μa with exact data; (c) reconstruction of μa with noise level 0.1%; (d) reconstruction of μa with noise level 1%.
(a) (b) (c) (d)


## 5.1. Simulation 1: Reconstruction with Different Anisotropic Factors

In this simulation, we reconstruct one circle inclusion with the radius 0.5 mm centered at (7.0 mm, 0.0 mm). In this circle, the absorption coefficientμa=0.02mm-1 and μs=20mm-1. Our goal is to reconstruct μa. For solving the inverse problem, we use a mesh of 772 nodes and 1484 triangular elements. To avoid the inverse crime, we use a mesh of 1296 nodes and 2440 triangular elements for the synthetic data. The simulated absorption distributions and the meshes can be seen in Figure 1. We solve the minimizer of J(μa) by using the split Bregman method for TV-L1 method with exact data, that is, δ=0, in this example. Then we compare the reconstruction results for different anisotropic factors g. For the comparison purpose, we use the same parameters α, β, η, and ε for different g. Noticing that the so-called exact data in fact contains noise, so the regularization parameter α,β should not be too small. Since the inclusion is very small in this example, we enhance more weight of L1 penalty than the weight of TV penalty. Here we take α=10-4 and β=10-3. The other two parameters are taken as η=10-5 and ϵ=10-6, respectively.Figure 1
The forward mesh, inverse mesh, and the true distribution ofμa in simulation 1.
(a)
Forward mesh (b)
Inverse mesh (c)
TrueμaThe reconstruction results can be seen in Figure2. The first row of Figure 2 is the reconstruction results for g=0.1 and g=0.4 from the left to the right. The second row of Figure 2 is the reconstruction results for g=0.7 and g=0.9 from the left to the right. As it can be seen, the reconstructions for the bigger g are more clear than that for smaller g. There are some apparent blurred parts in the left side of the inclusion in Figures 2(a) and 2(b). The blurred parts are much smaller in Figure 2(c) than in Figure 2(b) and almost disappeared in Figure 2(d).Figure 2
Reconstruction results with and without scaling strategy. (a) Reconstruction ofμa for g=0.1; (b) reconstruction of μa for g=0.4; (c) reconstruction of μa for g=0.7; (d) reconstruction of μa for g=0.9. The reconstructions are done by using the split Bregman algorithm for TV-L1 regularization with noise-free synthetic data.
(a) (b) (c) (d)From this example, we can find that the proposed algorithm works better for the bigger anisotropic factorg. The following simulations in this paper are all performed with g=0.9.


## 5.2. Simulation 2: Reconstruction with TV-L1 Regularization and the Standard TV Regularization

There are two groups of experiments in this section. In the first group of experiment, we design three types of inclusions to compare the performance of theTV-L1 regularization and the TV regularization with noise-free synthetic data. In the first case, a small inclusion with radius 0.5 mm centered at (7.0 mm, 0.0 mm) is designed for μa=0.02mm-1 and μs=20.0mm-1, respectively. In the second case, a middle inclusion with radius 2 mm centered at (5.0 mm, 0.0 mm) is designed for μa=0.02mm-1 and μs=20mm-1, respectively. In the third case, a middle inclusion with radius 4 mm centered at (3.0 mm, 0.0 mm) is designed for μa=0.02mm-1 and μs=20mm-1, respectively. The simulated true absorption coefficient can be seen on the first column of Figure 3. For the small inclusion case, we use a mesh of 772 nodes and 1484 triangular elements for the inverse problem and a mesh of 1267 nodes and 2440 triangular elements for the forward problem, the simulated true absorption coefficient and reconstruction results with TV-L1 regularization and TV regularization are shown in Figure 3. For the middle inclusion case, we use a mesh of 575 nodes and 1080 triangular elements for the inverse problem and a mesh of 813 nodes and 2642 triangular elements for the forward problem, the simulated true absorption coefficient and reconstruction results with TV-L1 regularization and TV regularization are shown in Figure 4. For the big inclusion case, we use a mesh of 583 nodes and 1096 triangular elements for the inverse problem and a mesh of 871 nodes and 3124 triangular elements for the forward problem, the simulated true absorption coefficient and reconstruction results with TV-L1 regularization and TV regularization are shown in Figure 5. The values of the parameters are shown in Table 1. Since we use the same ε=10-6 as the last section, we will not present it repeatedly in Table 1.Table 1
Parameters corresponding to various regularization.
Cases Small inclusion Middle inclusion Big inclusion ( α , β , η ) for TV-L1 ( 0.0003,0.003,0.0003 ) ( 0.003,0.003,0.0003 ) ( 0.003,0.0003,0.0003 ) ( α , β , η ) for TV ( 3 × 1 0 - 4 , 0,0 ) ( 3 × 1 0 - 4 , 0,0 ) ( 3 × 1 0 - 4 , 0,0 )Figure 3
Reconstruction results with TV-L1 regularization and TV regularization for one small inclusion. (a) True μa; (b) reconstruction of μa with TV-L1 penalty; (c) reconstruction of μa with only TV penalty.
(a) (b) (c)Figure 4
Reconstruction results with TV-L1 regularization and TV regularization for one middle inclusion. (a) True μa; (b) reconstruction of μa with TV-L1 penalty; (c) reconstruction of μa with only TV penalty.
(a) (b) (c)Figure 5
Reconstruction results with TV-L1 regularization and TV regularization for one big inclusion. (a) True μa; (b) reconstruction of μa with TV-L1 penalty; (c) reconstruction of μa with only TV penalty.
(a) (b) (c)Observing from the results, we can see that, for small and middle inclusions, theTV-L1 regularization performs better than the TV regularization in both localizing the location and quantifying the values; see Figures 3(b) and 4(b) for TV-L1 regularization and Figures 3(c) and 4(c) for TV regularization. In fact, these results are reasonable and can be interpreted. The TV-L1 regularization imposes both TV penalty and L1 penalty on μa. The TV penalty tends to find edges and the L1 penalty tends to find the sparse details of the inclusion. As can be seen from Table 1, for the big inclusion, we enhance the weight of TV penalty. When the inclusion gets smaller, we enhance the weight of L1 penalty.For the big inclusion,TV-L1 regularization and TV regularization perform no big differences, but we still can see that the blurred parts in Figure 5(b) with TV-L1 regularization are smaller than that in Figure 5(c) with TV regularization. We can see from Figures 5(b) and 5(c) that there are large area of blurs in the big inclusion case. It is reasonable in the sense that large area of absorption inclusion means more photons are absorbed in propagation process. Hence, we can alleviate this phenomenon by increasing the source-detector pairs. In other words, by increasing the measurements, one can alleviate the effect of absorption to some extent.In the second group of experiments of this section, three types of complicated and multi-inclusions are designed for comparing the reconstruction results byTV-L1 regularization, TV regularization, and L1 regularization.In the first case, we reconstruct two small circle inclusions centered at (7 mm, 0 mm) and (0 mm, 0 mm) with the same radius 0.5 mm. In both of the two inclusions,μa=0.02mm-1 and μs=20mm-1. We use a mesh of 1277 nodes and 2488 triangular elements for the inverse problem and a mesh of 1861 nodes and 3616 triangular elements for the forward problem. The true distributions of μa and reconstruction results with various regularizations can be seen in Figure 6.Figure 6
Reconstruction results with various regularization for two small inclusions. (a) Trueμa; (b) reconstruction of μa with TV-L1 penalty; (c) reconstruction of μa with TV penalty; (d) reconstruction of μa with L1 penalty.
(a) (b) (c) (d)In the second case, we reconstruct three small circles centered at (7 mm, 0 mm), (0 mm, −7 mm), and (0 mm, 0 mm) with the same radius, 0.5 mm. In all the three inclusions,μa=0.02mm-1 and μs=20mm-1. We use a mesh of 1671 nodes and 3264 triangular elements for the inverse problem and a mesh of 2141 nodes and 4176 triangulations for the forward problem. The true distributions of μa and reconstruction results with various regularizations can be seen in Figure 7.Figure 7
Reconstruction results with various regularization for three small inclusions. (a) Trueμa; (b) reconstruction of μa with TV-L1 penalty; (c) reconstruction of μa with TV penalty; (d) reconstruction of μa with L1 penalty.
(a) (b) (c) (d)In the third case, we design three inclusions. One big circle centered at (−5 mm, 2 mm) with radius 2 mm. Two small circles centered at (0 mm, −7 mm) and (7 mm, 0 mm) with the same radius 0.5 mm. In all the three cases, we take the same coefficient valuesμa=0.02mm-1 and μs=20mm-1. We use a mesh of 1323 nodes and 2568 triangulations for the inverse problem and a mesh of 1809 nodes and 3512 triangulations for the forward problem. The true distributions of μa and reconstruction results with various regularizations can be seen in Figure 8.Figure 8
Reconstruction results with various regularization for one big and two small inclusions. (a) Trueμa; (b) reconstruction of μa with TV-L1 penalty; (c) reconstruction of μa with TV penalty; (d) reconstruction of μa with L1 penalty.
(a) (b) (c) (d)In Table2, for all the three cases, we present the parameter settings for the iterative procedure (α, β, η), the steps of the iterations Iters, the relative residual error Eresi, and the relative solution error Eμa in the sense of L2 norm corresponding to each regularization method. Eresi and Eμa are defined as(40)Eresi=Fμan-FμatrueFμatrue,Eμa=μan-μatrueμatrue.Table 2
Regularization parameters, error estimates, and iteration number for different cases in the second example in Section5.2.
Cases Parameters TV-L1 TV L 1 Two small inclusions ( α , β , η ) ( 1 0 - 4 , 1 0 - 3 , 1 0 - 5 ) ( 1 0 - 4 , 0,0 ) ( 0,1 0 - 3 , 1 0 - 5 ) Iters 11 39 5 E r e s i 7.3 × 1 0 - 5 6.4 × 1 0 - 5 3.3 × 1 0 - 4 E μ a 13.29% 14.78% 15.76%  Three small inclusions ( α , β , η ) ( 1 0 - 4 , 1 0 - 3 , 1 0 - 5 ) ( 1 0 - 3 , 0,0 ) ( 0,1 0 - 3 , 1 0 - 5 ) Iters 8 79 6 E r e s i 8.6 × 1 0 - 5 7.0 × 1 0 - 5 6.3 × 1 0 - 4 E μ a 15.02% 15.80% 15.88%  One big and two small inclusions ( α , β , η ) ( 1 0 - 3 , 1 0 - 3 , 1 0 - 5 ) ( 1 0 - 3 , 0,0 ) ( 0,1 0 - 4 , 1 0 - 6 ) Iters 32 53 6 E r e s i 9.5 × 1 0 - 4 6.8 × 1 0 - 4 0.0019 E μ a 15.75% 15.47% 15.17%We say that in the first and the second case, compared with the inclusions near the boundary, the inclusion deep in the domain is difficult to identify. In the third case, it is more difficult to reconstruct the small inclusion than the big inclusion. From the reconstruction results, we can see that theTV-L1 regularization can preferably identify the location and the value of the inclusions. In other words, the advantages with TV-L1 regularization are indeed apparent over TV regularization and L1 regularization; see Figures 6(b), 7(b), and 8(b) for the reconstructions results with TV-L1 regularization. The LM type algorithm for the TV regularization converges slowly especially when the inclusions are very small; see Figures 6(c) and 7(c). When big inclusion is included, the LM for TV regularization can well identify it, but it still cannot perfectly reconstruct the value of the small inclusions; see Figure 8(c). The split Bregman algorithm for L1 regularization is an efficient algorithm with high convergent rate; we can see this from the number of iterations in Table 2. But, in our simulations, although the algorithm converges quickly, it loses validity when the inclusions are very small; see Figures 6(d) and 7(d). In Figure 8(d), the big inclusion is identified, but we only can see two blurred circles at where the two small inclusions are located. Moreover, we can find that the reconstruction may be oversparsified by using the L1 regularization.


## 5.3. Simulation 3: High Absorbing and Low Scattering Inclusions with TV-L1

As the last simulation, our purpose is to show the validity of our algorithm in the situation that high absorbing inclusion and low scattering inclusion are contained in the domainX. Moreover, we investigated the reconstruction under different noise level to show the robustness of the proposed algorithm. The noises added to the exact data Mitrue are by the following rule: Mi=Mitrue(1+δ∗random), where δ is the signal-to-noise ratio; random is a Gaussian random variable with zero mean and unity variation. The background values of the absorption and the scattering are 0.05mm-1 and 5mm-1, respectively. The absorbing inclusion for which we set the absorption coefficient is μa=0.1mm-1 which is centered at (3.5 mm, 3.5 mm) with radius 2 mm. The scattering inclusion for which we set the scattering coefficient is μs=10mm-1 which is centered at (0 mm, −5 mm) with radius 2 mm. For the discretization, we use a mesh of 463 nodes and 856 triangulations for the inverse problem and a mesh of 681 nodes and 1288 triangulations for the forward problem. See the first row of the Figure 9(a) for the true distribution of μa. In Figure 9(b), the reconstruction of μa with exact data is presented. The parameter value in this case is chosen as α=0.0005, β=0.0005, and η=10-6. Figure 9(c) presents the reconstruction with noise level 0.1%, in which the parameter is taken as α=0.0005, β=0.0005, and η=10-6. In Figure 9(d), the reconstruction with noise level 1% is shown, in which the parameter value is taken as α=0.005, β=0.0005, and η=10-6. From the four reconstructions under four different noise levels, we can find that as the noise level increased, the quality of the reconstruction gets worse. Here, we take the noise level under 1%, since the iteration will not converge if the noise level is bigger than 1% which reflects the severe ill-posedness of the problem. We also find that in order to obtain the reconstruction results with the relative error in L2 norm no more than 20% the noise level on the synthetic data should be less than 1%.Figure 9
Reconstruction ofμa with Algorithm 1 for high absorption and low scattering situation. (a) True μa; (b) reconstruction of μa with exact data; (c) reconstruction of μa with noise level 0.1%; (d) reconstruction of μa with noise level 1%.
(a) (b) (c) (d)


## 6. Summary

In this paper, forward and inverse problems of the radiative transfer equation are considered. An image reconstruction method based on theTV-L1 regularization is proposed. The forward problem is discretized with the discontinuous Galerkin method on the spatial space and the finite element method on the angular space which both are implemented on the piecewise linear basis. We discretize the absorption and scattering coefficients on the piecewise constant basis. The minimization problem is solved by a Levenberg-Marquardt type method which is equipped with a split Bregman algorithm for the L1 penalty. The adjoint method is used to compute the Jacobian matrix which is the discretized Fréchet derivative of the forward operator. We numerically compare the proposed reconstruction method with the other imaging reconstruction methods based on TV and L1 regularizations. The simulation results show the validity and robustness of the proposed method.



---

*Source: 2953560-2017-02-09.xml*