## Abstract

Inverse source reconstruction is the most challenging aspect of bioluminescence tomography
(BLT) because of its ill-posedness. Although many efforts have been devoted to this problem, so
far, there is no generally accepted method. Due to the ill-posedness property of the BLT
inverse problem, the regularization method plays an important role in the inverse
reconstruction. In this paper, six reconstruction algorithms based on
*l*_{p} regularization are surveyed. The effects of the permissible
source region, measurement noise, optical properties, tissue specificity and source locations
on the performance of the reconstruction algorithms are investigated using a series of single
source experiments. In order to further inspect the performance of the reconstruction
algorithms, we present the double sources and the *in vivo* mouse experiments to
study their resolution ability and potential for a practical heterogeneous mouse experiment. It
is hoped to provide useful guidance on algorithm development and application in the related
fields.

© 2012 OSA

## 1. Introduction

In the past decade, bioluminescence imaging (BLI) has been one of the hot topics in optical imaging and has been successfully applied in tumorigenesis studies, cancer diagnosis, metastasis detection, drug development, gene therapy, etc. [1–5]. However, the intrinsic planar imaging mode of BLI is just a qualitative imaging tool, which provides no absolute quantitative information and detects limited depths of a few millimeters inside tissues [4,5]. The bioluminescence tomography (BLT), the three-dimensional (3D) tomographic counterpart of BLI, has received much attention because it can provide localization and quantitative information of the bioluminescent source inside a living small animal [6–9].

In BLT, the tumor cells of the living mouse model are usually transfected with a luc-gene. Then, the oxidation of the substrate and luciferin results in a red-shifted light emission with a wavelength in the range from 500 nm to about 760 nm [3]. Based on the surface bioluminescent images and the heterogeneous tissue structure information captured by, for example, the micro-CT imaging system, BLT can reconstruct the location and distribution of the internal embedded bioluminescent light source and acquire quantitative volumetric reconstruction information [8].

The inverse source reconstruction is crucial for recovering three dimensional information of
the light source. Many efforts have been devoted to studies on efficient inverse reconstruction
[10–28].
Generally, some prior knowledge, such as the permissible source region and multispectral
information, are used to reduce the ill-posedness of inverse reconstruction [10–14,16,17]. Meanwhile, the
development of a fast and robust algorithm is another challenging problem. The common schemes
include the regularization method, iterative method, stochastic method and so on [15,18–26]. Since the compressive sensing (CS) theory was introduced
in 2006, which is a signal processing technique for efficiently acquiring and reconstructing
images with sparse representation from far fewer measurements or signals below the Nyquist
sampling theorem demands, it has been utilized by various areas of applied mathematics, computer
science, electrical engineering and medicine [29,30]. As a new imaging technique, the sparse bioluminescent
source and insufficient measuring signal are the characteristics of BLT, which can benefit by
using the CS technique. Along with the development of the compressive sensing (CS) theory, the
*l*_{p} regularization method has become the mainstream strategy to
obtain the optimal solution of the BLT inverse problem [17,22,24–26].

Generally speaking, regularization is used to construct the approximation of an ill-posed
problem by a family of neighboring well-posed problems [31]. So far, the regularization methods have been used for solving various types of
inverse problems in mathematics, statistics, machine learning, signal and image processing, heat
conducting, and so on. The most commonly used form of regularization is the Tikhonov-type
regularization. This method achieves a compromise between the minimization of the residual norm
and the penalty term [32]. The
*l*_{2}-norm is the integrant penalty term for the Tikhonov-type
regularization. The *l*_{2}-norm penalty term imposes a small weight on
small residuals, but a large weight on others. As a result, the
*l*_{2}-norm approximation has many modest residuals and relatively few
larger ones, which generates smooth solutions and brings multipseudo sources to surround the
true source in BLT reconstruction [22]. Considering the
source sparsity characteristics and the insufficiency of the measured data, the
*l*_{1}-norm penalty term has been adopted in the regularization method
and was successfully applied to BLT [24–26,33]. Unlike
*l*_{2}-norm regularization, *l*_{1}-norm
regularization tends to produce a sparse solution with a large number of components equal to
zero. However, in many practical applications, the solutions of the
*l*_{1} regularizer are often less sparse than those of the
*l*_{0} regularizer. The *l*_{0} regularizer aims
to find the number of non-zero elements in the solving domain, but it is an impossible mission
in practical applications [34]. To find more sparse
solutions than the *l*_{1} regularizer, *l*_{p} (0
< p <1) regularization methods have been developed and employed to solve the minimization
problems [34–36] and the BLT reconstruction problem [28].

