Sparse reconstruction for quantitative bioluminescence tomography based on the incomplete variables truncated conjugate gradient method

Xiaowei He; Jimin Liang; Xiaorui Wang; Jingjing Yu; Xiaochao Qu; Xiaodong Wang; Yanbin Hou; Duofang Chen; Fang Liu; Jie Tian

doi:10.1364/OE.18.024825

1. Introduction

Bioluminescence tomography (BLT) is a promising optical imaging technique that offers the possibilities to monitor physiological and pathological processes at cellular and molecular levels in vivo [1–4]. Based on a photon propagation model, the anatomical structure information and the associated tissue optical properties, BLT attempts to reconstruct 3D distribution of the probe inside of a small animal from the detection of light emission on the body surface [5]. Due to the insufficient measurements and highly diffusive nature of the photon propagation within biological tissues, BLT reconstruction is still a very challenging ill-posed inverse problem.

The ill-posedness is the crucial issue that reconstruction algorithms have to address. In most existing reconstructions so far, multi-spectral measurement [6–10] and a priori information of the permissible source region [11–15] are two common ways to reduce ill-posedness of BLT by means of increasing the measurement and/or restricting the region of interest respectively. Although these techniques improve reconstruction qualities to a certain degree, they in turn impose some limitations on practical applications. For example, multi-spectral or spectrally-resolved BLT usually need more signal acquisition time and have a high computational cost. Additionally, the size of the permission region has a serious impact on the reconstruction results [11] and a relatively small permission region is adopted to obtain meaningful results, which is not available in most cases. Recently, an iteratively shrunk permissible region strategy has been incorporated in a two-step reconstruction algorithm to improve the reconstruction quality [16].

No matter which approach is adopted, it is essential to combine regularization techniques to overcome the ill-posedness when trying to recover source distribution from noisy measurements. Most regularization methods used in BLT, such as Tikhonov regularization, try to stabilize the problem by achieving a trade-off between a loss term and an ℓ ₂-norm regularization term [11–14]. However, due to the inherent characteristic of the ℓ ₂ norm, these methods usually produce smooth solutions which will result in the loss of high-frequency parts of the original signal and introduce some artifacts in the reconstruction images.

Over the past years, compressed sensing (CS) has attracted more interests especially in signal processing due to its ability of capturing salient information of sparse or compressible signals at a rate significantly below the demands of the Nyquist sampling theorem [17,18]. Coincidently, in practical BLT applications, the bioluminescent source distribution is usually sparse and only insufficient boundary measurements are available. In view of the characteristics in BLT scenarios, the CS techniques potentially bring benefits in spatial resolution and algorithm stability to BLT reconstruction. Recently, the CS inspired methods have emerged in BLT with different photon propagation models. Lu et al. presented a reconstruction algorithm for a spectrally-resolved BLT based on the diffusion approximation (DA) model where the ℓ ₂ norm in Tikhonov regularization was replaced by the ℓ ₁ norm and a limited memory variable metric optimization method was used to solve the bound constrained BLT inverse problem [19]. Cong et al. proposed a reconstruction method of bioluminescence sources based on the phase approximation model, in which the CS technique is applied to regularize the inverse source recovery [20]. Gao et al. developed an ℓ ₁-regularized multilevel FEM approach for BLT based on the radiative transfer equation (RTE) [21]. Based on the DA model, we applied a truncated Newton interior-point method (ℓ1-ℓs) on multilevel adaptive meshes for solving the BLT inverse problem, which can achieve simultaneous recovery of density and power as well as accurate source location with the gradually reduced permissible region [22,23]. These works demonstrated the feasibility and potential of the CS techniques with numerical simulation and phantom experiments. However, all of these previous reconstructions required either permissible region constraint or multi-spectral measurements. Furthermore, these methods were demonstrated with only regular phantoms or homogeneous phantom simulations and presented no in vivo experiment validation.

In this paper, we propose an incomplete variables truncated conjugate gradient (IVTCG) algorithm for whole domain BLT reconstruction. In view of the characteristics of bioluminescent source distribution in most BLT applications, the solution should be in a sparse format where most of its elements are zeros. We first reformulate the BLT inverse problem into minimizing an objective function which integrates a sparseness-inducing (ℓ ₁ norm) regularization term with a quadratic error term. The IVTCG algorithm is then applied to solve this ℓ ₁-norm minimization problem.

Numerical simulations on the mouse atlas model and in vivo mouse experiments validate the proposed scheme and the performance is assessed in terms of location accuracy, relative error of reconstructed power, reconstruction time, robustness to noise, and robustness to optical parameters perturbation.

The paper is organized as follows. In section 2, we present the DA model and its finite element solution. Section 3 elaborates on the IVTCG based sparse reconstruction method for BLT. We then validate the proposed method with a mouse atlas in section 4 and show an in vivo experiment and results in section 5. Finally, we provide the discussion and conclusion.

2. Bioluminescence tomography framework

2.1 Photon propagation model

BLT reconstruction must take into account photon propagation through the underlying biological tissues. The RTE is a precise mathematical model describing light transport in turbid media [20,21,24,25], but a high computation burden limits its application in practice. Under the assumption that scattering dominates absorption for near-infrared light in biological tissues, the RTE can be simplified by diffusion approximation [26]. In this paper, we model the photon propagation with the steady diffusion equation complemented by a Robin boundary condition [11,24–27]:

- \nabla \cdot (D (r) \nabla Φ (r)) + μ_{a} (r) Φ (r) = S (r) (r \in Ω),

