Multilevel bioluminescence tomography&#x2028;based on radiative transfer equation&#x2028;Part 2: total variation and l1 data fidelity

Hao Gao; Hongkai Zhao

doi:10.1364/OE.18.002894

1. Introduction

In bioluminescence tomography (BLT) [1–3] with proper discretization, boundary measurements $Χ^{m} : = {x_{j}^{m}, 1 \leq j \leq M}$ are linearly dependent on the bioluminescent sources $Q : = {q_{i}, 1 \leq i \leq N_{s}}$ . For conventional reconstruction methods for BLT, people use diffusion approximation (DA) [4] as forward model, and then formulate it as a standard least-square problem, i.e.,

\min_{Q} | | J Q - X^{m} | |_{2}^{2} + \frac{λ_{2}}{2} | | Q | |_{2}^{2},

in whichJis the system matrix usually called Jacobian or sensitivity matrix corresponding to the forward model and

λ_{2}

is the l2-regularzation parameter.

However the reconstructed images based on DA approximation from (1) are not satisfactory mainly for two reasons. First physically DA approximation may not be accurate, especially in the non-diffusive regime, e.g., close to the source, near the boundary, forward-peaking scattering, discontinuous absorption or scattering, or when absorption is smaller or comparable with the scattering. Second, mathematically BLT with DA as forward model is severely ill-posed. Since the boundary measurement provides only angular-averaged data whose dimension does not match with that of the interior region, i.e., number of measurements Mis much smaller than the total number of unknowns $N_{s}$ numerically, it was shown mathematically [3] that there is no unique solution for this inverse source problem based on DA. Moreover, Jcomes from a smoothing diffusion operator which destroys image singularities and features. In practice all the above means that one has to solve a very ill-conditioned and underdetermined system for DA-based BLT. In general, one does not expect a good image quality with clean background from l2-regularized BLT with DA forward model alone without either more data or reduction of unknowns [5–12].

As a remedy for poor reconstruction due to forward modeling by DA, we use radiative transfer equation (RTE) [13–16] as the forward model. Physically, RTE is an accurate model for photon migration in tissues. Moreover, with RTE one can utilize much richer data, such as boundary angular-resolved measurements, whose dimension can match that in the interior region. Mathematically, it was shown [17,18] that the RTE-based inverse source problem with angular-resolved boundary measurements almost always has a unique solution with certain stability. Numerically it means that a better conditioned and less underdetermined system resulted from RTE-based BLT than from DA-based BLT. So in those applications in non-diffusive regime, it is more appropriate to use RTE model and better image reconstruction is expected.

A powerful and general framework in dealing with inverse problem is through variational formulation. It usually contains two key ingredients in the energy functional, data fidelity, which is to match the solution with measurements based on an appropriate forward model, and regularization, which is to regularize an ill-posed problem. With probability interpretation, the form of data fidelity should be determined by noise model and the form of regularization should be determined by a priori knowledge. The balance between this two is usually determined by signal to noise ratio of the data and the smallest scale one wants to resolve in the reconstruction.

In general l2 regularization tends to penalize large elements and create spurious small elements which smears the reconstructed image and produces noisy background when the inverse problem is severely ill-posed. On the other hand, l1 regularization weights all elements equally and tends to produce a sparse solution with clean background which would be a suitable choice for sparse sources [19–22]. As a result, we studied l1-regularization strategy in a multilevel framework and showed its superiority over the conventional l2 regularization for sparse sources in our previous work [23] in which we minimize

\min_{Q} | | J Q - X^{m} | |_{2}^{2} + λ_{1} | | Q | |_{1} .

However, l1 regularization tends to find sparse solutions, which on one hand, makes it successful for sparse source reconstruction, but on the other hand, may shrink the support for non-sparse sources. Therefore, it is desirable to design more appropriate formulations for non-sparse sources that also inherit the merits of l1 regularization. In this paper we propose and compare a few different choices of data fidelity and regularization combinations for RTE-based BLT and justify their performances in different scenarios.