Although the *l*_{p} (0 < p ≤ 1, p = 2) regularization
methods have been applied to the BLT reconstruction problem, so far, there is no generally
accepted method which can be suitable for all of the reconstruction cases. This paper intends to
fill the gap in the existing studies to systematically benchmark the performance of the
*l*_{p}-regularization-based BLT reconstruction algorithms. Our group
has been conducting research in the BLT reconstruction method for years. In this paper, we
evaluated six reconstruction algorithms, including the Truncated Total Least Square method
(TTLS) [27], the Incomplete Variables Truncated Conjugate
Gradient method (IVTCG) [25], the Truncated Newton
Interior-Point method (TNIPM) [24], the Primal-Dual
Interior-Point method (PDIP) [26], and the Weighted
Iterative Shrinkage/Thresholding algorithm (WISTA) [28]
that were developed by our group, and the Tikhonov regularization that we used for the BLT
reconstruction in its standard way [31]. In order to
investigate the response of these algorithms to the permissible source region, measurement
noise, optical properties, and tissue specificity, we conducted a series of single source
numerical phantom experiments. Then, the double source numerical phantom experiment and the
*in vivo* mouse experiment were carried out to further test their
performance.

This paper is organized as follows. The next section describes the diffusion approximation
model of BLT and the inverse reconstruction formula. The essence of the six regularization
algorithms are then introduced briefly. In Section 3, the response of these algorithms to
several influence factors are investigated and discussed based on the numerical phantom and
*in vivo* mouse experiments to demonstrate the performance of the different
algorithms. Finally, the discussion and conclusions are addressed in Section 4.

## 2. Methodology

#### 2.1. Inverse reconstruction formula of BLT

By assuming light as photons, the propagation of light in biological tissues is an extremely complex process that involves a series of optical behaviors, including the absorption, scattering, internal reflection and transmission. The radiative transfer equation (RTE) is often used to analytically model the process of light propagation. Although the accuracy of the RTE model is indubitable, it is generally difficult to directly solve it in complex biological tissues. Furthermore, the complicated implementation in the numerical settings and the large numbers of the resulting equations makes its efficiency low for BLT [37]. In infrared optical imaging, most of the biological tissues exhibit the characteristics of high scattering and low absorption. As the first order approximation model of RTE, the diffusion approximation (DA) is widely used to depict the process [8]:

where $S\left(r\right)$ is the source distribution, *r* is the position vector, $\Phi \left(r\right)$ is the photon flux density, ${\mu}_{a}(\text{r})$ is the absorption coefficient, $D(\text{r})=1/(3({\mu}_{a}(\text{r})+(1-g){\mu}_{s}(\text{r})))$ is the diffusion coefficient, ${\mu}_{s}(\text{r})$ is the scattering coefficient, *g* is the anisotropy parameter, and $\Omega \subset {R}^{3}$ is the domain of interest. Equation (2) is the Robin boundary condition, and $v\left(r\right)$ is the unit outer normal to $\partial \Omega $ at r, ${A}_{n}(r)=(1+R(r))/(1-R(r))$ is the boundary mismatch factor and $R(r)$ depends on the refractive index *n* of the
surrounding medium that can be approximated by $R\approx -1.4399{n}^{-2}+0.7099{n}^{-1}+0.6681+0.0636n$ [38]. In the optical
imaging experiment, the outgoing flux density on $\partial \Omega $ is measured by the highly sensitive CCD camera and can be
expressed as follows [8]:

To solve Eqs. (1) and (2), the finite element method (FEM) was adopted and
the equations were discretized by the finite element basis functions. A linear relationship
between the measurable boundary data *b* and the unknown source distribution *S* is then established and the BLT forward problem can be obtained
as [39]:

where the weight matrix A establishes the relationship between the measurements and the unknowns.

Because the measurement is insufficient compared with the whole solving domain, the inverse
reconstruction is an ill-posed problem. As a result, it is difficult and impossible to solve
Eq. (4) directly. Usually, the regularization
method is employed to convert Eq. (4) into the
following optimization problem and *S* is replaced by *x* to make the notation more general:

where ${\Vert x\Vert}_{\text{p}}={\left({\displaystyle \sum _{i\text{=}1}^{n}{\left|{x}_{i}\right|}^{\text{p}}}\right)}^{1\text{/p}}$ and p denotes the parameter used in the
*l*_{p} regularization and satisfies 0 < p ≤ 1 or p = 2.
Selecting an appropriate optimization method, the localization and distribution of the
bioluminescent probes can be obtained from Eq.
(5).

#### 2.2. Regularization methods

### 2.2.1. Tikhonov regularization

Tikhonov regularization method is one of the most commonly used methods to solve the ill-posed inverse problems. When p = 2, Eq. (5) is converted into the following minimization problem:

where $\text{\lambda >0}$ is a properly chosen regularization parameter and controls the
weight given to the minimization of the regularization penalty term relative to the
minimization of the residual norm term. In the BLT application, the right-hand side *b*
is always contaminated by various types of errors, such as the
measurement errors, approximation errors, and rounding errors. Assuming the right-hand side *b*
satisfies the discrete Picard condition, an explicit solution of (6)
is given by

By using the singular value decomposition (SVD) explicitly, the system matrix A can be expressed as $A=U\Sigma {V}^{T}$ and the Tikhonov regularized solution is given by

