Multilevel bioluminescence tomography based on radiative transfer equation Part 1: l1 regularization

Hao Gao; Hongkai Zhao

doi:10.1364/OE.18.001854

1. Introduction

Bioluminescence tomography (BLT) [1–3] is an emerging molecular imaging modality, which can be used to monitor physiological and pathological activities at molecular levels, specially visualize primary tumor growth and tumor cell metastasis. In BLT [4–13], one tries to recover the location and intensity of isotropic bioluminescent sources from boundary measurements. However, it remains challenging in this area that it is highly desirable but extremely difficult to improve the image resolution as well as reconstruction stability and efficiency in biomedical and clinical practice.

Regarding the forward model of in vivo photon migration, radiative transfer equation (RTE) [14–17] is a well-accepted golden standard and needs to be appropriately discretized numerically to handle complex geometry [18]. In practice, due to its computational burden as a multi-dimensional problem involving both angular directions and spatial variables, the most popular forward model is diffusion approximation (DA) [19], a lower-order approximation of RTE. However, the validation of DA breaks down in small animal imaging, for example, in the region close to light sources or with relatively low scattering. Recently, we have developed a fast forward solver to reduce the computational complexity so that numerical RTE solver is accessible in practice. We will present numerical examples to show that BLT based on RTE can achieve satisfactory resolution, which is usually not the case for BLT based on DA.

In this work, we assume that the optical property of the underlying medium is known. So BLT becomes a linear inverse source problem. In practice the main challenge is due to severe ill-posedness and large degrees of freedom for inverse problem. The ill-posedness may come from two sources: (1) insufficient measurements, such as angular-averaged data, that cause an underdetermined system; (2) ill-conditioning of the system matrix whose columns correspond to Green’s functions of RTE. The ill-posedness makes inversion or reconstruction sensitive to small perturbation in data, such as measurement error and noise. Moreover, the computation cost for computing the Green’s function for RTE increases dramatically and the inversion becomes even more ill-posed when degrees of freedom are increased. This type of problem is typical in many inverse problems. Construction of efficient, accurate and stable reconstruction algorithms is problem-dependent and challenging.

In our application we will explore sparsity structure in the inverse problem. For many BLT applications, the bioluminescent source is sparsely distributed in space. We propose a multilevel imaging method to fully take advantage of the sparsity to improve both efficiency and stability. In this multilevel framework, we first reconstruct the sparse source efficiently on a coarse mesh with good localization, which provides a good initial guess as well as a small region of interest (ROI) for the fine reconstruction. And next reconstruction on fine mesh is carried out within the localized ROI with dramatically reduced degrees of freedom. Through such coarse-to-fine reconstructions, one can stably extract high-resolution images while substantially reducing the computational time. As another consequence, the amount of measurements can also be reduced by utilizing source sparsity in our multilevel reconstruction.

The first crucial step in our multilevel approach is efficient and stable reconstruction on the coarse mesh with good localization property, which is equivalent to finding the sparse solution to the linear system defined on the coarse mesh. How to find the sparse solution efficiently, accurately and stably with as few measurements as possible is at the heart of compressive sensing methods. Compressive sensing has been studied for decades and recently has received a lot of attention due to theoretical justification and many practical implementations, see [20–25] and references therein. The key idea is that, given the sparsity of the signal, an efficient and stable algorithm through l1 regularization can be applied to recover or approximate the original signal with appropriate measurements that have cardinality comparable to the sparsity of the signal. In statistics community, a large body of recent work on least absolute shrinkage and selection operator (LASSO) is also on l1 regularization for sparse model selection, see [26,27] and references therein. However, in order to recover the sparse signal accurately and stably, current theory requires some nice property of the underdetermined system or the measurement matrix. Essentially it requires near-orthogonality or incoherence among columns of the system matrix. For example, the commonly used restricted isometry property (RIP) requires well-conditioning for all submatrices of size comparable to the sparsity or near-orthogonality among columns. RIP is used for both recoverability and stability results in [22–25]. Although a non-RIP result is established in [28] by using characterization of null space of the measurement matrix, the stability results still depend on the conditioning of the system. In recent work on LASSO from [26], an irrepresentable condition is proposed for consistent model selection by LASSO with the right amount of l1 regularization. However, the irrepresentable condition is again related to angle condition between columns of the measurement matrix that are used in the true model and the rest of the columns that are irrelevant. These requirements cannot be satisfied by our problem as well as most inverse problems due to the ill-posedness. In this work we show that by combining l1 regularization with the multilevel approach we can still take advantage of the sparsity and improve computation efficiency, stability and reconstruction results while requiring fewer measurements.

The inverse source problem for RTE has been studied in several papers in mathematics and engineering communities [29–38]. Here we approach this inverse source problem directly by taking advantage of the source sparsity with both l1 regularization and multilevel reconstruction. As we finish up this paper, we recognize a similar l1-regularized strategy for BLT was studied recently in [39], in which BLT based on DA with regularization by differentiable approximation of l1 norm instead of l1 norm was studied due to the difficulty of dealing with non-differentiable l1 norm. Also adaptive mesh rather than multi-level mesh has been studied in the context of adaptive finite element method for BLT [57]. However, our l1-regularized multilevel approach is entirely different although we both use “multilevel”.

Here is the outline of this paper. In section 2, we shall first introduce notations and discretization methods for RTE forward solver, formulate BLT as a linear inverse source problem, describe l1 regularization strategy, and then summarize our algorithms in the multilevel framework; 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. Forward problem

In this paper, we model light propagation by 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} \in Ω, \hat{s} \in {\hat{S}}^{n - 1} \\ Ψ (\vec{r}, \hat{s}) = 0, \vec{r} \in \partial Ω, \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 bioluminescent sourceq, the Henyey-Greenstein scattering function fwith anisotropic factorg, the absorption coefficient

μ_{a}

, the scattering coefficient

μ_{s}

, the outer normal

\hat{n}

at boundary

\partial Ω

of domainΩ and the angular space

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

, i.e., the unit circle in 2D or the unit sphere in 3D. The methods and algorithms in this paper can be extended easily to RTE with reflection boundary condition. We shall only describe RTE with vacuum boundary condition for the presentation purpose.

An efficient solver of (1) was proposed in [18] with angular discretization by finite element method and spatial discretization by piecewise-linear Discontinuous Galerkin method. We design multigrid methods in both angular and spatial space with improved source iteration as relaxations to achieve fast convergence in various regimes. We shall briefly introduce notations and formulate forward problem as follow.