For nearly piecewise-constant source distribution in BLT, a natural choice is to use total variation (TV) regularization, which was first proposed in [24], and has wide applications in various aspects in image processing. Now we consider the following minimization problem

\min_{Q} | | J Q - X^{m} | |_{2}^{2} + λ_{T V} | | Q | |_{T V},

in which

λ_{T V}

is the TV-regularization parameter and

| | \cdot | |_{T V}

denotes the TV norm that we shall specify in the next section.

It is well known that TV regularization [25] tends to be effective for piecewise-constant reconstruction and can preserve the boundary of large object well while removing small features and thin or small objects such as small inclusions or sparse sources. Therefore we propose a combination of l1 and TV regularization which has an interesting geometric balance between them which will be explained in the next section. We shall call it l1 + TV regularization, so that it combines the merits from both l1 and TV regularization in the sense that it successfully reconstructs piecewise-constant image that preserve both geometry and details with little background noise for both sparse and non-sparse sources. The optimization problem becomes

\min_{Q} | | J Q - X^{m} | |_{2}^{2} + λ_{1} | | Q | |_{1} + λ_{T V} | | Q | |_{T V} .

However, if the source distribution is smooth, it is well known that H1 seminorm, which is l2 norm of the gradient, is a good regularization. For example, H1 regularization with l2 data fidelity was used for RTE-based optical tomography in [26]. For comparison purpose, we present mathematical formulation and simulation results with H1 regularization in Section 3.

Most often in practice the data may be contaminated by various types of noises. For Gaussian type of noise, the l2 data fidelity term is a good and natural choice. However, l2 data fidelity term does not handle outliers well since it is biased towards large residuals. On the other hand, nonsmooth data fidelity [27,28] can reconstruct correct images under some assumptions. As a successful example of nonsmooth data fidelity, l1 data fidelity was first introduced in [29], and later found applications in medical imaging [30,31], computer vision [32] and image processing [33]. Therefore, we consider the following optimization problem with l1 data fidelity to deal with practical noise involving outliers

\min_{Q} | | J Q - X^{m} | |_{1} + λ_{1} | | Q | |_{1} + λ_{T V} | | Q | |_{T V} .

Besides, from our simulations, we find that l1 data fidelity with l1 + TV regularization tends to reconstruct as good results as l2 data fidelity with H1 regularization for the smooth non-piecewise-constant sources.

Finally, our multilevel reconstruction developed in [23] can be implemented for all the above formulations. If the source is sparse, coarse mesh reconstruction, which costs much less computation effort and is better conditioned, not only provides a good initial guess but also a good localization for the fine mesh. Hence, the reconstruction on fine mesh can be restricted to the localized region, which again reduces the computation cost. As another consequence, the amount of measurements can be reduced. On the other hand, if the source is not sparse, coarse mesh solution can still provide a good initial guess for the fine mesh and improves overall computational efficiency and stability. Moreover, coarse mesh reconstruction can also help to choose better rto balance between l1 and TV norm for the fine mesh reconstruction.

Here is the outline of this paper. In section 2, we shall first briefly introduce RTE-based BLT, define TV norm on triangular mesh and describe our interior-point methods for solving (4) and (5); in section 3, we shall validate our algorithms with simulations; in section 4, we shall conclude and comment on our future work.

2. Methods

2.1. RTE-based BLT

In this paper, we use RTE as forward model for reasons described in the introduction. The details of an efficient solver of RTE can be found in [34]. Next we shall introduce notations for presentation purpose.

For forward problem, in vivo light propagation is modeled by the following RTE with the vacuum boundary condition