where $U=({u}_{1},\mathrm{....},{u}_{n})$ and $V=({v}_{1},\mathrm{....},{v}_{n})$ are matrices with orthonormal columns, ${U}^{T}U={V}^{T}V={I}_{n}$, and where $\Sigma =diag({\sigma}_{1},\mathrm{...},{\sigma}_{n})$ has non-negative diagonal elements appearing in non-increasing order such that ${\sigma}_{1}\ge \cdots \ge {\sigma}_{n}\ge 0$ [40].

### 2.2.2. Truncated Total Least Squares Method

The Tikhonov regularization method only considers the measurement errors on the right-hand
side *b*, but the system errors also exist in the computed coefficient
matrix of the system equation. In fact, whether the system errors come from the coefficient
matrix or the measurement noise, they are inevitable in a real BLT experiment, so He
*et al.* proposed the truncated total least squares (TTLS) method combined
with a practical parameter-choice scheme termed as improved generalized cross-validation
(IGCV) to cope with the BLT reconstruction problem including both the measurement noise and
the system errors [27].

The TTLS method takes its basis in an SVD of the augmented matrix $(A\text{}b)=\overline{U}\overline{\Sigma}{\overline{V}}^{T}={\displaystyle \sum _{i=1}^{n+1}{\overline{u}}_{i}}{\overline{\sigma}}_{i}{\overline{v}}_{i}^{T}$, $\overline{U}\in {\text{R}}^{m\times (n+1)}$, $\overline{V}\in {\text{R}}^{(n+1)\times (n+1)}$, where ${\overline{U}}^{T}\overline{U}={\overline{V}}^{T}\overline{V}={I}_{n+1}$, and $\Sigma =diag({\overline{\sigma}}_{1},\mathrm{...},{\overline{\sigma}}_{n+1})$ have singular values ${\overline{\sigma}}_{1}\ge \cdots \ge {\overline{\sigma}}_{n+1}\ge 0$. The matrix $\overline{V}$ is partitioned such that $\overline{V}=\left(\begin{array}{cc}{\overline{V}}_{11}& {\overline{V}}_{12}\\ {\overline{V}}_{21}& {\overline{V}}_{22}\end{array}\right)$, ${\overline{V}}_{11}\in {\text{R}}^{n\times k}$, and ${\overline{V}}_{22}\in {\text{R}}^{1\times (n+1-k)}$, where *k* is the truncation parameter. Then, the TTLS solution is given by
${x}_{\text{TTLS}}=-{\overline{V}}_{12}{\overline{V}}_{22}^{\u2020}=-\frac{{\overline{V}}_{12}{\overline{V}}_{22}^{T}}{{\Vert {\overline{V}}_{22}\Vert}_{2}^{2}}$. Based on the modified GCV (MGCV) criterion for TTLS
$G=\frac{{\Vert Ax-b\Vert}^{2}}{{(m-{\text{enp}}_{k})}^{2}}$, where ${\text{enp}}_{k}$is the sum of the filter factors of TTLS, the truncation
parameter *k* is obtained by the IGCV scheme.

### 2.2.3. Incomplete Variables Truncated Conjugate Gradient (IVTCG) Method

In view of the insufficient measurement, the source sparse characteristic and the
oversmoothing of *l*_{2}-norm regularization in BLT, linear system
Eq. (4) is covert to the following
*l*_{1}-norm regularization problem based on the CS theory:

where ${\Vert x\Vert}_{1}={\displaystyle {\sum}_{i}{x}_{i}}$ is the *l*_{1} norm of
*x*. Introducing vectors *u* and *v*, and the
substitution $x=u-v$, Eq. (9) can be
reformulated as the following convex quadratic program with nonnegative constrained
conditions:

where $u\ge 0$, $v\ge 0$, $z=[u{\text{}v]}^{T}$, $c=\lambda {1}_{2N}+[-q{\text{}q]}^{T}$, ${1}_{2N}={\left[1,1,\cdots ,1\right]}^{T}$, $q={A}^{T}x$, and $B=[A{\text{}-A]}^{T}[A\text{}-A]$. Then, the incomplete variables truncated conjugate gradient
(IVTCG) algorithm was presented by He *et al.* to solve Eq. (10) [25]. First, we define the working set ${\Gamma}_{k}$ by the violation of the Karush-Kuhn-Tucker (KKT) optimal
conditions at the *k*th iteration,

where ${z}^{k}$ is the current iteration point. Then, we define the other two working sets ${I}_{k}$ and ${J}_{k}$ as follows:

where $\delta >0$, ${N}_{s}>0$and ${N}_{\mathrm{max}}>{N}_{s}$are three constants, ${\widehat{J}}_{k}={\Gamma}_{k}\backslash {I}_{k}$ and ${\widehat{I}}_{k}=\left\{i\in \left\{1,\cdots ,2N\right\}|\left[{\left({z}^{k}\right)}_{i}>0,{\left({z}^{k}\right)}_{i}/{\left(\nabla F({z}^{k})\right)}_{i}>\delta \right]\right\}$. The variables with the set will be updated ${I}_{k}$and ${J}_{k}$ at the current iteration.

Finally, we need to decide the search direction. For the search direction ${d}_{{I}_{k}}^{k}$of ${I}_{k}$, it can be obtained by solving the sub-problem