Φ (r) + 2 A (r; n, n^{'}) D (r) (v (r) \cdot \nabla Φ (r)) = 0 (r \in \partial Ω),

where

r \in R^{3}

is the position vector in domain Ω,

S (r)

represents the power density of the internal bioluminescence source,

Φ (r)

denotes the photon fluence rate, ν denotes the unit outer normal at boundary

\partial Ω

,

A (r; n, n^{'})

is the boundary mismatch factor [25],

D (r) = 1 / (3 (μ_{a} (r) + μ_{s}^{'} (r)))

is the optical diffusion coefficient with

μ_{a} (r)

being the optical absorption coefficient and

μ_{s}^{'} (r)

is the reduced scattering coefficient respectively. In a bioluminescent imaging experiment, the measured quantity is the outgoing photon density on ∂Ω [5]:

Q (r) = - D (r) (v (r) \cdot \nabla Φ (r)) = {(2 A (r; n, n^{'}))}^{- 1} Φ (r) (r \in \partial Ω) .

2.2 BLT inverse model

Given the optical properties of the underlying medium, the linear relationship between the interior source density and the boundary photon flux rate are derived by solving Eq. (2) and (3) with the finite element method [11,28], and the matrix equation that connects the discretized fluence rate Φ and the discretized source distribution S can be expressed as:

M Φ = F S,

where M is a positive definite matrix and F is the source weight matrix. Thus, the photon fluence rateΦis derived by

Φ = M^{- 1} F S = \bar{A} S .

Since only partial photons on the boundary are captured by CCD in BLT experiments by retaining those rows associated with the boundary measurements $Φ^{m}$ in the coefficient matrix $\bar{A}$ , $\bar{A}$ becomes A, and then we obtain the inverse model of BLT:

A S = Φ^{m},

where

A \in R^{M \times N} (M < N)

. It is known that this linear system of equations is under-determined and typically ill-conditioned.

3. Sparse reconstruction method

3.1 Sparse regularization

In view of the insufficient measurement and sparse characteristic of the light source in BLT applications, source reconstruction is equivalent to finding a sparse solution to Eq. (6). Inspired by the CS theory, we convert the BLT inverse model to the following minimization problem with ℓ ₁ regularization:

\min_{S} \frac{1}{2} {‖ A S - Φ^{m} ‖}_{2}^{2} + τ {‖ S ‖}_{1},

where

S \in R^{N}

,

Φ^{m} \in R^{M}

,

A \in R^{M \times N} (M < N)

,

τ > 0

is a regularization parameter,

{‖ \cdot ‖}_{2}

denotes the Euclidean norm, and

{‖ S ‖}_{1} = \sum_{i} | s_{i} |

is the ℓ ₁ norm.

3.2 IVTCG method

Since the objective function in Eq. (7) is convex but not differentiable, we reformulate it as a convex quadratic program with nonnegative constrained conditions by a similar method used in gradient projection for sparse reconstruction (GPSR) [29], an efficient algorithm for large scale sparse problems in compressed sensing and other inverse problems in signal processing and statistics. Specifically, we introduce vectors u and v, and make the substitution $S = u - v,$ $u \geq 0,$ $v \geq 0$ . Hence, Eq. (7) can be rewritten as the following program:

\begin{array}{l} \min_{z} c^{T} z + \frac{1}{2} z^{T} B z \equiv F (z), \\ s . t . z \geq 0, \end{array}

where

z = {[\begin{matrix} u & v \end{matrix}]}^{T}

,

c = τ 1_{2 N} + {[\begin{matrix} - b & b \end{matrix}]}^{T}

,

1_{2 N} = {[1, 1, \dots, 1]}^{T}

,

b = A^{T} Φ^{m}

, and

B = [\begin{matrix} A^{T} A & - A^{T} A \\ - A^{T} A & A^{T} A \end{matrix}] = {[A, - A]}^{T} [A, - A]

. The equivalence of the problems (8) and (7) has been proven in [29].

Assuming $z * = {[{(u *)}^{T}, {(v *)}^{T}]}^{T}$ is the optimal point of programming (8), since bioluminescence source distribution is sparse in space, $z *$ should be in a sparse format where most of its elements are zeros. Let $N_{s} > 0$ be a constant, and we assume that the number of the nonzero (NNZ) entries of $z *$ is less than $N_{s}$ .

Let the index set of nonzero components in $z *$ be $I * = u * \cup v *,$ where $u * ={i \in {1, \dots, N} | u_{i}^{*} > 0}$ , $v * = {i \in {N + 1, \dots, 2 N} | v_{i}^{*} > 0}$ . If we obtain $I *$ , programming (8) can be solved by optimizing a small sub-problem only involving variables in $I *$ .

In this paper, we design an iterative algorithm named IVTCG to seek a better $I^{k}$ as an approximation of $I^{*}$ in each iteration, and then optimize the sub-problem formed by setting the other variables to zero.

Unlike most optimization methods which update the entire solution vector in each step of the iterative process, the IVTCG algorithm modifies only partial variables per iteration. Consequently, we need to select a working set first.

According to the optimization theory, Karush-Kuhn-Tucker (KKT) conditions are necessary and sufficient for optimality. The KKT system of Eq. (8) is

\nabla F (z *) = Ι υ * = υ *,

υ * \geq 0, z * \geq 0, {(υ *)}^{T} z * = 0,

where

υ *

is a length-2N vector constituted by KKT multipliers.

If for all $i \in {1, \dots 2 N}$ there is one and only one of ${(υ *)}_{i}$ and ${(z *)}_{i}$ , which is 0, then programming (8) satisfies the strict complementary slackness at $z *$ . It follows that the KKT-system is equivalent to

\nabla F (z *) \geq 0, z * \geq 0, {(\nabla F (z *))}^{T} z * = 0, \nabla F (z *) + z * > 0.

Following a similar method used in [30] and [31], a working set $Γ^{k}$ can be derived by the violation of the KKT optimal conditions at the k-th iteration.

Γ^{k} = {i | i \in {1, \dots, 2 N}, [{(z^{k})}_{i} > 0, {(\nabla F (z^{k}))}_{i} \neq 0] o r [{(z^{k})}_{i} = 0, {(\nabla F (z^{k}))}_{i} < 0]},

where

z^{k}

is the current iteration point.

It is easy to see that condition (11) is equivalent to $\min {z *, \nabla F (z *)} = 0$ .Consequently, if $w ≜ \min {z *, \nabla F (z *)}$ , ${‖ w ‖}_{2} = 0$ is another necessary and sufficient condition for optimality. Moreover, since the strict complementary slackness holds, when $z^{k} \to z *$ it follows that,

\frac{{(z^{k})}_{i}}{| {(\nabla F (z *))}_{i} |} \to + \infty \forall i \in I *,

\frac{{(z^{k})}_{j}}{| {(\nabla F (z *))}_{j} |} \to 0 \forall j \notin I * .

Equation (13) and (14) can help us to form the $I^{k}$ . First, we define an index set

{\hat{I}}^{_{k}} = {i \in {1, \dots, 2 N} | {(z^{k})}_{i} > 0, {(z^{k})}_{i} / {(\nabla F (z^{k}))}_{i} > δ}

where

δ > 0

is a constant. Then, we sort

{(z^{k})}_{i} / {(\nabla F (z^{k}))}_{i}

in descending order as

\frac{{(z^{k})}_{i 1}}{| {(\nabla F (z *))}_{i 1} |} \geq \frac{{(z^{k})}_{i 2}}{| {(\nabla F (z *))}_{i 2} |} \geq \dots \geq δ .

The index set $I^{k}$ is obtained

I^{k} = {i_{l} \in {\hat{I}}^{k} | l \leq \min {| {\hat{I}}^{k} |, N_{s}}},

where

| {\hat{I}}^{k} |

denotes the cardinality of

{\hat{I}}^{k}

. It is obvious that the constant δ and

N_{s}

control the size of

I^{k}

.

Let ${\hat{J}}^{k} = Γ^{k} \ I^{k}$ . For the other variables corresponding to ${\hat{J}}^{k}$ , we sort $| {(\nabla F (z *))}_{j} |$ in descending order and define it as:

J^{k} = {j_{l} \in {\hat{J}}^{k} | l \leq \min {| {\hat{J}}^{k} |, N_{\max} - N_{s}}},

where

N_{\max} > N_{s}

is a constant that controls the size of

J^{k}

.

The variables related to $I^{k}$ and $J^{k}$ will be updated at the current iteration. However, different search directions and methods are used for $I^{k}$ and $J^{k}$ . Specifically, for $I^{k}$ the search direction $d_{I^{k}}^{k}$ is determined by solving a sub-problem that optimizes only the variables corresponding to $I^{k}$ , i.e.

\begin{array}{l} \min_{d_{I_{k}}^{k}} F (z_{I^{k}}^{k} + d_{I^{k}}^{k}, z_{{\bar{I}}^{k}}^{k} + d_{{\bar{I}}^{k}_{k}}^{k}), \\ s . t . z_{I^{k}}^{k} + d_{I^{k}}^{k} \geq 0, \end{array}

where

{\bar{I}}^{k} = {1, \dots, 2 N} \ I^{k}

. By omitting the constant in Eq. (16), the sub-problem can be rewritten as

\begin{array}{l} \min_{x} b_{s u b}^{T} x + \frac{1}{2} x^{T} B_{s u b} x \equiv F_{s u b} (x), \\ s . t . z_{I^{k}}^{k} + x \geq 0, \end{array}

where

b_{s u b}^{T} = {(\nabla F (z^{k}))}_{I^{k}}

and

B_{s u b} = B_{I^{k} I^{k}}

. For

J^{k}

we define the direction

d_{J^{k}}^{k} = - w_{J^{k}}

. The other entries of

d^{k}

are set to 0. After synthesizing the descent direction vector with

d_{I^{k}}^{k}

,

d_{J^{k}}^{k}

and 0, the next iteration point can be obtained by backtracking according to the Armijo rule [30].

The main algorithm of IVTCG is summarized in Algorithm1.

Algorithm 1: main algorithm of IVTCG

Step 0 (initialization): Given a length-2N initial vector

z^{0} = 0

, choose parameters

δ > 0

,

0 < β < 1 / 2,

0 < γ < 1, ε \geq 0, N_{s} > 0, N_{\max} > N_{s}

, and set

k = 0

. Step 1: Calculate

w^{k} = \min {z^{k}, \nabla F (z^{k})}

. Step 2: If

{‖ w^{k} ‖}_{2} \leq ε

, stop; otherwise, go to Step 3. Step 3: Update

I^{k}

and

J^{k}

according to Eq. (16) and Eq. (17). Step 4: Calculate

d_{I_{k}}^{k}

by solving the sub-problem of Eq. (19) with the truncated conjugate gradient method and synthesize in a descending direction

d^{k} = {[{(d_{I_{k}}^{k})}^{T}, {(- w_{J_{k}})}^{T}, 0]}^{T}

. Step 5: Calculate the minimum nonnegative integer q such that

F (z^{k} + γ^{q} d^{k}) \leq F (z^{k}) + β γ^{q} {(g^{k})}^{T} d^{k}

and let

z^{k + 1} = z^{k} + γ^{q} d^{k}, k = k + 1,

go to Step 1.

The truncated conjugate gradient method used for solving the sub-problem in Eq. (19) is summarized in Algorithm 2 [32].

Algorithm 2: the truncated conjugate gradient method for the sub-problem

Step 0: Let

x^{0} = 0

,

\nabla F_{s u b} (x^{0}) = b_{s u b}

,

d_{s u b}^{0} = - \nabla F_{s u b} (x^{0})

,

α_{\max} > 0

,

ε_{s u b} \geq 0

, and

i t e r_{\max} > 0

; set

t = 0

. Step 1: If

{‖ \nabla F_{s u b} (x^{t}) ‖}_{2}^{2} \leq ε_{s u b}

or

t > i t e r_{\max}

, then go to step 4; otherwise, let

α_{t} = {\begin{matrix} α_{\max}, & i f α_{\max} {(d_{s u b}^{t})}^{T} B_{s u b} d_{s u b}^{t} \leq {‖ \nabla F_{s u b} (x^{t}) ‖}_{2}^{2} \\ {‖ \nabla F_{s u b} (x^{t}) ‖}_{2}^{2} / ({(d_{s u b}^{t})}^{T} B_{s u b} d_{s u b}^{t}) & i f α_{\max} {(d_{s u b}^{t})}^{T} B_{s u b} d_{s u b}^{t} > {‖ \nabla F_{s u b} (x^{t}) ‖}_{2}^{2} \end{matrix}

. Step 2: If

z_{I_{k}}^{k} + s^{t} + α_{t} d_{s u b}^{t} \geq 0

, then set

x^{t + 1} = x^{t} + α_{t} d_{s u b}^{t}

,

\nabla F_{s u b} (x^{t + 1}) = \nabla F_{s u b} (x^{t}) + α_{t} B_{s u b} d_{s u b}^{t}

,

ζ_{t} = \frac{{‖ \nabla F_{s u b} (x^{t + 1}) ‖}_{2}^{2}}{{‖ \nabla F_{s u b} (x^{t}) ‖}_{2}^{2}}

,

d_{s u b}^{t + 1} = - \nabla F_{s u b} (x^{t + 1}) + ζ_{t} d_{s u b}^{t}

,

t = t + 1

, and go to step 1; otherwise, go to step 3. Step 3: Calculate the maximum nonnegative scalar

α^{*}

such that

z_{I_{k}}^{k} + x^{t} + α_{t} d_{s u b}^{t} \geq 0

. Let

x^{t} = x^{t} + α^{*} d_{s u b}^{t}

and go to step 4. Step 4:

d_{I_{k}}^{k} = x^{t}

and stop.

For all of the experiments in this paper, IVTCG adopts the following parameter settings: $N_{s} = ⌊ M / 10 ⌋$ , $N_{\max} = N_{s} + ⌊ N_{s} / 8 ⌋$ , $δ = 7, β = 0.01, γ = 0.9, α_{\max} = 1 e 10, ε_{s u b} = 1 e - 10$ , and $i t e r_{\max} = N_{s}$ .

We can make a concise analysis of the time complexity of the IVTCG algorithm as follows. For the gradient $\nabla F (z^{k}) = c + B z^{k} = \nabla F (z^{k - 1}) + γ^{q} B d^{k - 1}$ , the multiplication operations focus on $B d^{k - 1} = {[A, - A]}^{T} [A, - A] d^{k - 1}$ . Since the NNZ entries of $d^{k - 1}$ are less than $N_{\max}$ , it need $M \times N_{\max}$ multiplications to compute $A d^{k - 1}$ .The amount of multiplication $A^{T} A d^{k - 1}$ is $M \times N$ . Therefore, the total cost of $g^{k}$ is $O (M N_{\max}) + O (M N) = O (M N)$ . As for the back-tracking step, we find that the cost is $O (M N_{\max})$ . In addition, suppose m is the size of the sub-problem (i.e. $m = | I_{k} |$ , and m<N_s), since it needs at most m iterations for solving the sub-problem and the cost per iteration is $O (M m)$ , the computational cost regarding the sub-problem is $O (M m^{2})$ . Thus, the overall cost for each iterative step of IVTCG is $O (M m^{2}) + O (M N)$ .

From the above analysis, one can find that parameter N_s controls the size of the sub-problem and the computational cost may increase sharply for a non-very sparse problem. In fact, by adjusting the parameters of IVTCG with slight modification, it can also be applied to other inverse problems in which this type of ℓ ₁-regularized least squares minimization appears, which we will discuss in another paper [33].

4. Simulation studies in the mouse atlas

In our numerical simulations, a mouse atlas of CT and cryosection data was employed to provide anatomical information [34], and only the torso section with a height of 40 mm was selected as the volume to be investigated. The optical properties of this model are listed in Table 1 [35–37].

Table 1. Optical properties for the mouse organs

View Table | View all tables in this article

4.1 Quantitative reconstruction in a single source

In the first set of experiments, a cylindrical source with a 0.5 mm radius and 1 mm height was placed in the left kidney with the center at (10.2 mm, 15.5 mm, 24.5 mm) as shown in Fig. 1(a) .

Fig. 1 Reconstruction model in a single source case. (a) The torso of the mouse atlas model with a cylindrical source in the left kidney, (b) the mesh for reconstruction and the simulated photon distribution on the surface.

Download Full Size | PDF

The synthetic measurements on the boundary were generated by using the FEM to solve the forward model [11]. Specifically, the atlas model with the interior bioluminescent source was discretized into a tetrahedral-element mesh consisting of 142920 elements and 26587 nodes. In addition, Lagrange quadratic shape functions were adopted to make the surface flux density more accurate. Figure 1(b) presents the simulated photon distribution on the boundary and the mesh used for reconstruction. The actual source density was set at 1 nanoWatts/mm³, the actual source volume was 0.5404 mm³ after tetrahedralization with FEM, and thus we could figure out the initial source power was 0.5404 nanoWatts. The reconstruction mesh consisted of 14776 elements and 3098 nodes with 1155 nodes on the surface; hence, the size of the system matrix was 1155 × 3098.

Previous studies in [19–22] demonstrated the superiority of ℓ ₁-norm regularization to ℓ ₂-norm for sparse reconstruction in BLT, thus no such comparisons were presented in this paper. We compared the IVTCG method with ℓ1-ℓs and GPSR to exhibit the advantages of the proposed scheme, the former has been used in BLT [21,22] and the latter is a gradient projection algorithm that is very efficient for large scale sparse inverse problems [29]. All of the algorithms were coded in Matlab^TM and carried out on a personal computer with 2.6 GHz Intel^® Pentium^TM E5300 processor and 2GB RAM.

In order to quantitatively evaluate the reconstruction quality, besides the source location error (LE) and reconstructed density, the reconstructed source power was also considered. This is because reconstructed total power is very important for practical BLT, for example, it can reflect the total tumor cell number, which maybe crucial for longitudinal monitoring using BLT. Additionally, the reconstructed density is usually related to mesh dimension and it is difficult to differentiate the influence of a small source with high density and a large one with low density [8]. In this paper, LE was the Euclidean distance between centers of the reconstructed and the actual source; the source power was computed with the source integral method [11,15] and the relative error (RE) of power was calculated by $| {Power}_{r} - {Power}_{a} | / {Power}_{a}$ , where ${Power}_{a}$ denoted the actual value and ${Power}_{r}$ was the reconstructed value.

Because the regularization parameter plays an important role in regularization methods, we performed three groups of comparison experiments with different regularization parameters using the three methods. The detailed information about parameters and the final quantitative reconstruction results are summarized in Table 2 .

Table 2. Reconstruction results in single source case with different parameters

View Table | View all tables in this article

We found that the reconstructed positions by IVTCG and ℓ1-ℓs were identical in all cases considered in this set of experiments. Specifically, the reconstructed center was (9.91 mm, 15.87 mm, 24.45 mm) with an LE of 0.47 mm from the actual source, whereas the LE by GPSR was up to 2.23 mm. For each τ, GPSR performed faster than IVTCG and ℓ1-ℓs, but the quantitative results were inferior to that of the other two methods in terms of location, density and power. As for IVTCG, it performed slightly slower than GPSR but much faster than ℓ1-ℓs; for example, for $τ = 1e - 5$ although ℓ1-ℓs and IVTCG produced comparative quantitative results, the reconstruction time of ℓ1-ℓs was about 319 times of that of IVTCG.

The reconstruction results in Table 2 demonstrate that the regularization parameter does affect the final reconstruction quality. Although the reconstructed positions do not vary with τ, the quantitative values of density and power descend when τdecreases. The figures in Table 2 witness an increase in the relative error of reconstructed power by IVTCG i.e., from 15.80% at $τ = 1e - 3$ to 44.1% at $τ = 1e - 5$ . In this set of experiments, $τ = 1e - 3$ was the optimal regularization parameter and the final reconstruction results in this case are shown in Fig. 2 . Although the quantitative results by ℓ1-ℓs and IVTCG are comparative, there are some artifacts in the results of ℓ1-ℓs as shown in Fig. 2(d), which indicates the solution is not sparse enough. In this case, IVTCG produced a very sparse solution with only 3 nonzero entries, whereas the NNZ of the counterparts were 7 for GPRS and 15 for ℓ1-ℓs respectively.

Fig. 2 Reconstruction results at $τ = 1e - 3$ . (a), (c) and (e) are respectively the isosurface views of the results by IVTCG, ℓ1-ℓs and GPSR, where the small red cylinder denotes the actual source; (b), (d) and (f) are the corresponding transverse views at z = 24.5mm where the small black circle denotes the actual source.

Download Full Size | PDF

4.2 Reconstruction robustness to measurement noise

It is well known that the system matrix A is typically ill-conditioned which make the BLT reconstruction very sensitive to measurement noise. A set of simulations were performed to test the reconstruction robustness by adding different levels of noise to boundary measurements. We conducted 100 independent reconstructions for each noise level. For all of the noise levels considered, the source locations were identical to those without noise. Furthermore, we found that the reconstructed power varied slightly with the increase in noise level, and the maximum deviation of the power occurred at 30% noise level, which possessed a maximum 9.4% deviation to the actual power. The error bar chart of the reconstructed power under different noise level is shown in Fig. 3 . The results show that the proposed method is robust to measurement noise.

Fig. 3 Error bar chart of power under different noise levels.

Download Full Size | PDF

4.3 Reconstruction robustness to optical parameters perturbation

In most BLT applications, there is inevitable discrepancy between the optical properties used in the reconstruction and the actual physiological values which would affect the reconstruction accuracy significantly [27,38]. In this section, we investigated the robustness of the proposed method against inaccurate optical parameters by adding $\pm 10 %$ perturbation to absorption and reduced scattering coefficients of different organs. The simulations were performed on the same model in the single source case. We considered all of the 9 cases of the combination of different deviations in $μ_{a}$ and $μ_{s}^{'}$ .

The reconstructed source center was not affected by the optical parameter perturbation. The results in Table 3 show that the proposed method is robust against optical parameter perturbation especially in source location. We can observe an interesting phenomenon that the reconstructed results in case 1 to case 4 are superior to those without any perturbation. It is noted that the common feature of the four cases is the ratio of the reduced scattering coefficient to the absorption coefficient increase, which means such cases better satisfy the hypothesis of the DA model, i.e. scattering predominates over absorption.

Table 3. Reconstruction results with optical parameter perturbation

View Table | View all tables in this article

4.4 Double source case evaluation

We also investigated the proposed method by a double source experiment. Two sources were located in the left kidney with their centers at (10.2 mm, 15.5 mm, 24.5 mm) and (10.2 mm, 16 mm, 28.5 mm) respectively. The size and density of each source were the same as in the single source case, but the initial power of source-1 was 0.5404 nanoWatts, whereas that of source-2 was 0.5378 nanoWatts, mainly due to the influence of the mesh. The final reconstruction results presented in Table 4 and Fig. 4 indicate that the sources can be accurately distinguished, although there is some difference in the reconstructed density and power.

Table 4. Reconstruction results in the double-source case

View Table | View all tables in this article

Fig. 4 Reconstruction results in the double source case (a) The isosurface view (b) The corresponding coronal view at Y = 15.8mm. The small red cylinders in (a) and small black squares in (b) denote actual sources.

Download Full Size | PDF

5. In vivo experiment validation

To further test the proposed method, an in vivo experiment was performed on an adult BALB/C mouse. The animal procedures were in accordance with the Fourth Military Medical University (FMMU) approved animal protocol.

A luminescent catheter filled with 5μl luminescent liquid was implanted into the abdomen to serve as the testing source in this experiment. The luminescent solution was extracted from a red luminescent light stick (Glow products, Canada) and the generated luminescent light had an emission peak wavelength of about 644 nm. The initial total power was 300 nanoWatts (the total power = luminescent solution volume × luminescent solution flux density = 5μl × 60 nanoWatts/μl).

The experimental data were acquired by a dual-modality BLT/micro-CT system developed in our lab [39]. The anesthetized mouse was first photographed and luminescent images were taken by a calibrated CCD camera from four directions at 90 degree intervals with different exposure times. Specifically, the luminescence images were acquired with a binning value of 2, integration time of 12 seconds and an aperture number f 8. The CCD camera was calibrated by an integrating sphere 12 inches in diameter (USS-1200V-LL Low-Light Uniform Source, Labsphere, North Sutton, NH) and the calibration formula is the same as used in [39]. The multi-view superimposed photographs and luminescent images are shown in Figs. 5(a)-(d) . After the optical data were acquired, the intact mouse was scanned using the Micro-CT. The volume data were reconstructed by GPU-accelerated FDK algorithm [40]. Because of the limited field of view, only the torso section was scanned. The center coordinates (22.08 mm, 22.24 mm, 18.08 mm) of the actual luminescent source were obtained from the CT slices.

Fig. 5 Multi-view superimposed images of photographs and luminescent images. (a)-(d) Anterior -posterior, right-lateral, posterior-anterior, and left-lateral views respectively.

Download Full Size | PDF

For implement reconstruction, the CT slices were segmented into major anatomical components, including lungs, heart, liver, kidneys, and muscle as shown in Fig. 6(a) . The optical parameters for different organs were calculated based on the literature [36] as listed in Table 5 .

Fig. 6 In vivo model (a) The 3D view of the segmented micro-CT slices of the imaged mouse with a luminescent source implanted beneath the liver. (b) Surface view of the reconstruction mesh, with the detected photon distribution mapping on it.

Download Full Size | PDF

Table 5. Optical properties of the mouse model

View Table | View all tables in this article

Note that bioluminescent images did not contain any structural information, thus we needed to match optical data to the coordinate system of the volume data. In the experiment, each optical data acquisition consisted of a visible-light photograph and a corresponding luminescent emission image. Because these photographs were acquired with the same camera, these images were expressed in the same coordinate space, and could therefore be used for registering the optical data with the CT data.

By using the rotation stage and the mouse holder, there was no substantial deformation of the mouse during the entire process. Therefore, a landmarks-based rigid-body registration method was adopted in this work. Four identical polyethylene balls placed on every side of the mouse holder served as landmarks for subsequent registration. The positions of the landmarks were read from the volume data. Meanwhile, 2D coordinates of the same marks on the planar optical photographs could be obtained. In addition, the third coordinate could be calculated through the optical imaging principle with the basic parameters, including the perpendicular distance between the image and the lens system, the distance between the object and the lens system, the focus of the lens system, and the magnification of the optical system. After registering the bioluminescence images and the volume data of the micro-CT, the absolute irradiance distribution was mapped onto the mouse surface using a 3D surface flux reconstruction algorithm based on the hybrid radiosity-radinace theorem as shown in Fig. 6(b) [41].

In the reconstruction process, the segmented mouse torso was discretized into a tetrahedral-element mesh containing 11085 elements and 2390 nodes for subsequent reconstruction. It took about 6 seconds to complete the reconstruction using the IVTCG method without any a priori information of the permissible source region. The final result is presented in Fig. 7 , where the source center is (23.58 mm, 24.17 mm, 18.05 mm) with a deviation of 2.44 mm to the actual center. The reconstructed power was 276.7 nanoWatts possessing a relative error of 7.8%.

Fig. 7 The reconstruction result (a) the isosurface view of the reconstruction result. (b) the transverse view of the result and the comparison with the corresponding CT slices. The cross of the green lines denotes the actual source center and the cross of the red lines denotes the reconstructed center.

Download Full Size | PDF

6. Conclusion and discussion

In this paper, we present a novel BLT sparse reconstruction method without any a priori information of the permissible source region. Based on the simulation experiments on a mouse atlas model we demonstrate that the proposed reconstruction method is able to accurately localize and quantify light source distribution of the entire region even with noisy measurements and inaccurate optical parameters. The in vivo experiment further indicates its capability of quantitative reconstruction on whole bodies of mice.

It is well known that the ill-posedness of the reconstruction problem and large scale numerical problems involved in tomography are two vital issues that reconstruction algorithms need to address. Our simulations indicate that sparseness-inducing (ℓ ₁ norm) regularization can effectively reduce the ill-posedness of BLT and thus produce a very sparse and stable solution without permissible region constraint. On the other hand, by limiting the size of the variables updated in the current iteration and locating the nonzero entries of the sparse solution more efficiently, the IVTCG based reconstruction algorithm can decrease the computational costs and converge after finite iterations, which make it more appropriate for practical BLT applications.

It can be noted that in vivo experiments are not as accurate as simulation ones, but they are reasonable for such a newly developed imaging system. The reconstruction quality could be enhanced by improving the imaging system and the experimental procedures. For the location error, one of the error sources may originate from 3D surface mapping. In addition, although the reconstruction is based on a heterogeneous model, some organs are ignored for simplicity which would also lead to errors.

Since bioluminescence signals are generally very weak and the designation of volume of interest needs human intervention, reconstruction algorithms that need no explicit knowledge of the permissible region are thus very meaningful. Recently, a generalized graph cuts based approach has been proposed for localizing the bioluminescent source in the entire region [42]. Nevertheless, it cannot quantify the source distribution, and that is why we did not compare it with our method.

Mathematically, the goal of bioluminescence source reconstruction can be considered to find sparse solutions from large underdetermined systems of linear equations, which is essentially the same as that of compressive sensing tasks. Although verifying the restricted isometry property (RIP) condition or incoherence condition of the CS theory is very difficult [43], the experiment results in this paper demonstrate that the CS techniques enhance the numerical stability of the BLT reconstruction.

In general, the proposed scheme is an efficient quantitative reconstructions method for complex BLT applications. Although it can recover the interior source efficiently, the error caused by the limitation of the DA model is still inevitable, more precise models should be considered for the practical BLT problem. In the future, we will further investigate this method for BLT reconstruction based on a more accurate forward model.

Acknowledgments

This work is supported by the Program of the National Basic Research and Development Program of China (973) under Grant Nos. 2006CB705700, 2011CB707702, CAS Hundred Talents Program, the National Natural Science Foundation of China under Grant Nos. 81090272, 81000632, 30900334, the Shaanxi Provincial Natural Science Foundation Research Project under Grant No. 2009JQ8018, the Fundamental Research Funds for the Central Universities, and the Science Foundation of Northwest University under Grant No. 09NW34.

References and links

1. V. Ntziachristos, J. Ripoll, L. V. Wang, and R. Weissleder, “Looking and listening to light: the evolution of whole-body photonic imaging,” Nat. Biotechnol. 23(3), 313–320 (2005). [CrossRef] [PubMed]

2. C. H. Contag and M. H. Bachmann, “Advances in in vivo bioluminescence imaging of gene expression,” Annu. Rev. Biomed. Eng. 4(1), 235–260 (2002). [CrossRef] [PubMed]

3. J. K. Willmann, N. van Bruggen, L. M. Dinkelborg, and S. S. Gambhir, “Molecular imaging in drug development,” Nat. Rev. Drug Discov. 7(7), 591–607 (2008). [CrossRef] [PubMed]

4. R. Weissleder and V. Ntziachristos, “Shedding light onto live molecular targets,” Nat. Med. 9(1), 123–128 (2003). [CrossRef] [PubMed]

5. 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,” Curr. Med. Imaging Rev. 2(4), 453–457 (2006). [CrossRef]

6. 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]

7. 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]

8. H. Dehghani, S. C. Davis, and B. W. Pogue, “Spectrally resolved bioluminescence tomography using the reciprocity approach,” Med. Phys. 35(11), 4863–4871 (2008). [CrossRef] [PubMed]

9. G. Wang, H. Shen, W. Cong, S. Zhao, and G. W. Wei, “Temperature-modulated bioluminescence tomography,” Opt. Express 14(17), 7852–7871 (2006), http://www.opticsinfobase.org/oe/viewmedia.cfm?URI=oe-14-17-7852&seq=0. [CrossRef] [PubMed]

10. S. Ahn, A. J. Chaudhari, F. Darvas, C. A. Bouman, and R. M. Leahy, “Fast iterative image reconstruction methods for fully 3D multispectral bioluminescence tomography,” Phys. Med. Biol. 53(14), 3921–3942 (2008). [CrossRef] [PubMed]

11. W. Cong, G. Wang, D. Kumar, Y. Liu, M. Jiang, L. V. Wang, E. A. Hoffman, G. McLennan, P. B. McCray, J. Zabner, and A. Cong, “Practical reconstruction method for bioluminescence tomography,” Opt. Express 13(18), 6756–6771 (2005), http://www.opticsinfobase.org/jdt/viewmedia.cfm?id=85344&seq=0. [CrossRef] [PubMed]

12. X. He, J. Liang, X. Qu, H. Huang, Y. Hou, and J. Tian, “Truncated total least squares method with a practical truncation parameter choice scheme for bioluminescence tomography inverse problem,” Int. J. Biomed. Imaging 2010, 291874 (2010). [CrossRef] [PubMed]

13. Y. Lv, J. Tian, W. Cong, G. Wang, J. Luo, W. Yang, and H. Li, “A multilevel adaptive finite element algorithm for bioluminescence tomography,” Opt. Express 14(18), 8211–8223 (2006), http://www.opticsinfobase.org/jot/viewmedia.cfm?id=97939&seq=0. [CrossRef] [PubMed]

14. J. Feng, K. Jia, G. Yan, S. Zhu, C. Qin, Y. Lv, and J. Tian, “An optimal permissible source region strategy for multispectral bioluminescence tomography,” Opt. Express 16(20), 15640–15654 (2008), http://www.opticsinfobase.org/oe/viewmedia.cfm?uri=oe-16-20-15640&seq=0. [CrossRef] [PubMed]

15. H. Huang, X. Qu, J. Liang, X. He, X. Chen, D. Yang, and J. Tian, “A multi-phase level set framework for source reconstruction in bioluminescence tomography,” J. Comput. Phys. 229(13), 5246–5256 (2010). [CrossRef]

16. M. A. Naser and M. S. Patterson, “Algorithms for bioluminescence tomography incorporating anatomical information and reconstruction of tissue optical properties,” J. Biomed. Opt. Express 1(2), 512–526 (2010). [CrossRef]

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

18. E. Candès, “Compressive sampling,” in Proceedings of the International Congress of Mathematicians, (Madrid, Spain, 2006), pp. 1433–1452.

19. 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), http://www.opticsinfobase.org/abstract.cfm?URI=oe-17-10-8062. [CrossRef] [PubMed]

20. W. Cong and G. Wang, “Bioluminescence tomography based on the phase approximation model,” J. Opt. Soc. Am. A 27(2), 174–179 (2010). [CrossRef]

21. H. Gao and H. Zhao, “Multilevel bioluminescence tomography based on radiative transfer equation Part 1: l1 regularization,” Opt. Express 18(3), 1854–1871 (2010), http://www.opticsinfobase.org/oe/viewmedia.cfm?URI=oe-18-3-1854&seq=0. [CrossRef] [PubMed]

22. X. He, Y. Hou, D. Chen, Y. Jiang, M. Shen, J. Liu, Q. Zhang, and J. Tian, “Sparse regularization-based reconstruction for bioluminescence tomography using a multilevel adaptive finite element method,” Int. J. Biomed. Imaging 2011, 203537 (2011). [CrossRef]

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

24. 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]

25. L. V. Wang, and H. Wu, Biomedical Optics: Principles and Imaging (New York: John Wiley, 2007), chap.5.

26. H. Dehghani, D. T. Delpy, and S. R. Arridge, “Photon migration in non-scattering tissue and the effects on image reconstruction,” Phys. Med. Biol. 44(12), 2897–2906 (1999). [CrossRef]

27. Q. Z. Zhang, L. Yin, Y. Y. Tan, Z. Yuan, and H. B. Jiang, “Quantitative bioluminescence tomography guided by diffuse optical tomography,” Opt. Express 16(3), 1481–1486 (2008), http://www.opticsinfobase.org/josab/viewmedia.cfm?id=149907&seq=0. [CrossRef] [PubMed]

28. M. Schweiger, S. R. Arridge, M. Hiraoka, and D. T. Delpy, “The finite element method for the propagation of light in scattering media: boundary and source conditions,” Med. Phys. 22(11 Pt 1), 1779–1792 (1995). [CrossRef] [PubMed]

29. M. Figueiredo, R. Nowak, and S. Wright, “Gradient projection for sparse reconstruction: Application to compressed sensing and other inverse problems, IEEE J. Select,” IEEE J. Sel. Top. Signal Process. 1(4), 586–597 (2007). [CrossRef]

30. E. Osuna, R. Freund, and F. Girosi, “An improved algorithm for support vector machines,” in Proceedings of IEEE Signal Processing Society Workshop (Amelia Island, USA, 1997), pp. 276–285.

31. R. Fan, P. Chen, and C. Lin, “Working Set Selection Using Second Order Information for Training Support Vector Machines,” J. Mach. Learn. Res. 6, 1889–1918 (2005).

32. H. Qi, L. Qi, and D. Sun, “Soving Karush-Kuhn-Tucker systems via trust region and the conjugate gradient methods,” SIAM J. Optim. 14(2), 439–463 (2003). [CrossRef]

33. X. Wang and F. Liu “Incomplete variables truncated conjugate gradient method for signal reconstruction in compressed sensing,” in preparation.

34. B. Dogdas, D. Stout, A. F. Chatziioannou, and R. M. Leahy, “Digimouse: a 3D whole body mouse atlas from CT and cryosection data,” Phys. Med. Biol. 52(3), 577–587 (2007). [CrossRef] [PubMed]

35. 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), http://www.opticsinfobase.org/as/viewmedia.cfm?id=97670&seq=0. [CrossRef] [PubMed]

36. 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]

37. X. Song, D. Wang, N. Chen, J. Bai, and H. Wang, “Reconstruction for free-space fluorescence tomography using a novel hybrid adaptive finite element algorithm,” Opt. Express 15(26), 18300–18317 (2007), http://www.opticsinfobase.org/VJBO/viewmedia.cfm?uri=oe-15-26-18300&seq=0. [CrossRef] [PubMed]

38. A. Cong, W. Cong, Y. Lu, P. Santago, A. Chatziioannou, and G. Wang, “Differential evolution approach for regularized bioluminescence tomography,” IEEE Trans. Biomed. Eng. 57(9), 2229–2238 (2010). [CrossRef] [PubMed]

39. J. Liu, Y. Wang, X. Qu, X. Li, X. Ma, R. Han, Z. Hu, X. Chen, D. Sun, R. Zhang, D. Chen, D. Chen, X. Chen, J. Liang, F. Cao, and J. Tian, “In vivo quantitative bioluminescence tomography using heterogeneous and homogeneous mouse models,” Opt. Express 18(12), 13102–13113 (2010), http://www.opticsinfobase.org/abstract.cfm?URI=oe-18-12-13102. [CrossRef] [PubMed]

40. G. Yan, J. Tian, S. Zhu, Y. Dai, and C. Qin, “Fast cone-beam CT image reconstruction using GPU hardware,” J. XRay Sci. Technol. 16, 225–234 (2008).

41. X. Chen, X. Gao, D. Chen, X. Ma, X. Zhao, M. Shen, X. Li, X. Qu, J. Liang, J. Ripoll, and J. Tian, “3D reconstruction of light flux distribution on arbitrary surfaces from 2D multi-photographic images,” Opt. Express 18(19), 19876–19893 (2010), http://www.opticsinfobase.org/abstract.cfm?URI=oe-18-19-19876. [CrossRef] [PubMed]

42. K. Liu, J. Tian, Y. Lu, C. Qin, X. Yang, S. Zhu, and X. Zhang, “A fast bioluminescent source localization method based on generalized graph cuts with mouse model validations,” Opt. Express 18(4), 3732–3745 (2010), http://www.opticsinfobase.org/oe/viewmedia.cfm?URI=oe-18-4-3732&seq=0. [CrossRef] [PubMed]

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

Material	Muscle	Lungs	Heart	Liver	Kidneys	Stomach
$μ_{a}$ [mm⁻¹]	0.23	0.35	0.11	0.45	0.12	0.21
$μ_{s}^{'}$ [mm⁻¹]	1.00	2.30	1.10	2.00	1.20	1.70

τ	Recon. Method	Recon. Time(s)	Recon. Position(mm)	LE (mm)	Recon. Density (nW/mm³)	Recon. Power (nW)	RE of Power
1e-3	ℓ1-ℓs	87.37	(9.91,15.87,24.45)	0.47	0.0191	0.4402	18.54%
	GPSR	0.19	(8.46,16.90,24.55)	2.23	0.0056	0.0670	87.60%
	IVTCG	0.38	(9.91,15.87,24.45)	0.47	0.0201	0.4550	15.80%
1e-4	ℓ1-ℓs	357.22	(9.91,15.87,24.45)	0.47	0.0137	0.3299	38.95%
	GPSR	0.17	(8.46,16.90,24.55)	2.23	0.0040	0.0622	88.49%
	IVTCG	1.16	(9.91,15.87,24.45)	0.47	0.0145	0.3468	35.83%
1e-5	ℓ1-ℓs	1340.76	(9.91,15.87,24.45)	0.47	0.0132	0.3188	41.01%
	GPSR	0.17	(8.46,16.90,24.55)	2.23	0.0037	0.0593	89.03%
	IVTCG	4.21	(9.91,15.87,24.45)	0.47	0.0123	0.3021	44.10%

Case	Deviation in $μ_{a}$	Deviation in $μ_{s}^{'}$	Density (nW/mm³)	Power (nW)	RE of Power
1	−10%	−10%	0.0255	0.5660	4.74%
2	−10%	0	0.0240	0.5354	0.93%
3	−10%	+ 10%	0.0213	0.4790	11.36%
4	0	−10%	0.0235	0.5224	3.33%
5	0	0	0.0201	0.4550	15.80%
6	0	+ 10%	0.0167	0.3860	28.57%
7	+ 10%	−10%	0.0194	0.4393	18.71%
8	+ 10%	0	0.0154	0.3600	33.38%
9	+ 10%	+ 10%	0.0118	0.2846	47.34%

Source	Actual Position(mm)	Recon. Position(mm)	LE (mm)	Density (nW/mm³)	Power (nW)	RE of Power
Source-1	(10.2,15.5,24.5)	(9.91,15.87,24.45)	0.47	0.0260	0.5866	8.55%
Source-2	(10.2,16,28.5)	(9.33,16.46,28.63)	0.99	0.0232	0.4375	18.65%

Material	Muscle	Heart	Lungs	Liver	Kidneys
$μ_{a}$ [mm⁻¹]	0.0086	0.1382	0.4596	0.8291	0.1550
$μ_{s}^{'}$ [mm⁻¹]	1.2584	1.0769	2.2651	0.7356	2.5329

Sparse reconstruction for quantitative bioluminescence tomography based on the incomplete variables truncated conjugate gradient method

Abstract

1. Introduction

2. Bioluminescence tomography framework

2.1 Photon propagation model

2.2 BLT inverse model

3. Sparse reconstruction method

3.1 Sparse regularization

3.2 IVTCG method

4. Simulation studies in the mouse atlas

4.1 Quantitative reconstruction in a single source

4.2 Reconstruction robustness to measurement noise

4.3 Reconstruction robustness to optical parameters perturbation

4.4 Double source case evaluation

5. In vivo experiment validation

6. Conclusion and discussion

Acknowledgments

References and links

Cited By

Figures (7)

Tables (5)

Equations (20)

Optics Express