\begin{matrix} \hat{s} \cdot \nabla Ψ (\vec{r}, \hat{s}) + (μ_{a} + μ_{s}) Ψ (\vec{r}, \hat{s}) = μ_{s} \int_{{\hat{S}}^{n - 1}} f (\hat{s}, {\hat{s}}^{'}) Ψ (\vec{r}, {\hat{s}}^{'}) d {\hat{s}}^{'} + q (\vec{r}, \hat{s}) \\ Ψ (\vec{r}, \hat{s}) = 0, if \hat{s} \cdot \hat{n} < 0 \end{matrix}

where the quantities in the equation are the photon fluxΨwhich depends on both space

\vec{r}

and angle

\hat{s}

, the light sourceq, the Henyey-Greenstein scattering function fwith anisotropic factorg, the absorption coefficient

μ_{a}

, the scattering coefficient

μ_{s}

, the outer normal

\hat{n}

at boundary and the angular space

{\hat{S}}^{n - 1}

, i.e., the unit circle in 2D or the unit sphere in 3D.

After Discontinuous Galerkin discretization of (6), we numerically denote piecewise-constant isotropic bioluminescent source by $Q = [q_{i}]$ and piecewise-linear photon flux by $Φ = [ϕ_{i j, m}]$ with $1 \leq i \leq N_{s}, 1 \leq j \leq n_{d} + 1, 1 \leq m \leq N_{a}$ , where $N_{a}$ , $N_{s}$ and $n_{d}$ are the number of angles, the number of spatial mesh elements and problem dimension respectively. On the other hand, we use ${P_{j}, 1 \leq j \leq M}$ to denote measure operators for reconstruction. In this paper, we shall use two types of measurements, i.e., boundary angular-resolved measurementΨand boundary angular-averaged measurement $\int_{\hat{s} \cdot \hat{n} > 0} \hat{s} \cdot \hat{n} Ψ d \hat{s}$ . For RTE-based BLT, which is a linear inverse source problem, the Jacobian matrixJis composed of Green’s function of RTE. We refer readers to [23] for the computation of Jacobian J. Just emphasize here that since we solve RTE without forming the system matrix, the computation for adjoint solver can be realized in three steps: first reverse the directions for detectors and use them as sources due to reciprocity, then compute it in the same way as for direct solver, and last reverse the directions of photon fluxes to get the adjoint fluxes.

2.2. Regularization by TV

Although l1 regularization performs well for sparse sources and is capable of recovering high-resolution details as demonstrated in [23], it tends to shrink the support of large or non-sparse sources, especially in the isotropic scattering regime or diffusive regime. To overcome this, we introduce the regularization by TV norm $| | Q | |_{T V}$ for nearly piecewise-constant source reconstruction.

Since our algorithm is discretized on triangular mesh instead of Cartesian grid, we shall use a discretized version of TV coarea formula [35] for piecewise-constantQ. Let ${Γ_{k}, 1 \leq k \leq N_{e}}$ be edges, then we have

| | Q | |_{T V} = \sum_{k = 1}^{N_{e}} L_{k} | q_{}^{k, l} - q_{}^{k, r} |,

in which

L_{k}

denotes the length of the kth edge, and

q_{}^{k, l}

and

q_{}^{k, r}

represent the left element and the right element with respect to the kth edge with the convention that the spatial index of the left element

q_{}^{k, l}

is always smaller than that of right element

q_{}^{k, r}

to avoid the repetitive counting of edges.

2.3. Optimization Algorithm for l1 and TV regularization

In this section, we shall describe an interior-point method [36] for solving (4), which can be easily modified for minimization problem (2) and (3).

First we recast (4) as a constrained convex optimization with differentiable functional

\min_{Q, u, v} | | J Q - X^{m} | |_{2}^{2} + λ_{1} (\sum_{i = 1}^{N_{s}} A_{i} u_{i} + \sum_{k = 1}^{N_{e}} L_{k} v_{k}),

with constraints

- u_{i} \leq q_{i} \leq u_{i}, 1 \leq i \leq N_{s}

and

- v_{k} \leq q_{}^{k, l} - q_{}^{k, r} \leq v_{k}, 1 \leq k \leq N_{e}

. Here

A_{i}

is area of the ith element and

r : = λ_{T V} / λ_{1}

. Here we incorporate area weights

A_{i}

and edge weights

L_{k}

into the scheme to account for non-uniformity of the mesh, which is not necessary on a uniform mesh.

Next we transform (8) into an sequence of t-dependent unconstrained convex optimization problems

Φ_{t} (Q, u, v) = | | J Q - X^{m} | |_{2}^{2} + [λ_{1} \sum_{i = 1}^{N_{s}} A_{i} u_{i} + \frac{1}{t} ϕ_{1} (Q, u)] + [λ_{T V} \sum_{k = 1}^{N_{e}} L_{k} v_{k} + \frac{1}{t} ϕ_{T V} (Q, v)]

by penalizing the element and edge constraints corresponding with logarithmic barrier functions

ϕ_{1} (Q, u) = - \sum_{i = 1}^{N_{s}} (\log u_{i}^{+} + \log u_{i}^{-})

and

φ_{T V} (Q, v) = - \sum_{k = 1}^{N_{e}} (\log v_{k}^{+} + \log v_{k}^{-})

with

u_{i}^{+} = u_{i} + q_{i}

,

u_{i}^{-} = u_{i} - q_{i}

,

v_{k}^{+} = v_{k} + q_{}^{k, l} - q_{}^{k, r}

and

v_{k}^{-} = v_{k} - q_{}^{k, l} + q_{}^{k, r}

.

The standard interior point method for solving constrained convex problem (8) is through minimizing a sequence of unconstrained convex problems (9), for which the solution converges as $t \to \infty$ . Actually the solution of (9) is no more than $2 (N_{s} + N_{e}) / t$ -suboptimal. This algorithm is what so called barrier method (algorithm 11.1) [36].

For subproblems (9), a simple gradient descent method (algorithm 9.3) [36] is implemented with the efficient computation of the search direction by a preconditioned conjugate gradient method (PCG) [37–39] and backtracking line search algorithm [36,40].

Next we shall describe our PCG method for computing the search direction by solving

\nabla^{2} Φ_{t} [\begin{matrix} Δ Q \\ Δ u \\ Δ v \end{matrix}] = - \nabla Φ_{t} .

Before that, we first compute $\nabla Φ_{t}$ and the Hessian $\nabla^{2} Φ_{t}$ , i.e.,

\nabla Φ_{t} = 2 J^{T} (J Q - X^{m}) + λ_{1} {[A]}_{u} + λ_{T V} {[L]}_{v} + \frac{1}{t} ([\nabla ϕ_{1}] + [\nabla ϕ_{T V}]),

\nabla^{2} Φ_{t} = 2 J^{T} J + \frac{1}{t} (\nabla^{2} ϕ_{1} + \nabla^{2} ϕ_{T V}) .

Note the equality above should be understood in the sense the left hand side is assembled from the matrices on the right hand side with correct ordering of unknowns. Here ${[A]}_{u}$ denotes the element area vector for unknownuand ${[L]}_{v}$ for the edge length vector for unknownv. (11) and (12) can be easily implemented numerically.

For the convenience of presentation, we define $f^{+} (a, b) : = 1 / a^{2} + 1 / b^{2}$ , $f^{-} (a, b) : = 1 / a^{2} - 1 / b^{2}$ , $g^{+} (a, b) : = 1 / a + 1 / b$ and $g^{-} (a, b) : = 1 / a - 1 / b$ . And then we assemble the following local vectors and matrices corresponding to each pair of constraints for either some element or edge into both sides of (10).

First, with the first column or row for $q_{i}$ and the second for $u_{i}$ , the local vector and matrix for the pair of constraints corresponding to the ith element are ${[\nabla ϕ_{1}]}_{i} = - [\begin{matrix} α_{i}^{1} \\ α_{i}^{2} \end{matrix}]$ and ${[\nabla^{2} ϕ_{1}]}_{i} = [\begin{matrix} σ_{i}^{1} & σ_{i}^{2} \\ σ_{i}^{2} & σ_{i}^{1} \end{matrix}]$ , $1 \leq i \leq N_{s}$ , with $α_{i}^{1} = g^{-} (u_{i}^{+}, u_{i}^{-})$ , $α_{i}^{2} = g^{+} (u_{i}^{+}, u_{i}^{-})$ , $σ_{i}^{1} = f^{+} (u_{i}^{+}, u_{i}^{-})$ and $σ_{i}^{2} = f^{-} (u_{i}^{+}, u_{i}^{-})$ .

Second, with the first column or row for $q^{k, l}$ , the second for $q^{k, r}$ and the third for $v_{k}$ , the local vector and matrix for the pair of constraints corresponding to thekth edge are ${[\nabla ϕ_{T V}]}_{k} = - [\begin{matrix} α_{k}^{3} \\ - α_{k}^{3} \\ α_{k}^{4} \end{matrix}]$ and ${[\nabla^{2} ϕ_{T V}]}_{k} = [\begin{matrix} σ_{k}^{3} & - σ_{k}^{3} & σ_{k}^{4} \\ - σ_{k}^{3} & σ_{k}^{3} & - σ_{k}^{4} \\ σ_{k}^{4} & - σ_{k}^{4} & σ_{k}^{3} \end{matrix}]$ , $1 \leq k \leq N_{e}$ , with $α_{k}^{3} = g^{-} (v_{k}^{+}, v_{k}^{-})$ , $α_{k}^{4} = g^{+} (v_{k}^{+}, v_{k}^{-})$ , $σ_{k}^{3} = f^{+} (v_{k}^{+}, v_{k}^{-})$ and $σ_{k}^{4} = f^{-} (v_{k}^{+}, v_{k}^{-})$ .

It is easy to see that the Hessian (12) is positive symmetric definite. In order to make PCG method work, we need to find a precondition matrix Pwhich is not only fast to invert, but also a good approximation of the Hessian. Here we choose the sparse matrix $2diag (J^{T} J) + \frac{1}{t} (\nabla^{2} ϕ_{1} + \nabla^{2} ϕ_{T V})$ asP, which is also symmetric positive definite, and thus $M^{- 1} g$ can be efficiently solved by Gauss-Seidel method for any vectorg.

We summarize our minimization algorithm below. With εbetween 0.00001 and 0.001 as stopping criterion, $λ_{1}$ between 0.01 and 1, andrbetween 0.1 and 10 as the ratio of TV norm over l1 norm, our barrier method gives good results from our numerical tests. Also one needs to increase values of εand $λ_{1}$ when noise level is high. Note in backtracking line search, we check the constraint condition for (8) for feasibility.

Algorithm: barrier method for optimization problem (8)
given tolerance $ε > 0$ , $λ_{1}$ and r.
Initialize $μ = 2$ , $t = 1 / λ_{1}$ , $λ_{T V} = r λ_{1}$ , $q = 0$ , $u = 1$ , $v = 1$ .
Repeat optimization of subproblems (9)
- 1. compute the search direction $(Δ Q, Δ u, Δ v)$ by (12) through PCG with preconditionerP, the relative tolerance $ε_{pcg} = \min (0.1, 0.01 (N_{s} + N_{e}) / t)$ and the relative tolerance $ε_{gs} = 0.0001$ for Gauss-Seidel iterations in PCG.
- 2. compute step size sby backtracking line search with $α = 0.01$ , $β = 0.5$ , $s_{0} = 1$ .
- 3. update: $(Q, u, v) = (Q, u, v) + s (Δ Q, Δ u, Δ v)$
- 4. quit if $2 (N_{s} + N_{e}) / t < ε$
- 5. $t = μ t$

Next, we consider the balance between l1 and TV regularization. For source distribution that is piecewise constant with uniform intensityI,TV norm is approximately $I C$ when Cis the perimeter of boundary of the support and l1 norm is approximately $A I$ , whereAis the area of the support. So in this case a good choice of ratior, should be roughly proportional to $A / C$ . Note $A / C$ can be independent of the length unit since one can always nondimensionalize RTE (6) so that $A / C$ is dimensionless, to whichrshould be equal. This implies that ratiorgives a geometric balance of l1 norm and TV norm. For example, l1 regularization should weight more than TV regularization for sources with smaller ratio $A / C$ , e.g., sources with small or thin support, which corresponds to sparse sources. As a result, with sparse sources as a priori, one should put more weights on l1 regularization by choosing smallrin the algorithm to get better results. We use numerical experiments to demonstrate our rule of optimal choices for different configurations in Section 3. On the other hand, our tests seem to suggest that the results are not very sensitive to the choice ofr. From our numerical experiments, the value of ris typically between 0.1 and 10.

Last but not the least, we normalize Jby columns to get $J^{*}$ , and then actually solve the following instead of (8)

\min_{Q^{*}} | | J^{*} Q^{*} - X^{m} | |_{2}^{2} + λ_{1} \sum_{i = 1}^{N_{s}} A_{i} u_{i} + λ_{T V} \sum_{k = 1}^{N_{e}} L_{k} v_{k}

This is crucial for the successful reconstruction. Physically this means the total energy of measurements with respect to each unknown is on the same level, while mathematically $diag (J^{*}^{T} J^{*})$ is equal to identity matrix, which makes problem isotropic to all unknowns. After solving (13), we need to change correspondingly from $Q^{*}$ toQ, which reflects the sensitivity of each unknown on measurements. This is also equivalent to optimization in weighted norm.

2.4. Optimization Algorithm for l1 Data Fidelity

To simplify the discussion, we will only present the treatment of l1 data fidelity term in the optimization problem (5), and the l1 or TV regularization can be easily adapted to the problem (5) in the same way as what we have done for (4) by (8) in the previous subsection.

Following the standard interior-point method, we first transform the problem

\min_{Q} | | J Q - X^{m} | |_{1}

into the following constrained convex optimization problem

\min_{Q} \sum_{j = 1}^{M} y_{j}

with constraints

- y_{j} \leq (\sum_{i = 1}^{N_{s}} J_{i j} q_{i} - X_{j}^{m}) \leq y_{j}, 1 \leq j \leq M

.

Similarly, (15) can be recasted as a sequence of unconstrained convex optimization problems as

Φ_{t} (Q, y) = \sum_{j = 1}^{M} y_{j} + \frac{1}{t} ϕ (Q, y) .

by penalizing the constraints with the logarithmic barrier function

ϕ (Q, y) = - \sum_{j = 1}^{M} (\log y_{j}^{+} + \log y_{j}^{-})

with

y_{j}^{+} = y_{j} + \sum_{i = 1}^{N_{s}} J_{i j} q_{i} - X_{j}^{m}

and

y_{j}^{-} = y_{j} - \sum_{i = 1}^{N_{s}} J_{i j} q_{i} + X_{j}^{m}

..

To find the search direction, we need to compute

\nabla^{2} Φ_{t} [\begin{matrix} Δ Q \\ Δ y \end{matrix}] = - \nabla Φ_{t},

for which we have

\nabla Φ_{t} = - [\begin{matrix} J^{T} diag (α_{}^{5}) / t \\ {[α_{}^{6}]}_{y} / t - {[1]}_{y} \end{matrix}]

and the Hessian

\nabla^{2} Φ_{t} = \frac{1}{t} [\begin{matrix} J^{T} diag (σ^{5}) J & J^{T} diag (σ^{6}) \\ diag (σ^{6}) J & diag (σ^{5}) \end{matrix}]

with

α_{j}^{5} = g^{-} (y_{j}^{+}, y_{j}^{-})

,

α_{j}^{6} = g^{+} (y_{j}^{+}, y_{j}^{-})

,

σ_{j}^{5} = f^{+} (y_{j}^{+}, y_{j}^{-})

and

σ_{j}^{6} = f^{-} (y_{j}^{+}, y_{j}^{-})

,

1 \leq j \leq M

. Here

{[1]}_{y}

is the unit vector for y,

{[α_{}^{6}]}_{y}

is the constraint vector foryand

diag (\cdot)

is the diagonal matrix with diagonal elements in the bracket.

To invert (17) efficiently, we use a similar PCG method with a sparse symmetric positive definite preconditioner $P = [\begin{matrix} diag (J^{T} diag (σ^{5}) J) & 0 \\ 0 & diag (σ^{5}) \end{matrix}]$ .

Last, we want to emphasize in our algorithm we actually minimize (5) with normalized Jacobian instead of (14), i.e.: we solve a sequence of subproblems coming from the combination of (9) and (16) with l2 data fidelity replaced by l1 data fidelity; in the barrier method, we compute $\nabla^{2} Φ_{t} Δ x = - \nabla Φ_{t}$ with. $Δ x : = {[\begin{matrix} Δ Q & Δ u & Δ v & Δ y \end{matrix}]}^{T}$ for search direction.

3. Simulations

For comparison purpose, besides (1) for l2 regularization, we also consider H1 regularization by a discrete version of, but not exactly, H1 seminorm $| Q |_{H_{1}}^{2} : = \sum_{k = 1}^{N_{e}} {(q_{}^{k, l} - q_{}^{k, r})}^{2}$ , with the same definition as in (7), i.e.,

\min_{Q} | | J Q - X^{m} | |_{2}^{2} + \frac{λ_{2}}{2} | Q |_{H_{1}}^{2},

or by both

| Q |_{H_{1}}^{2}

and the usual discrete l2 norm

| | Q | |_{2}^{2} : = \sum_{i = 1}^{N_{s}} q_{i}^{2}

, i.e.,

\min_{Q} | | J Q - X^{m} | |_{2}^{2} + λ_{2} (| | Q | |_{2}^{2} + r_{2} | Q |_{H_{1}}^{2}),

where

r_{2}

is the parameter to balance between l2 norm and H1 norm.

(1) and (18) can be solved through Levenberg-Marquardt method [41–43], in which we iteratively solve

Q^{k + 1} = Q^{k} + {(J^{T} J + λ_{2, k}^{} D)}^{- 1} (J^{T} (J Q^{k} - X^{m})) .

Din (20) comes from l2 or H1 regularization. For l2-regularized problem (1), simply $D = I$ ; for H1-regularized problem (18), Dcomes from H1 norm and can be assembled from the local matrix ${[D]}_{k} = [\begin{matrix} 1 & - 1 \\ - 1 & 1 \end{matrix}]$ for each edge. The detailed implementation can be found in Algorithm 3.16 in [43] with one modification of the gain predicted by linear model as $\frac{1}{2} {(Q_{k + 1} - Q_{k})}^{T} [λ_{2, k}^{} D (Q_{k + 1} - Q_{k}) - J^{T} X^{m}]$ to guarantee a successful automatic adjust of $λ_{2}$ adaptively. Also (20) can be inverted with a similar PCG method in the previous section.

In this section we use various numerical simulations to validate our novel l1 + TV regularization strategy for BLT and compare l1 data fidelity with l2 data fidelity. In all figures shown in this section, reconstructed images are plotted as what they originally are from reconstructions without further image processing, such as thresholding. As pointed out at the beginning, radiative transport equation (RTE) is used as the underlying model for both forward and inverse problem. Data for forward problem is generated on a fine mesh that is different from the fine mesh used for inverse problem. The fast forward solver developed in [34] is used to compute both the forward problem and the Green’s functions for inverse problem. The methods for solving l2-regularized BLT can be found in [23]. In the following figures for l2 regularization except Fig. 14 and 15 , we only plot results from (1) since l2 regularization (18) and (19) give the similar results. In particular we design different simulations to demonstrate the following points.

Fig. 14 Reconstructions with angular-resolved data for non-piecewise-constant source with g = 0.9. Figure (a) is the true source. Figure (b), (c), (d), (e) and (f) are for l2 data fidelity with l2, H1, l1, TV, l1 + TV regularization respectively. Figure (g), (h) and (i) are for l1 data fidelity with l1, TV, l1 + TV regularization respectively.

Abstract

1. Introduction

2. Methods

2.1. RTE-based BLT

2.2. Regularization by TV

2.3. Optimization Algorithm for l1 and TV regularization

2.4. Optimization Algorithm for l1 Data Fidelity

3. Simulations

4. Conclusions and discussions

Acknowledgement

References and links

Cited By

Figures (16)

Equations (20)

Optics Express