For ${J}_{k}$, the search direction ${d}_{{J}_{k}}^{k}\text{=}-{w}_{{J}_{k}}$. The other entries of ${d}^{k}$ are set to 0. After determining the search direction, the next iteration point is derived by the backtracking method.

In the IVTCG method, the parameters ${N}_{s}$and ${N}_{\mathrm{max}}$are set as ${N}_{s}=\lfloor \raisebox{1ex}{$M$}\!\left/ \!\raisebox{-1ex}{$4$}\right.\rfloor $and ${N}_{\mathrm{max}}={N}_{s}+\lfloor \raisebox{1ex}{${N}_{s}$}\!\left/ \!\raisebox{-1ex}{$8$}\right.\rfloor $. The time complexity of the IVTCG algorithm is $O(M{m}^{2})+O(MN)$.

### 2.2.4. Truncated Newton interior-point method (TNIPM)

Using the truncated Newton interior-point method (TNIPM), the object function (9) is transformed to a convex quadratic problem with inequality constraints:

Then, a logarithmic barrier function for the constraints $\left|{x}_{i}\right|\le {u}_{i}$ is added to the minimization function (14), and Eq. (14) is converted to the following problem:

where parameter $t\in (0,\infty )$.

In order to compute the search direction, a preconditioned conjugate gradient method is adopted to compute the Newton system

Then, we can obtain a feasible point by using the initial value and the search direction. The implementation of TNIPM in BLT reconstruction is presented in [24].

### 2.2.5. Primal-dual interior-point method (PDIP)

As an efficient method for solving the optimization problems, interior-point methods have
been applied in many areas in the past twenty years. In [26], the primal-dual interior-point method was adopted to solve the BLT inverse
problem. Equation (4) is transformed to the
following *l*_{1}-norm minimization problem:

and then converted to the standard primal and dual forms in linear programming:

By introducing a logarithmic barrier term, (P) is converted to the following logarithmic barrier problems:

where *θ* is the barrier parameter. Then, we can obtain the
Karush-Kuhn-Tucker conditions for ${\text{P}}_{\theta}$:

where $e={\left(1,\mathrm{...},1\right)}^{T}$, the matrix *X* and *S* are the $n\times n$ diagonal matrix whose diagonal entries are precisely the
components of x and *s* respectively.

Hence, a new feasible point and the Newton direction $\left(\Delta x,\Delta y,\Delta s\right)$ at the current value $\left({x}^{k},{y}^{k},{s}^{k}\right)$can be obtained by solving the following equation system:

The general framework of PDIPM is presented in [26].

### 2.2.6. Weighted iterative shrinkage/thresholding algorithm

In some cases, the *l*_{1}-norm regularization methods are often less
sparse than *l*_{p}-norm (0 < p < 1) regularization methods. In
the paper [28], we proposed a weighted iterative
shrinkage/thresholding algorithm (WISTA) to solve the *l*_{p}-norm
minimization problem (5) (0 < p < 1). Because we cannot solve the
*l*_{p}-norm minimization problem directly, we reformed it by the
weighted *l*_{1}-norm:

In order to obtain an ideal solution of (23), we used the iterative form as follows:

where the weight vector *ω* is calculated from the last iterate ${x}^{i}$, soft(*z*, *q*) is the
soft-threshold function with a threshold *q*, and $\text{soft(}x\text{,}q\text{)=sign(}x\text{)max(}\left|x\right|\text{-}q\text{,0)}$. For the implementation of WISTA, we set two stopping
conditions. The first is that the quantity of descending in the objective function is less
than the stopping threshold. The second is that the algorithm reached the maximum number of
iterations.

The final algorithm of WISTA is formulated in [28].

## 3. Experiments and results

In this section, a series of numerical experiments and an *in vivo* experiment
were conducted to benchmark the performance of the six investigated reconstruction algorithms.
We used two indicators to evaluate the reconstruction results: the reconstruction time (Time) of
these algorithms and the distance (Loc_Err) between the actual source position and the
reconstruction position. We defined the reconstruction position as the point with the maximum
density of the reconstructed source. The regularization parameter of the Tikhonov method was
chosen by the classical L-curve method and the parameter-choice scheme of the TTLS method was
the improved generalized cross validation (IGCV). The regularization parameters of TNIPM and
IVTCG used in reconstruction were manually optimized and they were set as 1e-6 in this paper.
All experiments were performed in MATLAB on a personal computer with 2.66 GHz Intel Core 2 Quad
CPU and 4 GB RAM.

#### 3.1. Numerical Experiments

### 3.1.1. Reconstruction in a single source

In this part, a cylindrical heterogeneous phantom was employed to demonstrate the features of the investigated algorithms. The numerical phantom consisted of muscle, lungs, heart, bone and liver, with the dimensions of 30 mm in height and 20 mm in diameter, as shown in Fig. 1 [39]. The optical parameters of the five tissues were specified as listed in Table 1 . A solid sphere source with a 1 mm radius was centered at (3, 5, 0) inside the right lung. Meanwhile, 3874 nodes and 17763 tetrahedral elements were used in the inverse reconstruction. In order to compare the performance of the six algorithms among various conditions, we designed five experiments for the single source reconstruction.