Suppose RTE is discretized with $N_{a}$ angular directions, and let $T : = {τ_{i}, i = 1, \dots, N_{s}}$ be a triangulation of the domain with the dimension $n_{d}$ , ${φ_{i j} (k) = 1_{τ_{i}} δ_{j k}, i = 1, \dots, N_{s}, j, k = 1, \dots, n_{d} + 1}$ be spatial nodal linear bases in DG, and ${{φ^{'}}_{i} = 1_{τ_{i}}, i = 1, \dots, N_{s}}$ be spatial piecewise-constant bases for q. With our discretization, the photon flux $Φ_{m} = \sum_{i = 1}^{N_{s}} \sum_{j = 1}^{n_{d} + 1} ϕ_{i j, m} φ_{i j}, m = 1, \dots, N_{a}$ and can be represented as the column vector $Φ = [ϕ_{i j, m}]$ . Note in BLT, we assume the unknown bioluminescent sourceqis isotropic, i.e., $q = \sum_{i = 1}^{N_{s}} q_{i} {φ^{'}}_{i}$ , and can be represented by the column vector $Q = [q_{i}]$ .

Then we discretize (1) into the following matrix form

(A + B - C) Φ = D

where

{[A]}_{i j m, i^{'} j^{'} m^{'}} = δ_{m m^{'}} δ_{i i^{'}} (- \int_{τ_{i}} φ_{i j} (\hat{s} \cdot \nabla φ_{i^{'} j^{'}}) d \vec{r} + \int_{\partial τ_{i}} \hat{Φ} φ_{i^{'} j^{'}} \hat{s} \cdot \hat{n} d Γ)

with the upwind flux

\hat{Φ}

,

{[B]}_{i j m, i^{'} j^{'} m^{'}} = δ_{m m^{'}} δ_{i i^{'}} \int_{τ_{i}} (μ_{a} + μ_{s}) φ_{i j} φ_{i^{'} j^{'}} d \vec{r}

,

{[C]}_{i j m, i^{'} j^{'} m^{'}} = δ_{i i^{'}} w_{m m^{'}} \int_{τ_{i}} μ_{s} φ_{i j} φ_{i^{'} j^{'}} d \vec{r}

with the angular weight

w_{m m^{'}}

and

{[D]}_{i j m} = \int_{τ_{i}} q_{i} φ_{i j} d \vec{r}

. For details, we refer readers to [18]. Note that in practice we do not explicitly form matrices in (2) since our method is local in the sense that (2) can be solved by iteratively inverting

n_{d} + 1

by

n_{d} + 1

matrices with correct ordering of elements. Here we write RTE in the matrix form (2) for the convenience to introduce our reconstruction algorithms.

Last, we use column vectors ${P_{j}, j = 1, \dots, M}$ of the same structure asΦto denote measure operators for reconstruction. In this paper, we shall use three kinds of measurements, i.e., boundary angular-resolved measurementΨ, boundary angular-averaged measurement $\int_{\hat{s} \cdot \hat{n} > 0} \hat{s} \cdot \hat{n} Ψ d \hat{s}$ and internal angular-averaged measurement $\int Ψ d \hat{s}$ . Each type of measurement can be numerically represented by the vector product $P^{T} Φ$ , whereTdenotes transpose.

2.2. RTE-based BLT

In this section, we describe RTE-based BLT through Green’s functions.

Since the light source is assumed to be isotropic, let $G (\vec{r}, \hat{s}; {\vec{r}}^{'})$ be the Green’s function of RTE (1) for point source in space, then the solution of RTE is $Ψ (\vec{r}, \hat{s}) = \int G (\vec{r}, \hat{s}; {\vec{r}}^{'}) q ({\vec{r}}^{'}) d {\vec{r}}^{'}$ . Thus, for boundary angular-resolved data, at $({\vec{r}}_{j}, {\hat{s}}_{j})$ for each measurement, we have

\sum_{i = 1}^{N_{s}} G ({\vec{r}}_{j}, {\hat{s}}_{j}; {\vec{r}}_{i}) q_{i} = Ψ ({\vec{r}}_{j}, {\hat{s}}_{j})

Mathematically, we need enough angular-resolved data in (3) to match the dimension of unknowns. However, what is usually available in practice are boundary angular-averaged data, for which we have

\sum_{i = 1}^{N_{s}} [\int_{\hat{s} \cdot \hat{n} > 0} \hat{s} \cdot \hat{n} G ({\vec{r}}_{j}, \hat{s}; {\vec{r}}_{i}) d \hat{s}] q_{i} = \int_{\hat{s} \cdot \hat{n} > 0} \hat{s} \cdot \hat{n} Ψ ({\vec{r}}_{j}, \hat{s}) d \hat{s}

Due to dimension mismatch between unknowns and measurements in (4), no unique inversion can be defined mathematically. Numerically, insufficient data lead to a underdetermined linear system, especially on the fine mesh. In this work we apply the proposed multilevel reconstruction algorithm to both cases of (3) and (4) and show that computation efficiency, image resolution and reconstruction stability can be achieved when the source is sparse. In particular by exploring sparsity with multilevel approach one can alleviate the ill-posedness, improve the efficiency and reduce the amount of measurements.

Now we denote measurements byX, Green’s functions byGand write (3) and (4) in the following matrix form for BLT

G Q = X

In particular, the ith column of Gcorresponds to measurements from the unit source sitting onith element and thejth row of Gcorresponds to adjoint solutions of RTE fromjth measurement operator.

The main difficulty to invert the linear system (5) with boundary data is due to its ill-posedness. First (5) is underdetermined with insufficient data, i.e., $M < < N_{s}$ , especially in fine reconstructions with limited angular-averaged data. Second, matrix Gis ill-conditioned, i.e., columns of Gfor adjacent sources are almost parallel since adjacent sources, especially those in depth, produce the similar measurements at the boundary. This can also be seen through resolution analysis [32,40]. Moreover, Gis even more ill-conditioned on finer mesh.

As a result, in BLT, the reconstructed image is usually blurred with low resolution, and hence does not provide good localization of the source; lots of measurements are essential for a stable reconstruction, which requires costly experiments and intensive computation; the reconstruction is usually sensitive to noise.

To deal with the ill-posedness, proper regularization is crucial, but subtle for inverse problem. A good choice of regularization should take into account a priori knowledge of the inverse problem. Based on sparsity assumption in our case, we choose l1 regularization and combine with a multilevel approach to explore source sparsity. For comparison reason we first present the conventional l2 regularization.

Before doing that, since the analytic Green’s function is only available for some simple geometry and it has to be computed numerically in most cases, first we shall start from (2) and derive discrete version of (5) in the next subsection.

2.3. Computation of numerical Jacobian

Let $Χ^{m} = {x_{j}^{m}, j = 1, \dots, M}$ be the measured data for reconstruction and $Χ = {x_{j}, j = 1, \dots, M}$ be the corresponding measurements from our forward solver, i.e., $x_{j} = P_{j}^{T} Φ (Q)$ . Since the solutionΦlinearly depends onQ, we have $x_{j} = \sum_{i = 1}^{N_{s}} \frac{\partial P_{j}^{T} Φ}{\partial q_{i}} q_{i}$ for eachj. Thus $Χ = J Q$ , with matrix $J : = [P_{j}^{T} \frac{\partial Φ}{\partial q_{i}}]$ , which is so called Jacobian or sensitivity matrix. Next we shall compute Jacobian.

Differentiating both sides of (2) with respect to $q_{i}$ , we have $F \frac{\partial Φ}{\partial q_{i}} = \frac{\partial D}{\partial q_{i}}$ with $F : = A + B - C$ . Then,

J_{i j} = P_{j}^{T} (F^{- 1} \frac{\partial D}{\partial q_{i}})

for which, we have

{[\frac{\partial D}{\partial q_{i}}]}_{i^{'} j^{'} m^{'}} = δ_{i i^{'}} \int_{τ_{i^{'}}} q_{i^{'}} φ_{i^{'} j^{'}} d \vec{r}

that is sparse with onlyMnonzero entries.

Now we have two methods to compute $J_{i j}$ . The first method is the direct method by (6), for which we compute $F^{- 1} \frac{\partial D}{\partial q_{i}}$ , i.e., solve the forward problem $N_{s}$ times. The second method is the adjoint method by

J_{i j} = (P_{j}^{T} F^{- 1}) \frac{\partial D}{\partial q_{i}} = Ψ_{j}^{T} \frac{\partial D}{\partial q_{i}}

where

Ψ_{j}^{} = {(F^{T})}^{- 1} P_{j}

, i.e., the adjoint solution of RTE with

P_{j}

as the source. Thus for adjoint method by (7), we need to compute Mforward problems.

In the conventional setting for BLT with boundary angular-averaged data, Mis usually much smaller than $N_{s}$ and the adjoint method is preferred in this case. However, there are two cases that $N_{s}$ is smaller or much smaller thanM. First is when a large amount of data is available, such as boundary angular-resolved data or internal data. Second occurs in our multilevel approach for sparse source reconstruction since multilevel method allows us to reduce the number of unknowns dramatically on the fine mesh through good localization by coarse reconstruction so that the direct method may be favorable, and thus greatly reduces the number of forward problem computation

Regarding the computational time, the most time-consuming part is to compute the JacobianJ. However, it can be done in parallel, which we will implement in the future. On the other hand, it is obvious that the computation cost of (6) is proportional to the number of measurements. So it is highly desirable to have an accurate and stable reconstruction with as few measurements as possible.

2.4. l2 regularization

The popular way for BLT is to formulate it as a least-square minimization problem with an extra l2 regularization term, i.e.,

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

where

λ_{2}

is the parameter to balance the data fitting and l2 regularization.

In general, a good choice of $λ_{2}$ is not known as a priori. Therefore, we choose a nonlinear optimization method through Levenberg-Marquardt method [41,42] to automatically adjust $λ_{2}$ adaptively, in which (8) is reduced to a few iterations

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

Note that l2 regularization is embedded in the algorithm with adaptive $λ_{2, k}^{}$ [43].

Although (8) is relatively easy to solve, it does not take advantage of sparsity of source distribution. As a matter of fact the reconstructed image from BLT is usually blurred with low resolution, and does not localize well for sparse source. With insufficient data and noise, the method becomes more problematic, which we will see in simulations.

2.5. l1 regularization

When the source in BLT is sparse, l1 regularization is a natural and viable choice for finding a good approximate sparse solution. More importantly, finding a good sparse solution on a coarse mesh means good localization which is crucial in identifying the region of interest on finer mesh in our multilevel approach. Although the beautiful compressive sensing theory, which says that with both l1 regularization and appropriate measurements whose cardinality is comparable to the sparsity of the unknown one can stably recover the sparse signal, does not apply to our inverse problem due to severe ill-conditioning of the system matrix, we demonstrate that by combining l1 regularization with the multilevel approach we can still take advantage of the sparsity and improve computation efficiency, stability and reconstruction results while requiring fewer measurements.

Therefore, we regularize l1 norm instead of l2 norm as follow

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

where

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

is a spatially weighted l1 norm with area/volume

A_{i}

and takes into account of mesh nonuniformity.

The simulations in next sections show the improvement of l1 regularization (10) over l2 regularization (8) especially in terms of localization and stability when the source is sparse and the amount of measurements is small.

The numerical algorithm for problem (10) is not as easy as for problem (8) since l1 norm is non-differentiable. Since the main purpose of this paper is not about efficient solver of (10), we use a standard interior-point method [44–46] for minimizing (10) and point out that there are many recent works [47–52] in developing efficient algorithms for (10). All these methods can be directly used in our algorithms to improve the optimization efficiency further. Based on the interior-point method we choose for (10), around eighty to ninety percents of the computational time are spent on computing JacobianJ.

In this paper, we solve the convex optimization problem (10) by the barrier method [44]. We shall briefly formulate the method as follow.

First, we replace the l1 norm by inequality constraints as

\min_{Q, u} | | J Q - X^{m} | |_{2}^{2} + λ \sum_{i = 1}^{N_{s}} A_{i} u_{i},

with constraints

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

, and then add the constraints into the minimization problem by logarithmic barrier penalties as

\min_{Q, u} | | J Q - X^{m} | |_{2}^{2} + λ \sum_{i = 1}^{N_{s}} A_{i} u_{i} - \frac{1}{t} \sum_{i = 1}^{N_{s}} [\log (u_{i} + q_{i}) + \log (u_{i} - q_{i})] .

Next we solve a sequence of t-subproblems (12) with increasing t, where the solution $Q_{t}$ of (12) is no more than $2 N_{s} / t$ -suboptimal which implies $Q_{t}$ converges to the optimal point of (11) as $t \to \infty$ .

In general, the computational time for (11) is longer than (9). To reduce the optimization time by (11), a preconditioned conjugate gradient method [45,46,53] is implemented to compute the search direction for backtracking line search.

2.6. Multilevel reconstruction

It is typical for an inverse problem to become significantly more difficult to solve as mesh is refined. This is because not only degrees of freedom are increased as a more underdetermined linear system but also more importantly the system becomes much more ill-conditioned. Even one can afford significant amount of measurements and computation cost, it is often that using a fine mesh does not necessary lead to a high-resolution reconstruction in practice. The key points for our multilevel approach for sparse source reconstruction for BLT are as follow.

First, the inverse problem on the coarse mesh can be solved more efficiently and stably.

Second, the solution from the coarse mesh provides not only a good initial guess but more importantly a good localization for the sparse source with l1 regularization. The localization will allow us to solve an inverse problem restricted to a small region on fine mesh with the significantly reduced degrees of freedom and improved stability.

Last, the amount of measurements for this multilevel approach only needs to provide a good localization on coarse mesh and matches the degree of freedom in the confined region on fine mesh, which is determined by both sparsity of the source and resolution of the fine mesh.

In our multilevel method, we assume the spatial meshes can be refined from one to another. First we describe a two-level algorithm including three steps:

(1) minimize (10) on the coarse mesh $T_{c} : = {τ_{i}^{c}, 1 \leq i \leq N_{c}}$ to obtain coarse source distribution $Q_{c} : = {q_{i}^{c}, 1 \leq i \leq N_{c}}$ ;
(2) identify the support $S_{c}$ of $Q_{c}$ on the coarse mesh by simple threshold of $Q_{c}$ , e.g., with respect to the maximal values, i.e., $τ_{i}^{c} \in S_{c}$ if $q_{i}^{c} > ε | | Q_{c} | |_{\infty}$ , and interpolate the support $S_{c}$ from $T_{c}$ onto the fine mesh $T_{f} : = {τ_{i}^{f}, 1 \leq i \leq N_{f}}$ to find the support $S_{f} : = {τ_{i}^{f}, 1 \leq i \leq {N^{'}}_{f}}$ on the fine mesh;
(3) minimize (10) within the support $S_{f}$ on $T_{f}$ to obtain $Q_{f} : = {q_{i}^{f}, 1 \leq i \leq {N^{'}}_{f}}$ .

For sparse source reconstruction, (10) is preferred and usually we have ${N^{'}}_{f} < < N_{f}$ . In support detection, εis between 0.01 and 0.05, and we actually extend the support a little bit by one or two grids according to the distance function between centers of grids.

Algorithm: multilevel approach for BLT
FOR $i = N_{i}$ to $N_{f}$
Step 1: if $i = N_{i}$ , $S_{i} : = Ω$ is set to be the whole domainΩ, $Q_{i} = 0$ ;
otherwise, find $S_{i}$ by threshold of $Q_{i - 1}$ , and interpolate from $Q_{i - 1}$ to $Q_{i}$ .
Step 2: compute Jacobian with either (6) or (7).
Step 3: reconstruct $Q_{i}$ on $T_{i}$ by solving (8) or (10).
END

Next, we describe our multilevel algorithm in detail, which also incorporates the flexibility from the dual-mesh scheme [54], i.e., the mesh for unknowns is not necessarily the same as the mesh for solving RTE in the reconstruction. The advantage of the dual-mesh scheme is that we not only improve the reconstruction stability by doing coarse reconstruction and also keep the accuracy of forward solver by solving it on fine mesh.

We start to formulate our multilevel method by defining a sequence of meshes ${(n_{a}^{i}, n_{s}^{i}, n_{q}^{i}), N_{i} \leq i \leq N_{f}}$ , in which $n_{a}^{i}$ , $n_{s}^{i}$ , $n_{q}^{i}$ denotes angular mesh, spatial mesh for forward solver, spatial mesh for unknowns respectively with the convention that finer mesh has bigger value and $n_{q}^{i} \leq n_{s}^{i}$ . Also we denote $Q_{i}$ as the solution and $S_{i}$ as the support of the solution. Then we have the above multilevel algorithm.

Here we need to modify $\frac{\partial D}{\partial q_{i}}$ . Let $T_{q} : = {τ_{i}, i = 1, \dots, N_{q}}$ be the mesh for unknowns, $T_{f} : = {τ_{i}, i = 1, \dots, N_{f}}$ be the mesh for forward solver with $N_{q} \leq N_{f}$ , and ${I_{}^{i}, i = 1, \dots, N_{q}}$ be spatial mapping from $T_{q}$ to $T_{f}$ , for which $I_{}^{i}$ contains all elements on $T_{f}$ that belongs to some element iof $T_{q}$ . Then we have ${[\frac{\partial D}{\partial q_{i}}]}_{i^{'} j^{'} m^{'}} = 1_{i^{'} \in I_{}^{i}} \int_{τ_{i^{'}}} q_{i^{'}} φ_{i^{'} j^{'}} d \vec{r}$ .

3. Simulations

In this section we use various numerical simulations to validate our l1-regularized multilevel imaging algorithm for BLT. Most simulations are performed in two dimensions (2D), which serves the purposes of this paper. Although we give an example (Fig. 11 ) in three dimensions (3D), we will study 3D in vivo BLT further in more practical setting in our future study. 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 [18] is used to compute both the forward problem and the Green’s functions for inverse problem. In particular we design different simulations to demonstrate the following points.

Fig. 11 3D reconstruction. Figure (a) is the geometry. Figure (b) is the true source. Figure (c) is the reconstructed source.

Download Full Size | PDF

First, for sparse source, l1 regularization localizes the source better than l2 regularization, especially when there is noise and small amount of data. See Fig. 1 , 2 , 4 and 5 .

Fig. 1 l2 reconstructions for a small object with twelve measurements. Figure (a) is the mesh for forward problem. Figure (b) is the mesh for inverse problem. Figure (c) is the true source distribution. Figure (d) is the image from l2 reconstruction.

Download Full Size | PDF

Fig. 2 l1 reconstructions for a small object with twelve measurements. Figures (a), (b), (c), (d), (e) and (f) are reconstructed images by l1 regularization with 0%, 10%, 20%, 30%, 40% and 50% Gaussian noise respectively.

Download Full Size | PDF

Fig. 4 Reconstructions with boundary data when g = 0.9. Figures (a), (c) and (e) are reconstructed images by l1 regularization with 120 angular-averaged measurements, 121 angular-resolved measurements and 1320 angular-resolved measurements respectively. Figures (b), (d) and (f) are reconstructed images by l2 regularization with 120 angular-averaged measurements, 121 angular-resolved measurements and 1320 angular-resolved measurements respectively.

Download Full Size | PDF

Fig. 5 Reconstructions with boundary data when g = 0.0. Figures (a), (c) and (e) are reconstructed images by l1 regularization with 120 angular-averaged measurements, 121 angular-resolved measurements and 1320 angular-resolved measurements respectively. Figures (b), (d) and (f) are reconstructed images by l2 regularization with 120 angular-averaged measurements, 121 angular-resolved measurements and 1320 angular-resolved measurements respectively.

Download Full Size | PDF

Second, our multilevel algorithm is more effective and efficient than solving the inverse problem directly on the fine mesh. See Fig. 6 .

Fig. 6 Comparison of multilevel and single-mesh reconstructions with l1 regularization and boundary data when g = 0.9. Figures (a) and (b) are reconstructed images through only the fine mesh with 120 angular-averaged measurements and 720 angular-resolved measurements respectively. Figures (c) and (e) are reconstructed images on the coarse and fine mesh respectively through the multilevel approach with 60 angular-averaged measurements. Figures (d) and (f) are reconstructed images on the coarse and fine mesh respectively through the multilevel approach with 360 angular-resolved measurements.

Download Full Size | PDF

Third, our l1-regularized multilevel algorithm works well with both angular-resolved data and angular-averaged data. Besides, reconstruction with angular-resolved data has the potential for high-resolution imaging. See Fig. 4 and 5.

Fourth, RTE based model can reconstruct better than DA based model with the same amount of angular-averaged data. See Fig. 7 and 8 .

Fig. 7 Comparison of RTE-based and DA-based reconstruction with homogenous background. Figure (a) is the reconstructed image with RTE when g = 0.0. Figure (b) is the reconstructed image with DA when g = 0.0. Figure (c) is the reconstructed image with RTE when g = 0.9. Figure (d) is the reconstructed image with DA when g = 0.9.

Download Full Size | PDF

Fig. 8 Comparison of RTE-based and DA-based reconstruction with heterogeneous background. Figure (a) is the mesh for forward problem. Figure (b) is the true source. Figure (c) is the reconstructed image with DA when g = 0.0. Figure (d) is the reconstructed image with RTE when g = 0.0.

Download Full Size | PDF

Last, when interior measurements are available, better and more stable reconstruction can be achieved. See Fig. 9 and 10 .

Fig. 9 Comparison of l1-regularized reconstructions with boundary angular-averaged data and internal angular-averaged data when g = 0.9. Figure (a) is the true object. Figure (b) is the mesh to generate the forward data. Figure (c) is the result with 120 boundary angular-averaged data. Figure (d) is the result with 120 random internal angular-averaged data.

Download Full Size | PDF

Fig. 10 l1-regularized reconstructions with internal angular-averaged data when g = 0.9. Figures (a) and (b) are reconstructed images through the multilevel approach with 668 measurements on the coarse mesh and 324 measurements on the fine mesh respectively.

Download Full Size | PDF

The following 2D simulations are performed on a square domain with 20mm side length with homogeneous optical properties $μ_{a} = 0.01 m m^{- 1}$ and $μ_{s} = 1 m m^{- 1}$ unless otherwise indicated.

In the first example (Fig. 1), a small inclusion (c) with 1mm diameter is centered at (5mm, 5mm). The simulated angular-averaged measurements with $M = 12$ is generated with 64 angular directions and a spatial mesh (a) of 697 nodes and 1312 elements, while the reconstruction (b) mesh is of 32 directions, 697 nodes and 1312 elements.

The image from l2 reconstruction without noise is displayed in (d) in Fig. 1. While the images from l1 reconstruction without or with noise are displayed in Fig. 2 with the noise defined as $X^{m}^{'} = X^{m} (1 + σ N)$ , where σis the signal to noise ratio, andNis a Gaussian random variable with unity mean and unity variation. In Fig. 2, we vary σfrom 0 percent to 50 percents with 10 percents increment and plot l1-regularized image. Comparing images from l2 and l1 regularization, we conclude that l1 regularization can efficiently localize the sparse source with a few measurements, even in the presence of large noise, while l2 regularization fails to reconstruct the sparse source.

In the second example (Fig. 3 ), three letter inclusions (a) are put in the domain to evaluate our reconstruction algorithms; the forward mesh (b) has 1509 nodes and 2896 elements; the coarse mesh (c) for reconstructions has 365 nodes and 668 elements and the fine one (d) has 1397 nodes and 2672 elements.

Fig. 3 Meshes for forward and inverse problem. Figure (a) is the true source distribution. Figure (b) is the mesh for forward problem. Figures (c) and (d) are the coarse and fine mesh for inverse problem.

Download Full Size | PDF

In Fig. 4 and Fig. 5, the reconstructions with anisotropic factor g = 0.9 and 0.0 respectively are carried out to demonstrate that l1 regularization ((a), (c) and (e) in Fig. 4 and Fig. 0.5) in general reconstructs better images than l2 regularization ((b), (d) and (f) in Fig. 4 and Fig. 0.5), especially when inverse problem is more severely ill-posed due to fewer available data, e.g., angular-resolved data.

Comparing (a) and (c), we find that angular-resolved data are preferred over angular-averaged data for the similar amount of data, especially in the forward-peaking regime with big anisotropic factor g since anisotropic angular information improves the probing ability. When a rich amount of angular-resolved data is available as in (e), our algorithm with l1 regularization is able to reconstruct high-resolution images.

A justification of this is that boundary angular-resolved data are essentially three-dimensional in 2D, which is enough to match the dimensions of unknown sources. As a result, the reconstruction based on RTE has a great potential for high-resolution images even with only boundary data, while one does not expect a good image from DA-based reconstruction due to the dimension mismatch between the boundary data and unknowns since DA does not contain angular information.

In conclusion, one needs RTE model to alleviate the ill-posedness of standard BLT when only boundary data are available; on the other hand, efficient algorithms and proper regularization, e.g., l1 regularization for sparse source, are essential to recover high-resolution images. Besides, one may already notice a noise dot in the middle of reconstructed images in (a), (c) and (e) in Fig. 5. due to the ill-posedness of inverse problem. However, this dot has small contrast on the coarse mesh, thus can be removed through our multilevel approach by choosing proper threshold since its value is much smaller than that of other sources.

Next, we shall validate the superiority of proposed multilevel approach over the single-mesh reconstruction by two examples as shown in Fig. 6. The first example ((a), (c) and (e) in Fig. 6) is on the reconstructions with angular-averaged data wile the second example ((b), (d) and (f) in Fig. 6) is with angular-resolved data. In both examples, we performed the single-mesh reconstruction with images (a) and (b), and then the multilevel reconstruction with images (c) and (d) on coarse mesh and images (e) and (f) on fine mesh. From Fig. 6, we conclude that, with smaller amount of data, our multilevel approach gives better reconstruction than the single-mesh one. On the other hand, it implies that the multilevel approach is able to reduce the computational time for a given accuracy. Moreover, we may further accelerate the reconstruction by solving much less number of forward problems with the direct method (6) for Jacobian on the fine mesh when the sources are sparsely localized.

In Fig. 7 and 8, we compare the RTE-based reconstructions with DA-based reconstructions in diffusive or non-diffusive regime, for which l1-regularized reconstructions based on either RTE or DA are performed with 120 boundary angular-averaged measurements. To avoid the biased model error, for DA-based reconstructions, the simulation measurements are generated with DA as well.

In Fig. 7, we use the same meshes ((b) and (d) in Fig. 3) respectively for data generation and inverse problem for DA. In DA, we have the reduced scattering coefficient $μ_{s}^{'} = μ_{s} (1 - g)$ with g = 0.0 and 0.9 respectively. From reconstruction results, we conclude that BLT with DA does not produce as good images as with RTE. Actually with anisotropic g or within the domain of a few mean free paths, BLT with RTE produces much better results than with DA. Besides, RTE may utilize angular-resolved data for high-resolution images while DA is unable to incorporate the angular information.

In Fig. 8, we compare on a domain with 20mm side length with heterogeneous optical parameter. The optical parameters are as follow: $μ_{a} = 0.01 m m^{- 1}$ and $μ_{s} = 1 m m^{- 1}$ for the background, $μ_{a} = 1 m m^{- 1}$ and $μ_{s} = 1 m m^{- 1}$ for a circular inclusion centered at (4mm, 4mm) with diameter 6mm, $μ_{a} = 0.01 m m^{- 1}$ and $μ_{s} = 0.01 m m^{- 1}$ for a circular inclusion centered at (10mm, 10mm) with diameter 6mm, $μ_{a} = 0.01 m m^{- 1}$ and $μ_{s} = 0.1 m m^{- 1}$ for a circular inclusion centered at (16mm, 16mm) with diameter 6mm. The unit source ((b) in Fig. 8) is located at six smaller circular inclusions with diameter 2mm centered at (4mm, 10mm), (4mm, 14mm), (4mm, 18mm), (10mm, 4mm), (14mm, 4mm), (18mm, 4mm) respectively. g = 0.0 in this case. The forward mesh ((a) in Fig. 8) has 3792 elements and 1961 nodes; the reconstruction mesh ((d) in Fig. 3) is used.

Clearly the reconstruction is in non-diffusive regime containing both strong absorbing region and void-like region. From reconstruction results as in Fig. 8, we conclude that RTE-based BLT performs better DA-based BLT in the presence of non-diffusive regions. On the other hand, reconstruction with l1 norm tends to give sparse solution, which may have smaller support than the original support as shown in (d). In our next paper, we shall introduce the methods to resolve this issue.

As the last simulations in 2D, we compute two examples in Fig. 9 and 10 when internal angular-averaged data are available.

In Fig. 9, an L-shaped source inclusion (a) is embedded at the center of the square domain with the forward mesh (b) of 1469 nodes and 2816 elements. The fine mesh in the second example above ((d) in Fig. 3) is used for reconstruction. We randomly pick up 120 internal angular-averaged measurements and plot the reconstruction result in (d) while (c) is with 120 boundary angular-averaged measurements. Clearly, internal angular-averaged data give better reconstruction image than boundary angular-averaged data in reconstructing sources in depth.

In Fig. 10, we perform the multilevel method when full internal angular-averaged data are available. However, only partial data are necessary on the fine reconstruction. Again it indeed confirms the superiority of l1-regularized multilevel reconstruction: we first efficiently localize the support of unknowns through l1 minimization on coarse mesh, and then resolve the high-resolution image on fine mesh by optimization only in the region of interest while computing much less number of forward problems than the single-mesh reconstructions.

In conclusion, the internal data will better-pose the inverse problem, which motivates us to explore in the future the potential imaging modality to generate internal data, such as one similar to photoacoustic tomography [55].

Last, we give a 3D example with g = 0.0 in Fig. 11 on a cube (a) of side length 20mm with background $μ_{a} = 0.01 m m^{- 1}$ and $μ_{s} = 1 m m^{- 1}$ . A strong absorbing cubic inclusion of side length 5mm is centered at (15.5mm, 15.5mm, 10mm) with $μ_{a} = 1 m m^{- 1}$ and $μ_{s} = 1 m m^{- 1}$ ; a void-like cubic inclusion of side length 5mm is centered at (4.5mm, 4.5mm, 10mm) with $μ_{a} = 0.01 m m^{- 1}$ and $μ_{s} = 0.01 m m^{- 1}$ ; two isotropic sources are embedded into two cubes with side length 2mm centered at (4mm, 18mm, 7mm) and (10mm, 4mm, 13mm) respectively. The forward mesh for data generation has 258 angular directions and a spatial mesh (b) of 3500 nodes and 17683 elements; the reconstruction mesh has 66 angular directions and a spatial mesh (c) of 2837 nodes and 14592 elements; 128 boundary angular-averaged measurements are used for the reconstruction, and thus 128 adjoint RTE solvers are computed in the reconstruction, which takes roughly three hours in total on a desktop with Intel CPU E6850 (3.00GHz). Comparing true source (b) and reconstructed source (c), we conclude that our algorithm is able to reconstruct the source in a 3D non-diffusive regime, e.g., high absorbing and void-like region. As an extension of the proposed methods in this paper, we shall study 3D RTE-based BLT in details in the more practical setting in our future work. Besides, although our algorithm is pretty efficient even in 3D, it can be accelerated further by parallelizing RTE forward solvers in both forward and inverse problem.

4. Conclusions and discussions

In this paper, we first introduce l1-regularized reconstruction, which is proven to perform in general better than l2-regularized reconstruction, especially for sparse sources with less accessible data and large noise. Second, we propose a multilevel reconstruction approach, which reconstructs better and takes less computational time than the single-mesh reconstruction. With the combination of both, our RTE-based l1-regularized multilevel approach has the potential for high-resolution images with even only boundary data.

On the other hand, we compare reconstructions with boundary angular-average data, boundary angular-resolved data and internal angular-averaged data. With either rich angular-resolved data or internal data, we are able to reconstruct almost the exact light sources. As a result, we shall explore on the potential imaging modalities which make those data accessible.

As a similar problem to BLT, linearized optimal tomography based on perturbation of a priori background can be turned into an inverse source problem exactly as BLT, to which the algorithms here can be directly extended.

It is also worthy of mentioning another benefit from multilevel approach: that is coarse reconstruction is able to provide a good initial guess for fine reconstruction, which is crucial for nonlinear inverse problem, such as optical tomography [56].

Last, the l1 regularization can be very efficient for sparse source reconstruction. However, in many interesting cases, the total variation instead of unknowns itself is sparse. Thus we shall consider regularizing total variation instead in our next paper.

Acknowledgement

The work is partially supported by NSF grant DMS0811254 and ONR grant N00014-02-1-0090.

References and links

1. C. H. Contag and B. D. Ross, “It’s not just about anatomy: in vivo bioluminescence imaging as an eyepiece into biology,” J. Magn. Reson. Imaging 16(4), 378–387 (2002). [CrossRef] [PubMed]

2. G. Wang, E. A. Hoffman, G. McLennan, L. V. Wang, M. Suter, and J. Meinel, “Development of the first bioluminescent CT scanner,” Radiology 229(P), 566 (2003).

3. G. Wang, Y. Li, and M. Jiang, “Uniqueness theorems in bioluminescence tomography,” Med. Phys. 31(8), 2289–2299 (2004). [CrossRef] [PubMed]

4. G. Wang, W. Cong, K. Durairaj, X. Qian, H. Shen, P. Sinn, E. Hoffman, G. McLennan, and M. Henry, “In vivo mouse studies with bioluminescence tomography,” Opt. Express 14(17), 7801–7809 (2006). [CrossRef] [PubMed]

5. W. Cong, G. Wang, D. Kumar, Y. Liu, M. Jiang, L. Wang, E. Hoffman, G. McLennan, P. McCray, J. Zabner, and A. Cong, “Practical reconstruction method for bioluminescence tomography,” Opt. Express 13(18), 6756–6771 (2005). [CrossRef] [PubMed]

6. G. Wang, X. Qian, W. Cong, H. Shen, Y. Li, W. Han, K. Durairaj, M. Jiang, T. Zhou, J. Cheng, J. Tian, Y. Lv, H. Li, and J. Luo, “Recent development in bioluminescence tomography,” Current Medical Imaging Reviews 2(4), 453–457 (2006). [CrossRef]

7. G. Wang, H. Shen, K. Durairaj, X. Qian, and W. Cong, “The first bioluminescence tomography system for simultaneous acquisition of multi-view and multi-spectral data,” Int. J. Biomed. Imaging 2006, 1–8 (2006). [CrossRef]

8. X. Gu, Q. Zhang, L. Larcom, and H. Jiang, “Three-dimensional bioluminescence tomography with model-based reconstruction,” Opt. Express 12(17), 3996–4000 (2004). [CrossRef] [PubMed]

9. A. J. Chaudhari, F. Darvas, J. R. Bading, R. A. Moats, P. S. Conti, D. J. Smith, S. R. Cherry, and R. M. Leahy, “Hyperspectral and multispectral bioluminescence optical tomography for small animal imaging,” Phys. Med. Biol. 50(23), 5421–5441 (2005). [CrossRef] [PubMed]

10. H. Dehghani, S. C. Davis, S. Jiang, B. W. Pogue, K. D. Paulsen, and M. S. Patterson, “Spectrally resolved bioluminescence optical tomography,” Opt. Lett. 31(3), 365–367 (2006). [CrossRef] [PubMed]

11. C. Kuo, O. Coquoz, T. Troy, D. Zwarg, and B. Rice, “Bioluminescent tomography for in vivo localization and quantification of luminescent sources from a multiple-view imaging system,” Mol. Imaging 4, 370 (2005).

12. G. Alexandrakis, F. R. Rannou, and A. F. Chatziioannou, “Tomographic bioluminescence imaging by use of a combined optical-PET (OPET) system: a computer simulation feasibility study,” Phys. Med. Biol. 50(17), 4225–4241 (2005). [CrossRef] [PubMed]

13. Y. Lv, J. Tian, W. Cong, and G. Wang, “Experimental study on bioluminescence tomography with multimodality fusion,” Int. J. Biomed. Imaging 2007, 86741 (2007). [CrossRef]

14. K. M. Case, and P. F. P. F. Zweifel, Linear Transport Theory (Addison-Wesley Educational Publishers Inc., 1967).

15. S. Chandrasekhar, Radiative Transfer (Dover Publications, 1960).

16. E. E. Lewis, and W. F. Miller, Computational Methods of Neutron Transport (Wiley, 1984).

17. A. Ishimaru, Wave Propagation and Scattering in Random Media (Academic Press, 1978).

18. H. Gao and H. K. Zhao, “A fast forward solver of radiative transfer equation,” Transp. Theory Stat. Phys. 38(3), 149–192 (2009). [CrossRef]

19. A. R. Arridge, “Optical tomography in medical imaging,” Inverse Probl. 15(2), R41–R93 (1999). [CrossRef]

20. E. J. Candes and M. B. Wakin, “A introduction to compressive sampling,” IEEE Signal Process. Mag. 25(2), 21–30 (2008). [CrossRef]

21. D. Donoho, “Compresse sensing,” IEEE Trans. Inf. Theory 52(4), 1289–1306 (2006). [CrossRef]

22. E. J. Candes, J. Romberg, and T. Tao, “Robust uncertainty principles: exact signal reconstruction from highly incomplete frequency information,” IEEE Trans. Inf. Theory 52(2), 489–509 (2006). [CrossRef]

23. E. J. Candes, J. Romberg, and T. Tao, “Stable signal recovery from incomplete and inaccurate measurements,” Commun. Pure Appl. Math. 59(8), 1207–1223 (2006). [CrossRef]

24. E. J. Candes and T. Tao, “Decoding by linear programming,” IEEE Trans. Inf. Theory 51(12), 4203–4215 (2005). [CrossRef]

25. E. J. Candes and T. Tao, “Near optimal signal recovery from random projections: universal encoding strategies,” IEEE Trans. Inf. Theory 52(12), 5406–5425 (2006). [CrossRef]

26. P. Zhao and B. Yu, “On model selection consistency of Lasso,” J. Mach. Learn. Res. 7, 2541–2563 (2006).

27. R. Tibshirani, “Regression shrinkage and selection via the Lasso,” J. R. Stat. Soc., B 58, 267–288 (1996).

28. Y. Zhang, “Theory of compressive sensing via l1-minimization: a non-RIP analysis and extensions,” Rice University CAAM Technical Report TR08–11, (2008).

29. G. Bal and A. Tamasan, “Inverse source problems in transport equations,” SIAM J. Math. Anal. 39(1), 57–76 (2007). [CrossRef]

30. A. Charette, J. Boulanger, and H. K. Kim, “An overview on recent radiation transport algorithm development for optical tomography imaging,” J. Quant. Spectrosc. Radiat. Transf. 109(17-18), 2743–2766 (2008). [CrossRef]

31. A. D. Klose, V. Ntziachristos, and A. H. Hielscher, “The inverse source problem based on the radiative transfer equation in optical molecular imaging,” J. Comput. Phys. 202(1), 323–345 (2005). [CrossRef]

32. V. Markel and J. Schotland, “Fourier-laplace structure of the linearized inverse scattering problem for the radiative transport equation,” Inv. Prob. Imag. 1, 181–189 (2007). [CrossRef]

33. N. J. McCormick, “Inverse radiative transfer problems: a review,” Nucl. Sci. Eng. 112, 185–198 (1992).

34. A. N. Panchenko, “Inverse source problem of radiative transfer: a special case of the attenuated Radon transform,” Inverse Probl. 9(2), 321–337 (1993). [CrossRef]

35. C. E. Siewert, “An inverse source problem in radiative transfer,” J. Quant. Spectrosc. Radiat. Transf. 50(6), 603–609 (1993). [CrossRef]

36. P. Stefanov and G. Uhlmann, “An inverse source problem in optical molecular imaging,” Analysis and PDE 1, 115–126 (2008). [CrossRef]

37. Z. Tao, N. J. McCormick, and R. Sanchez, “Ocean source and optical property estimation from explicit and implicit algorithms,” Appl. Opt. 33(15), 3265 (1994). [CrossRef] [PubMed]

38. R. Weissleder and U. U. Mahmood, “Molecular imaging,” Radiology 219(2), 316–333 (2001). [PubMed]

39. Y. Lu, X. Zhang, A. Douraghy, D. Stout, J. Tian, T. F. Chan, and A. F. Chatziioannou, “Source reconstruction for spectrally-resolved bioluminescence tomography with sparse a priori information,” Opt. Express 17(10), 8062–8080 (2009). [CrossRef] [PubMed]

40. M. Boffety, M. Allain, A. Sentenac, M. Massonneau, and R. Carminati, “Analysis of the depth resolution limit of luminescence diffuse optical imaging,” Opt. Lett. 33(20), 2290–2292 (2008). [CrossRef] [PubMed]

41. K. Levenberg, “A method for the solution of certain nonlinear problems in least squares,” Q. Appl. Math. 2, 164–168 (1944).

42. D. W. Marquardt, “An algorithm for least-squares estimation of nonlinear parameters,” J. Soc. Ind. Appl. Math. 11(2), 431–441 (1963). [CrossRef]

43. K. Madsen, H. B. Nielsen, and O. Tingleff, Methods for Non-linear Least Squares Problems (Technical University of Denmark, 1999).

44. S. Boyd, and L. Vandenberghe, Convex Optimization (Cambridge university press, 2004).

45. S. J. Kim, K. Koh, M. Lustig, and S. Boyd, “An efficient method for compressed sensing,” IEEE International Conference on Image Processing 3, 117–120 (2007).

46. S. J. Kim, K. Koh, M. Lustig, S. Boyd, and D. Gorinevsky, “An Interior-Point Method for Large-Scale l1-Regularized Least Squares,” IEEE J. Sel. Top. Signal Process. 1(4), 606–617 (2007). [CrossRef]

47. W. Yin, S. Osher, D. Goldfarb, and J. Darbon, “Bregman iterative algorithms for l1-minimization with applications to compressed sensing,” SIAM J. Imaging Sciences 1(1), 143–168 (2008). [CrossRef]

48. Z. Wen, W. Yin, D. Goldfarb, and Y. Zhang, “A fast algorithm for sparse reconstruction based on shrinkage, subspace optimization and continuation,” Rice University CAAM Technical Report TR09–01, (2009).

49. E. Esser, X. Zhang, and T. Chan, “A general framework for a class of first order primal-dual algorithms for TV minimization,” UCLA CAM Report 09–67, (2009).

50. T. Goldstein and S. Osher, “The split bregman method for l1 regularized problems,” SIAM J. Imaging Sci. 2(2), 323–343 (2009). [CrossRef]

51. S. Osher, Y. Mao, B. Dong, and W. Yin, “Fast linearized Bregman iteration for compressed sensing and sparse denoising,” Commun. Math. Sci. (to be published).

52. J. Yang, Y. Zhang, and W. Yin, “An Efficient TVL1 Algorithm for Deblurring Multichannel Images Corrupted by Impulsive Noise,” SIAM J. Sci. Comput. 31(4), 2842–2865 (2009). [CrossRef]

53. J. Demmel, Applied Numerical Linear Algebra (Cambidge Univ. Press, 1997).

54. K. D. Paulsen, P. M. Meaney, M. J. Moskowitz, and J. R. Sullivan, “A dual mesh scheme for finite element based reconstruction algorithms,” IEEE Trans. Med. Imaging 14(3), 504–514 (1995). [CrossRef] [PubMed]

55. M. Xu and L. V. Wang, “Photoacoustic imaging in biomedicine,” Rev. Sci. Instrum. 77, 041101–1-041101–22 (2006). [CrossRef]

56. H. Gao, and H. Zhao, “A multilevel and multigrid optical tomography based on radiative transfer equation,” in Proceedings of SPIE (Munich, Germany, 2009), pp. 73690E–1-73690E–10.

57. Y. Lv, J. Tian, W. X. Cong, G. Wang, J. Luo, W. Yang, H. Li, and H. Li, “A multilevel adaptive finite element algorithm for bioluminescence tomography,” Opt. Express 14(18), 8211–8223 (2006). [CrossRef] [PubMed]

Multilevel bioluminescence tomography based on radiative transfer equation Part 1: l1 regularization

Abstract

1. Introduction

2. Methods

2.1. Forward problem

2.2. RTE-based BLT

2.3. Computation of numerical Jacobian

2.4. l2 regularization

2.5. l1 regularization

2.6. Multilevel reconstruction

3. Simulations

4. Conclusions and discussions

Acknowledgement

References and links

Cited By

Figures (11)

Equations (12)

Optics Express