### 1. Reconstruction using different permissible source regions

As *a priori* knowledge, the permissible source region (PSR) is commonly
used in the inverse reconstruction of BLT. Different PSR strategies may lead to different
reconstruction results. In Fig. 2
, we can see the surface flux distribution of the heterogeneous phantom in the coronal,
sagittal and translucency views respectively and six red lines indicate the PSR. We set four
different PSRs to verify its influence on the inverse reconstruction of the six algorithms.
The PSRs were defined as follows: PS1 = {(x, y, z)| the left lung}; PS2 = {(x, y, z)| 0<x,
y<7, −5< z <5}; PS3 = {(x, y, z)| −5<x, y<8, −5< z
<5}; and PS4 = {(x, y, z)| the whole phantom}. The volume ratios of the four PSR were
3.26%, 5.20%, 17.93% and 100% respectively. The reconstruction results are listed in Figs. 3
6. Since TTLS could only resolve the overdetermined equation, it did not provide results
for PS4. In Fig. 6(b), because the reconstruction time
of PDIP in PS4 was too long so as to make the time curves almost indistinguishable in the
Cartesian coordinate system, we used the log plot to depict the reconstruction time of
various methods.

We can draw the conclusion that the smaller the PSR, the better the reconstruction results
obtained. From Figs. 3–5, we found that the reconstruction results of *l*_{p} (0 < p
≤ 1) regularization methods (TNIPM, IVTCG, PDIP, WISTA) were better than those of
*l*_{2} regularization methods (TTLS, Tikhonov). From Fig. 3, the reconstruction results of
*l*_{2} regularization methods were scattered around the true source.
On the contrary, from Fig. 4
and Fig. 5, the reconstruction results of
*l*_{p} (0 < p ≤ 1) regularization methods were more
accurate and the reconstruction sources were smaller than those of
*l*_{2}, except for the results of PDIP.

From Fig. 6, the smallest location error was obtained by PDIP for PS1 and the second-smallest location error was obtained by IVTCG for PS2. As the permissible source region became larger, the location error of PDIP increased. The worst result of IVTCG was obtained for PS1 and it may be due to the fact that IVTCG couldn’t develop its capacity of sparse reconstruction in a small PSR. The reconstruction results of TNIPM and WISTA remained stable with different PSRs. The reconstruction times of WISTA and IVTCG were smaller than the remaining algorithms and PDIP was the most time consuming algorithm. On average, the reconstruction locations of IVTCG, TNIPM and WISTA outperformed the other algorithms.

### 2. Reconstruction at different measurement noise levels

In this experiment, random noise at different levels was added to the synthetic data to
evaluate the stability and robustness of the six algorithms using the permissible source
region of PS1. The reconstruction results under different noise levels (0%-50%) of the
Gaussian noise added to the synthetic data are shown in Fig.
7
and Table 2
[41]. The noise was generated by a
MATLAB function *randn* and each algorithm was tested with a different noise
sample. Tikhonov, TNIPM, IVTCG and WISTA were stable at all noise levels. The location error
of TTLS increased a little bit at the noise level of 10%, and then it was kept stable until
the noise level of 40%. The reconstruction results of PDIP behaved similarly. At the noise
level of 50%, both the TTLS and PDIP methods became unstable, so we ran them several times
and just showed the best results in Fig. 8(a)
and Fig. 8(e) respectively. In Fig. 7(b), TNIPM has a larger computation cost and IVTCG
methods are fast but produce larger localization errors.

### 3. Reconstruction using different optical parameters

In order to investigate the influence of optical property on inverse reconstruction, four
experiments with different deviations of absorption and scattering coefficients were
conducted by adding 10% perturbation to the absorption and scattering coefficients of
different organs. The corresponding results are shown in Table 3
and Fig. 9
. Therein, the symbols were defined as: P1 = {10% deviation of
${\mu}_{a}$ and ${\mu}_{s}$}; P2 = {-10% deviation of ${\mu}_{a}$ and ${\mu}_{s}$}; P3 = {10% deviation of ${\mu}_{a}$ and −10% deviation of ${\mu}_{s}$}; and P4 = {-10% deviation of ${\mu}_{a}$and 10% deviation of ${\mu}_{s}$}. From these results, we found that the results of P4 were
better than those of the others and the best reconstruction location was obtained by IVTCG in
P4. The reason should be that the optical parameters exhibit the characteristics of high
scattering and low absorption in P4, which conforms to the application scope of the diffusion
approximation model. From Fig. 9, IVTCG achieved the
smallest location error and WISTA had the least amount of computation time. The optical
property changes had little impact on the inverse reconstruction of WISTA and TNIPM. The
*l*_{p} (0 < p ≤ 1) regularization methods (TNIPM, IVTCG,
PDIP, WISTA) still obtained better reconstruction sources than the
*l*_{2} regularization methods (TTLS, Tikhonov).

### 4. Reconstruction in tissue specificity

In this part, the influence of tissue specificity on the inverse reconstruction was
observed in four experiments. In the first case of N1, the absorption, scattering and
anisotropy coefficients of all organs were set the same as that of muscle to simulate a 3D
homogeneous model. In the cases of N2-N4, we simulated three heterogeneous models (see Fig. 10
). In the case of N2, the absorption, scattering and anisotropy coefficients of the
lungs, heart and liver were set the same as those of muscle, thus the model was divided into
two parts: muscle and bone. In the case of N3, the absorption, scattering and anisotropy
coefficients of the lungs and heart were set the same, that is, the model was divided into
four parts. The case N4 was the PS1 in the PSR experiment and there were five tissues in the
model. The reconstruction results are presented in Table
4
and Fig. 11
. From these results, we found that the reconstruction results of the case N4 were
better than the other three cases for *l*_{p} (0 < p ≤ 1)
regularization methods (TNIPM, PDIP, WISTA), but this could not be observed for the
*l*_{2} regularization algorithms (TTLS, Tikhonov) and IVTCG. It may
be that there were more optical parameters in tissues in N4, so it was more approachable in
reality than the other three cases. As shown in Fig.
11(a), the PDIP’s performance was the best in the four cases and the
reconstruction results of WISTA ranked second. In Fig.
11(b), WISTA still used the minimal amount of time, followed by IVTCG and
Tikhonov.

### 5. Reconstruction at different source locations

In this experiment, we carried out three experiments to observe the influence of the source
locations on the inverse reconstruction. In case 1, the source location was (4,6,0) and the
distance between the source center and the surface of phantom was 2.789mm. In case 2, that is
PS1, the source location was (3,5,0) and the distance between the source center and the
surface of phantom was 4.169mm. In case 3, the source location was (2,4,0) and the distance
between the source center and the surface of phantom was 5.527mm. The reconstruction results
are listed in Table 5
and the error bar can be seen in Fig. 12
. From the results, we observed that the smaller the depth was, the better the
reconstruction results obtained. The standard deviation of WISTA was smallest indicating that
it’s robust to source locations. Except for PDIP, the *l*_{p}
(0 < p ≤ 1) regularization methods (TNIPM, IVTCG, WISTA) obtained smaller standard
deviations than *l*_{2} regularization methods (TTLS, Tikhonov).

### 3.1.2. Reconstruction of double sources

In this part, three experiments were carried out to verify the performance of these
algorithms on the double source resolution. In experiment 1, two sources were located in the
left lung with their centers at (3, 5, 1.25) and (3, 5, −1.25), and with a distance of
2.5mm. In experiment 2, two sources were located in the left lung with their centers at (3, 5,
1. 5) and (3, 5, −1.5), and with a distance of 3mm. In experiment 3, two sources were
located in the left lung with their centers at (3, 5, 2) and (3, 5, −2), and with a
distance of 4mm. The right lung was specified as the PSR in the three experiments. The
reconstruction results are presented in Figs. 13
–15
and Table 6
. From these results, we found that TTLS only reconstructs one source in case 1
and case 3, as shown in Fig. 13(a) and Fig. 13(c). In Fig.
13(f), the Tikhonov’s reconstruction of the two sources was almost overlapping.
On the contrary, *l*_{p} (0 < p ≤ 1) regularization methods
(TNIPM, IVTCG, PDIP, WISTA) obtained much better results. Comparing Fig. 14
with Fig. 15, we found that the reconstruction
sources of WISTA were smaller than that of the three *l*_{1}
regularization methods (TNIPM, IVTCG, PDIP), which confirm that WISTA was sparser than the
*l*_{1} regularization methods.

#### 3.2. In vivo mouse experiment

The *in vivo* experiment was presented here to further verify the performance
of these algorithms for the application of living animal studies. In the mouse experiment, an
implanted source filled with luminescent liquid was placed into the abdomen of the living
mouse. The source was 2.7 mm in diameter and 5.1 mm in length. The experimental details and the
optical parameters of each organ were presented in [26].
In this experiment, the whole mouse was specified as the PSR, which made the weight matrix A in
Eq. (5) undetermined. The reconstructed results
are presented in Fig. 16
and Table 7
. The TTLS did not provide results because it could not resolve the undetermined
equation. From Fig. 16(c), we found that the
reconstruction location of TNIPM was on the bottom left corner of the mouse model and we
concluded that its reconstruction failed. From Table
7, we can see that WISTA obtained the best reconstruction results and the run time of
PDIP was the longest.

## 4. Discussion and Conclusions

In this paper, we provided a comprehensive assessment of six regularization algorithms
developed by our group. In order to investigate the performance of these algorithms in the BLT
inverse reconstruction, we carried out five single source experiments, a double source
experiment and an *in vivo* experiment. Because the inverse reconstruction
equation is not well posed, the reconstruction results were influenced by various factors, such
as the permissible source region, measurement noise, optical properties, tissue specificity and
source positions. In this paper, we observed all of the five influential factors which affected
the results of BLT reconstruction in the single source experiments. In the double source case,
the ability of the six algorithms in resolving double sources was investigated. Finally, the
feasibility and practicability of these algorithms were further validated by the *in
vivo* experiment.

The extensive experiments showed that, under a wide range of conditions, there isn’t
always a winner that achieves the best performance for all of the investigated cases. Various
results for the six reconstruction algorithms under different conditions were obtained and some
interesting conclusions could be addressed. Consistent with our anticipation, the reconstruction
results became better as the PSR decreased, especially for the *l*_{2}
regularization methods (TTLS, Tikhonov), which could be mainly attributed to the fact that the
weight matrix was overdetermined and the sparsity characteristics of the inverse problem
didn’t pan out. In the noise experiment, four algorithms (Tikhonov, TNIPM, IVTCG, WISTA)
were robust to noise. The PDIP became unstable at a noise level of 50% and TTLS became unstable
at a noise level of 40%. The optical property variation experiments suggested that we could get
better reconstruction results if the mathematical model used to describe the light propagation
in the biological tissue fits the optical characteristics of the tissues. The tissue specificity
experiments demonstrated the advantage of the heterogeneous model over the homogeneous model,
which has been verified by many studies. But for *l*_{2} regularization
methods (TTLS, Tikhonov), especially the Tikhonov method, the reconstruction results
didn’t improve significantly, so we should use the *l*_{p} (0 <
p ≤ 1) regularization algorithms in the reconstruction of the heterogeneous model. In the
last single source experiments, WISAT was the most robust algorithm at different source
positions than the other five algorithms.

In the reconstruction of double sources, TTLS had two cases and Tikhonov had one case where
they couldn’t separate the two sources, which contrasted sharply with the reconstruction
results of the *l*_{p} (0 < p ≤ 1) regularization algorithms
(TNIPM, IVTCG, PDIP, WISTA). It is due to the fact that *l*_{2}
regularization is easy in generating smooth solutions, which leads to multi-pseudo sources that
can’t distinguish the two sources.

In the *in vivo* mouse experiment, except for the TNIPM method, other
algorithms performed normally. As in the double source experiment, WISTA still obtained the best
performance because it could maintain its preferable sparsity characteristics. We saw that the
reconstruction results of the *in vivo* mouse experiment were unsatisfactory
compared to the simulation experiment. Actually, many factors influenced the reconstruction
results of the *in vivo* experiment, such as the signal acquisition, image
registration of the CT data and optical data, the quality of mesh dissection, the accuracy of
the optical propagation model, the inverse reconstruction algorithm and so on.

For most experiments, the *l*_{p} (0 < p ≤ 1) regularization
algorithms (TNIPM, IVTCG, PDIP, WISTA) achieved better reconstruction results than
*l*_{2} regularization (TTLS and Tikhonov). In the four
*l*_{p} (0 < p ≤ 1) regularization algorithms, IVTCG and WISTA
were generally more effective than the other two methods.

In order to evaluate the *l*_{p}-regularization-based BLT
reconstruction algorithms, we limited the review to six algorithms (TTLS, Tikhonov, TNIPM,
IVTCG, PDIP, WISTA) that were developed by our group. This is mainly because our group has been
engaged in the research of bioluminescence tomography for years, we have the source code of the
six algorithms and we can show their best performance and give them a fair comparison. Some
interesting conclusions were obtained in this paper which hope to provide useful information for
the researcher in selecting and designing algorithms for BLT reconstruction.

## Acknowledgments

This work was supported by the Program of the National Basic Research and Development Program of China (973) under Grant No. 2011CB707702, the National Natural Science Foundation of China under Grant Nos. 81090272, 81101083, 81101084, 81101100, 81000632, 30900334, the National Key Technology Support Program under Grant No. 2012BAI23B06, the Natural Science Basic Research Plan in Shaanxi Province of China under Grant No. 2012JQ4015, and the Fundamental Research Funds for the Central Universities.

## References and links

**1. **R. Weissleder and M. J. Pittet, “Imaging in the era of molecular
oncology,” Nature **452**(7187), 580–589
(2008). [CrossRef] [PubMed]

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

**3. **M. Rodriguez-Porcel, J. Wu, and S. Gambhir,
“Molecular imaging of stem cells,” in *StemBook* [Internet]
(Harvard Stem Cell Institute,Cambridge, MA, 2008), available from http://www.ncbi.nlm.nih.gov/books/NBK27079/

**4. **F. S. Azar and X. Intes, *Translational
Multimodality Optical Imaging* (Artech House, 2008), Chap. 17.

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

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

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

**8. **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). [CrossRef] [PubMed]

**9. **G. Wang, W. Cong, H. Shen, X. Qian, M. Henry, and Y. Wang, “Overview of bioluminescence tomography--a new molecular
imaging modality,” Front. Biosci. **13**(13), 1281–1293
(2008). [CrossRef] [PubMed]

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

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

**12. **Y. Lv, J. Tian, W. Cong, G. Wang, W. Yang, C. Qin, and M. Xu, “Spectrally resolved bioluminescence tomography with
adaptive finite element analysis: methodology and simulation,”
Phys. Med. Biol. **52**(15), 4497–4512
(2007). [CrossRef] [PubMed]

**13. **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). [CrossRef] [PubMed]

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

**15. **J. Feng, K. Jia, C. Qin, G. Yan, S. Zhu, X. Zhang, J. Liu, and J. Tian, “Three-dimensional bioluminescence tomography based on
Bayesian approach,” Opt. Express **17**(19), 16834–16848
(2009). [CrossRef] [PubMed]

**16. **C. Qin, J. Tian, X. Yang, J. Feng, K. Liu, J. Liu, G. Yan, S. Zhu, and M. Xu, “Adaptive improved element free Galerkin method for
quasi- or multi-spectral bioluminescence tomography,” Opt.
Express **17**(24), 21925–21934
(2009). [CrossRef] [PubMed]

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

**18. **H. Gao and H. Zhao, “Multilevel bioluminescence tomography based on
radiative transfer equation Part 1: l1 regularization,” Opt.
Express **18**(3), 1854–1871
(2010). [CrossRef] [PubMed]

**19. **H. Gao and H. Zhao, “Multilevel bioluminescence tomography based on
radiative transfer equation part 2: total variation and l1 data
fidelity,” Opt. Express **18**(3), 2894–2912
(2010). [CrossRef] [PubMed]

**20. **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). [CrossRef] [PubMed]

**21. **B. Zhang, X. Yang, C. Qin, D. Liu, S. Zhu, J. Feng, L. Sun, K. Liu, D. Han, X. Ma, X. Zhang, J. Zhong, X. Li, X. Yang, and J. Tian, “A trust region method in adaptive finite element
framework for bioluminescence tomography,” Opt.
Express **18**(7), 6477–6491
(2010). [CrossRef] [PubMed]

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

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

**24. **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] [PubMed]

**25. **X. He, J. Liang, X. Wang, J. Yu, X. Qu, X. Wang, Y. Hou, D. Chen, F. Liu, and J. Tian, “Sparse reconstruction for quantitative bioluminescence
tomography based on the incomplete variables truncated conjugate gradient
method,” Opt. Express **18**(24), 24825–24841
(2010). [CrossRef] [PubMed]

**26. **Q. Zhang, H. Zhao, D. Chen, X. Qu, X. Chen, X. He, W. Li, Z. Hu, J. Liu, J. Liang, and J. Tian, “Source sparsity based primal-dual interior-point method
for three-dimensional bioluminescence tomography,” Opt.
Commun. **284**(24), 5871–5876
(2011). [CrossRef]

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

**28. **Q. Zhang, X. Qu, D. Chen, X. Chen, J. Liang, and J. Tian, “Experimental three-dimensional bioluminescence
tomography reconstruction using the *l*_{p}
regularization,” Adv. Sci. Lett. **16**(1), 125–129
(2012). [CrossRef]

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

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

**31. **H. W. Engl, M. Hanke, and A. Neubauer,
*Regularization of Inverse Problems* (Springer, 2000).

**32. **A. N. Tikhonov, “Solution of incorrectly formulated problems and the
regularization method,” Soviet Math. Dokl. **4**, 1035–1038
(1963).

**33. **D. Donoho, “For most large underdetermined systems of linear
equations the minimal l1-norm near solution is also the sparest
solution,” Commun. Pure Appl. Math. **59**(6), 797–829
(2006). [CrossRef]

**34. **Z. Xu, H. Zhang, Y. Wang, X. Chang, and L. Yong, “L1/2 regularization,”
Sci. China Inform. Sci. **53**(6), 1159–1169
(2010). [CrossRef]

**35. **X. Chen, F. Xu, and Y. Ye, “Lower bound theory of nonzero entries in solutions of
l_{2}-l_{p} minimization,” SIAM J. Sci.
Comput. **32**(5), 2832–2852
(2010). [CrossRef]

**36. **M. Wei, W. Scott, J. James, H. McClellan, and G. Larson, “Estimation of the discrete spectrum of relaxations for
electromagnetic induction responses using l_{p}-regularized least squares for 0
≤ p ≤ 1,” IEEE Geosci. Remote Sens.
Lett. **8**(2), 233–237
(2011). [CrossRef]

**37. **M. Chu, K. Vishwanath, A. D. Klose, and H. Dehghani, “Light transport in biological tissue using
three-dimensional frequency-domain simplified spherical harmonics
equations,” Phys. Med. Biol. **54**(8), 2493–2509
(2009). [CrossRef] [PubMed]

**38. **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), 1779–1792
(1995). [CrossRef] [PubMed]

**39. **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). [CrossRef] [PubMed]

**40. **A. Ribés and F. Schmitt, “Linear inverse problems in
imaging,” IEEE Signal Process. Mag. **25**(4), 84–99
(2008). [CrossRef]

**41. **R. Han, J. Liang, X. Qu, Y. Hou, N. Ren, J. Mao, and J. Tian, “A source reconstruction algorithm based on adaptive
hp-FEM for bioluminescence tomography,” Opt. Express **17**(17), 14481–14494
(2009). [CrossRef] [PubMed]