Computational and Mathematical Methods in Medicine

Volume 2012 (2012), Article ID 494808, 10 pages

doi:10.1155/2012/494808

Research Article

Sparse Reconstruction for Bioluminescence Tomography Based on the Semigreedy Method

Wei Guo,¹ Kebin Jia,¹ Qian Zhang,² Xueyan Liu,³ Jinchao Feng,¹ Chenghu Qin,² Xibo Ma,² Xin Yang,² and Jie Tian²

¹The College of Electronic Information & Control Engineering, Beijing University of Technology, Beijing 100124, China

²Intelligent Medical Research Center, Institute of Automation, Chinese Academy of Sciences, Beijing 100190, China

³Sino-Dutch Biomedical and Information Engineering School of Northeastern University, Liaoning, Shenyang 110004, China

Received 31 May 2012; Accepted 26 June 2012

Academic Editor: Huafeng Liu

Copyright © 2012 Wei Guo et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

Bioluminescence tomography (BLT) is a molecular imaging modality which can three-dimensionally resolve the molecular processes in small animals in vivo. The ill-posedness nature of BLT problem makes its reconstruction bears nonunique solution and is sensitive to noise. In this paper, we proposed a sparse BLT reconstruction algorithm based on semigreedy method. To reduce the ill-posedness and computational cost, the optimal permissible source region was automatically chosen by using an iterative search tree. The proposed method obtained fast and stable source reconstruction from the whole body and imposed constraint without using a regularization penalty term. Numerical simulations on a mouse atlas, and in vivo mouse experiments were conducted to validate the effectiveness and potential of the method.

1. Introduction

Due to its ability of monitoring physiological and pathological activities at the molecular level, small animal optical molecular imaging has become an important method for biomedical research. Bioluminescence imaging, as one of optical molecular imaging modalities, has attracted attention for its advantages in sensitivity, specificity, and cost effectiveness in cancer research and drug development [1–3]. Compared with planar bioluminescent imaging, BLT can three-dimensionally reconstruct the bioluminescent probes in small animals [4].

The generic BLT model is ill-posed. However, it has been theoretically proven that the solution uniqueness can be established under practical constraints using a priori knowledge [5]. In most existing reconstructions, multi spectral measurement [6–10], permissible source region (PSR) [9–12], and sparse reconstruction [13–15] are three common strategies to reduce ill posedness of BLT. Although multi spectral techniques improve reconstruction qualities to a certain degree by increasing the measurable data, they in turn impose some limitations in practical applications such as increased signal acquisition time and a high computational cost [16]. Besides, the PSR strategy can significantly improve the location accuracy of reconstructed source and reduce the computational cost by limiting the reconstruction region into a small area. However, in practical applications, both the size and position of the permissible region have significant impact on imaging results [5]. Additionally, since the bioluminescent source distribution is usually sparse in practical applications and only insufficient boundary measurements are available, the compressed sensing can bring benefits in spatial resolution and algorithm stability to BLT reconstruction. Recently, many sparse reconstruction methods have emerged in BLT [17–19]. The majority of them reformulate the BLT inverse problem into minimizing an objective functions that integrate a sparse regularization term with a quadratic error term and solve it via computationally tractable convex optimization methods, such as linear programming and gradient methods, However, the efficacy of the sparse regularization methods strongly depends on the choice of regularization parameter in practical applications [19].

The existing works have demonstrated that PSR can improve reconstruction qualities by reducing the number of unknown variables. Although the meaningful results can be obtained by using the PSR that is based on the bioluminescent signals and a priori knowledge available from a specific biomedical application [4, 11, 12], in most cases, it is rather difficultly to manually select such a small and appropriate region. Recently, some optimal permissible source region methods have emerged in BLT. Feng et al. presented a reconstruction algorithm for a spectrally resolved BLT based on an adaptive rough estimate of an optimal permissible source region and multilevel finite element method approach (FEM), where Tikhonov regularization was used to solve the constrained BLT inverse problem [9]. Naser and Patterson proposed a two-step reconstruction algorithm of bioluminescence, in which the permissible regions were shrunk by using an iterative minimization solution based on the L1 norm [10]. These works demonstrated the feasibility and potential of the optimal PSR techniques with numerical simulation. However, both of two previous reconstructions also needed the regularization methods to reconstruct the sources in the allowed region, which made the reconstruction results also depend on the choice of regularization parameter. Furthermore, they were demonstrated with only regular phantoms simulations and presented no in vivo experiment validation.

In this work, a novel BLT reconstruction algorithm based on the semi greedy method was proposed. The optimal PSR problem was cast into a search for the correct support of source distribution among a number of dynamically evolving candidate subsets, and the optimal PSR was chosen automatically by using an iterative search tree. Therefore, the columns of the system matrix were treated as the nodes for building up the search tree where each path from the root to a leaf node denoted a candidate. The search tree was initialized with some unspecific nodes. At each iteration, new nodes were appended to the most promising path, which were selected to minimize the cost function based on the residue. The permissible source region was expended by adding nodes with high a probability to contribute to the source. Among the system matrix, the columns that were corresponding to the nodes contained in the most promising path were selected to obtain the source distribution. By automatically choosing an optimal PSR, the method reduced the ill posedness of the problem and imposed constraint without using a regularization penalty term.

This paper is organized as follows. In Section 2, the forward photon propagation model, the inverse problem for BLT with FEM, and the proposed algorithm are introduced. In Section 3, the numerical simulations in a mouse atlas demonstrate the performance of the proposed method. In Section 4, an in vivo mouse experiment is conducted to further evaluate its reliability. Finally, we discuss the results and conclude this paper.

2. Method

Light propagation in biological media is essentially the transport of radiant energy. The radiative transfer equation (RTE) can rigorously describe light transport in turbid media [20]. Compared with the other approximations of RTE such as simple spherical harmonics, spherical harmonics and discrete ordinates, the following steady-state diffusion approximate equation (DA) in (1) is the most popular one as a result of its moderate computational efficiency and explicit physical meaning [4, 5, 11, 15, 21]. Assuming that the bioluminescence imaging experiment is performed in a totally dark environment and no photon travels into through the boundary , the equation is subject to the Robin boundary condition in (2) as follows:

where is the domain of the problem, donates the source energy distribution, represents the photo fluence rate, is the absorption coefficient, is the reduced scattering coefficient, indicates the optical diffusion coefficient, donates the boundary of the problem, and represents the mismatch coefficient between and its surrounding medium. The measured quantity on the boundary is given by the outgoing radiation as follows:

FEM is a powerful tool for solving the DA equation [4, 6–10]. By using FEM to discretely approach the solving domain and making a series of transformations and rearrangements, the linear relationship links the source distribution inside the heterogeneous medium, and the photon fluence rate on the surface is established as follows:

where is the source distribution of the interior nodes, is the measurable photon flux photon on the boundary nodes, is the positive definite matrix, and is the source weight matrix. The nonmeasurable entries in and corresponding rows in can be removed. Then a new linear relationship can be obtained as follows [4, 22]:

For BLT, the domains of the bioluminescent sources are usually very small and sparse compared with the entire reconstruction domain. That means that there are only fewer nonzero components in . Therefore, the system matrix can be seen as the dictionary, and has a less-term representation over the dictionary. As a consequence, the aim of the proposed method is imposing the constraint on the source space by choosing only the part that contributes to the source distribution. In the language of sparse approximation, greed pursuit algorithms are useful methods for solving this problem [23, 24]. For instance, Orthogonal Matching Pursuit (OMP) is to pick columns in a greedy fashion [25]. According to the introduction of OMP, the reconstruction starts with an empty index set. At each iteration, we choose the single column of that is most strongly correlated with the remaining part of . Then, we subtract its contribution to and iterate on the residual. The reconstruction is stopped after a number of iterations. Unfortunately, from experimental results shown in Figure 2 and Table 2, we found that OMP, as the single-path algorithm, could not achieve the desired expectations for reconstructing the bioluminescent source. In the experiments, when the computation of a single path selects a wrong column, the correct one is still in the set of candidate representations. Therefore, incorporation of a multipath search strategy is motivated to improve reconstruction. In this section, the semigreedy method was used to search for the correct support of among a number of dynamically evolving candidate subsets.

A general best first (GBF) is a search algorithm which constructs a tree by expanding the most promising node chosen according to a specified rule. Search algorithm , as one of the most studied versions of GBF, can find path in combinatorial search. It selects an optimal path by minimizing an additive evaluation function , where is the cost of the currently evaluated path from start node to , and is a heuristic estimate of the cost of the path remaining between and some goal nodes [26–29]. In our problem, the search tree was iteratively built up by nodes which represent the dictionary atoms. Each path from the root to a leaf node denoted a subset of dictionary atoms which was a candidate support for . Let us define the notation. and denote the atoms contained in path and the vector of corresponding coefficients obtained after orthogonal projection of the residue onto the set of selected columns [25]. Similarly, and represent the selected atom at the th node on path and corresponding coefficient.

The search tree starts with less unspecified nodes. A simple way is selecting the nodes that have the highest absolute inner product with . In order to find the fewest possible nodes, the search must constantly make an evaluation to decide which available paths should be expended next. Therefore, the evaluation function is defined as follows:

Beside the evaluation function, the auxiliary function is also needed to assess the cost brought by adding a preferred goal node to the path. Generally, according to the expectation on average equal contribution of unopened nodes, the auxiliary function can be built as follows:

where coefficient is defined as . is a ratio between the number of the nonzero entries and the zero entries in the solution . It is well known that the sparse solution has only less nonzero entries. Therefore, in most practical application, is very smaller than 1. Here, we selected . If the source distribution could be seen as the K-sparse signal, could be computed by . The cost function can be written as follows:

where is a constant. A lot of experiments for different reconstruction models including 2D and 3D experiments were performed to evaluate the impact of on the source reconstruction. We found that better results could be obtained when it varied in the interval [1.00, 1.25]. Therefore, in our experiments, its value was taken from the interval and set to be1.05.

In practice, if all of the children of the most promising partial path are added to the search tree at each iteration, the tree might have too many search paths. Therefore, the following pruning strategies are employed a guide on how the tree grows.

The first one is about extensions per path. At each step, it is not necessary to have all of the unopened atoms added to the current optimal path. We can expand the search tree only by the children which have the highest absolute inner-product with the residue to the selected path. This pruning strategy decreases the upper bound to on the number of paths. Practically, and are selected to be much smaller than , which can drastically decrease the paths involved in the search. Although the number of extensions per path is limited to , it is also necessary to control the size of path. That is because that adding new paths at each iteration continues increasing required memory. Therefore, to reduce memory requirements, we adopted the “beam search” strategy [30] and limit the maximum number of paths in the tree by the beam width . When this limit is exceeded, the paths with maximum cost are seen as the worst paths and are removed from the tree until paths remain. Here, and were set to be 4 and 200, respectively. Moreover, since order of nodes along a path is unimportant, amalgamating the equivalent path is also important to improve the search performance. For this purpose, we define a path equivalency notion; and are two paths with different length . If all atoms of can be found in and these composed the continuous subset of , we define the above two paths as being equivalent. Consequently, the insertion of into the tree is unnecessary.

After the growing of the search tree is finished, the linear relationship between the observation and the selected PSR can be established as follows:

Since the nodes contained in the optimal path can be much smaller than the number of all nodes , (9) is an overdetermined linear equation. Therefore, a limited memory variation of the Broyden Fletcher Goldfarb Shanno (LBFGS) [31] was used to directly solve (9).

3. Simulation Studies in the Mouse Atlas

In this subsection, heterogeneous simulations were presented to demonstrate the performance of the proposed method for mouse applications. The experimental data were acquired by a dual-modality BLT/micro-CT system developed in our lab [32, 33]. By using image processing and interactive segmentation technology, heterogeneous model including heart, lungs, liver, bone, spleen, and muscle was built. The optical coefficients for each organ are listed in Table 1 [34]. Here, the torso section with a height of 25 mm was selected as the reconstruction region.

Table 1: Optical parameters for each organ in the mouse atlas [34].

Muscle
Heart
Lungs
Liver
Spleen
Bone

(mm⁻¹)
0.032
0.022
0.071
0.128
0.075
0.002

(mm⁻¹)
0.586
1.129
2.305
0.646
2.178
0.935

Table 2: Quantitative comparisons of the reconstruction results.

Method
Recon. location center (mm)
Location error (mm)
Recon. time (s)
Maximum recon. value (nW/mm³)

OMP
(15.61, 32.92, 47.48)
2.9605
0.798
0.00412

FIST-L1
(18.26, 31.97, 47.28)
0.3995
16.94
0.01371

FIST-L1
(15.61, 32.92, 47.48)
2.9605
25.44
0.00414

The proposed method
(18.26, 31.97, 47.28)
0.3995
12.20
0.01889

3.1. Reconstruction in a Single Source Case

In the experiment, a spherical source with a 0.6 mm diameter was placed in the liver with the center at (18.24 mm, 31.58 mm, 47.29 mm) as shown in Figure 1(a). The source was modeled as isotropic point sources whose strength was set to be 2 nW/mm³. As for the forward problem, the FEM was used to obtain the synthetic measurements on the boundary. The atlas model was discretized into a tetrahedral-element mesh with 30892 nodes and 167841 elements. The generated simulated photon distribution on the boundary is presented in Figure 1(b). Then the forward solutions were projected onto a single coarse mesh consisting of 20068 elements and 3098 nodes, which was used for reconstructing the source.

(a)

(b)

Figure 1: Reconstruction model with a single source. (a) The torso of the mouse atlas model with one source in the liver. (b) The simulated photon distribution on the surface.

(a)

(b)

(c)

(d)

Figure 2: Comparison of reconstruction results. (a), (b), (c) and (d) are the reconstruction results with OMP, FIST-L1 (with regularization parameter set to be 4e-11), FIST-L1 (with regularization parameter set to be 4e-10), and the proposed method, respectively. The results are shown in the form of isosurfaces for 40% of the maximum value (left column) and slice images in  mm plane (right column). The small yellow sphere in the iso-surfaces view image and black circles in the slice images denote the real position of the bioluminescent source.

To better illustrate the performance of the proposed method, we compared the proposed method with OMP [25] and FIST-L1 [35, 36]. The former is a typical greedy pursuit method for sparse signal recovery. The latter, as a sparse regularization method, can be viewed as a standard approach to ill-posed linear inverse problems and has been used in fluorescence molecular tomography (FMT) and BLT. Here, the step size in FIST-L1 was computed by using the estimation algorithm introduced in [35]. Since the regularization parameter plays an important role in the regularization methods, we performed two experiments with different regularization parameters that were set to be 4e-11 and 4e-10, respectively. All of the reconstructions were carried out on a personal computer with 3.2 GHz Intel Core2 duo CPU and 2 GB RAM.

The qualities of the reconstruction were quantitatively assessed in terms of location error and the maximum reconstructed intensity. The location error was defined as Euclidean distance between and , where and were the real locations of the source center and the location of the node with the maximum reconstructed value, respectively. The visual effects of the reconstruction results are presented in the form of slice images and iso-surfaces, as shown in Figure 2. Additionally, the detailed information about parameters and the final quantitative reconstruction results are summarized in Table 2.

We found that the reconstructed positions by L1 regularization with an optimal regularization parameter and the proposed method were identical. Specifically, the reconstructed center was (18.26 mm, 31.97 mm, 47.28 mm) with a location error of 0.3995 mm from the actual source, whereas the location error by OMP was 2.9605 mm. The performance of OMP was inferior to the other two methods. L1 regularization performed well and obtained satisfactory source localizations and maximum reconstruction value. However, the selection of the regularization parameter had a great impact on the reconstruction results. As for the proposed method, it performed slightly better than L1 regularization with manually optimized regularization parameter in terms of maximum reconstruction value. Moreover, it was also an efficient reconstruction method.

The above experiments were performed without noise. In order to evaluate the sensitivity of the proposed method to various noise levels, six cases were carried out where the measurements were added to 5%, 10%, 15%, 20%, 25% and 30% Gaussian noise, respectively. We also made a comparison between the proposed method and FIST-L1 with the regularization parameter that was set to be 4e-11. The reconstruction results under different noise levels are compiled in Table 3 and Figure 3 showing that the proposed method was robust to measurement noise.

Table 3: Quantitative results in a single source case with different noise levels.

Method
Noise level
Recon. location

center (mm)
Location error (mm)
Recon. Time (s)
Maximum recon. value (nW/mm³)

FIST-L1
5%
(18.26, 31.97, 47.28)
0.3995
24.43
0.0122

10%
(18.26, 31.97, 47.28)
0.3995
24.99
0.0085

15%
(18.26, 31.97, 47.28)
0.3995
24.77
0.0050

20%
(15.61, 32.92, 47.48)
2.9605
24.82
0.0032

25%
(15.61, 32.92, 47.48)
2.9605
24.94
0.0025

30%
(15.61, 32.92, 47.48)
2.9605
24.60
0.0023

The proposed method
5%
(18.26, 31.97, 47.28)
0.3995
13.41
0.0192

10%
(18.26, 31.97, 47.28)
0.3995
14.15
0.0234

15%
(18.26, 31.97, 47.28)
0.3995
16.37
0.0204

20%
(18.26, 31.97, 47.28)
0.3995
18.09
0.0223

25%
(18.26, 31.97, 47.28)
0.3995
18.95
0.0116

30%
(18.26, 31.97, 47.28)
0.3995
46.61
0.0054

(a)

(b)

(c)

(d)

(e)

(f)

Figure 3: The proposed method for different noise levels of 0%, 5%, 10%, 15%, 20%, 25%, and 30%.

3.2. Double Source Case and MultiSource Case

Dual source setting was also considered in order to evaluate the proposed method. Two sources have the same size as one used in the single source case. Their strength and position were set to be 2 nW/mm³ and (18.24 mm, 31.58 mm, 47.29 mm), 1 nW/mm³ and (18.74 mm, 39.15 mm, 47.06 mm), respectively. The reconstruction results are shown in Figure 4 and Table 4.

Table 4: Reconstruction results in double source case.

Source number
Actual position (mm)
Recon. location center (mm)
Location error (mm)
Maximum recon. value (nW/mm³)

1
(18.24, 31.58, 47.29)
(18.26, 31.97, 47.28)
0.3995
0.0204

2
(18.74, 39.15, 47.06)
(18.75, 39.20, 46.86)
0.2064
0.0208

Figure 4: Reconstruction results in double source case. The results are shown in the form of iso-surface for 40% of the maximum value (left column) and slice image in  mm plane (right column). The small yellow sphere in the iso-surface view image and black circles in the slice image denote the real position of the bioluminescent source.

Multiple sources setting simulation experiment was also presented to further demonstrate the ability of the proposed method. Based on the setting in double source case, the third source with the same size and shape was added. Its strength and position were set to be 1.5 nW/mm³ and (23.60 mm, 37.94 mm, 47.45 mm). The final reconstruction results are presented in Figure 5 and Table 5. The result of two group experiments indicated that the sources can be accurately distinguished by using the proposed method.

Table 5: Reconstruction results in multiple source case.

Source number
Actual position (mm)
Recon. location center (mm)
Location error (mm)
Maximum recon. value (nW/mm³)

1
(18.24, 31.58, 47.29)
(18.26, 31.97, 47.28)
0.3995
0.0195

2
(18.74, 39.15, 47.06)
(18.75, 39.20, 46.86)
0.2064
0.0226

3
(23.60, 37.94, 47.45)
(23.60, 37.97, 47.65)
0.2023
0.0300

Figure 5: Reconstruction results in multisource case. The results are shown in the form of iso-surface for 40% of the maximum value (left column) and slice image in  mm plane (right column). The small yellow sphere in the iso-surface view image and black circles in the slice image denote the real position of the bioluminescent source.

4. In Vivo Experiment Validation

Besides the numerical simulations with mouse atlas, an in vivo experiment was carried out on a mouse to further test the proposed method. The experiment was also performed with a dual-modality BLT/micro-CT system developed in our lab [32, 33, 37]. A nude hairless mouse (Nu/Nu, Laboratory Animal Center, Peking University, China) was used in this experiment. To simulate a known bioluminescent source, a home-made cylindrical light source about 3 mm long and 1.5 mm in diameter was implanted into the abdomen of the mouse in this experiment. The source was made of a catheter filled with luminescent liquid and emitted a red luminescent light that had a similar emission spectrum with a firefly luciferase-based source.

Before the beginning of the experiment, the CCD (VersArray, Princeton Instruments, Trenton, NJ, USA) was cooled to 110° by using liquid nitrogen. The mouse was anesthetized and placed in a mouse holder. The mouse holder was set to rotate to 0°, 90°, 180° and 270°. At each of four positions, the mouse was photographed by the CCD camera. After the optical data were acquired, the mouse was scanned by using the micro-CT to obtain the anatomical maps which could provide structural information for the source reconstruction. Then the CT data were segmented into five regions represent muscle, lungs, heart, liver, and kidneys, respectively, as shown in Figure 6(). The heterogeneous model including five tissues was discretized into the mesh containing 11917 tetrahedral elements and 2557 nodes. The optical parameters for different tissues were calculated based on the literature as listed in Table 6 [34, 38]. The optical data was registered with the volumetric mesh, and measured data were mapped onto the surface of the mesh. The result of mapped photon distribution is shown in Figure 6().

Table 6: Optical parameters for each organ in the heterogeneous model [34, 38].

Muscle
Heart
Lungs
Liver
Kidneys

(mm⁻¹)
0.008
0.138
0.456
0.829
0.150

(mm⁻¹)
1.258
1.076
2.265
0.735
2.507

Figure 6: In vivo heterogeneous model. (a) The torso of the model. (b) The mapped photon distribution on the mouse surface.

It took about 8 seconds to complete the reconstruction using the proposed method. The final results are presented in Figure 7, where the reconstruction source center is (37.17 mm, 38.82 mm, and 20.92 mm) with a deviation about 2 mm to the actual center. As can be seen in the reconstruction results, the proposed methods could obtain satisfactory bioluminescent source localizations.

Figure 7: The reconstruction results. The results are shown in the form of isosurface for 40% of maximum value (left column) and slice images (right column). The transverse view of the results and the comparisons with the corresponding CT slices. The cross of the red lines denotes the actual source center.

5. Discussion and Conclusion

In this paper, we have proposed a new method based on the semigreedy for bioluminescence tomography. The reconstruction results of the simulations on a mouse atlas demonstrate that the proposed reconstruction method is able to accurately and stably localize bioluminescent source from whole body, even with noisy measurements. The in vivo experiment further shows its performance.

The PSR strategy can significantly improve reconstruction qualities. However, in most cases, empirically selecting such small and appropriate region is unconvenient, even available. In this study, The optimal PSR problem is cast into a search for the correct support of source distribution among a number of dynamically evolving candidate subsets. In view of the characteristics in BLT sparse distribution, only the columns that contribute to the source reconstruction are chosen automatically by using semi-greedy method. The constraint imposed on the source space reduces the ill posedness of the problem and computational cost.

It is noted that in vivo experiment is not as accurate as simulations. Some reasons can be explained for this phenomenon. First of all, the error was generated, when the energy distribution was mapped from 2D images to a 3D mouse surface. Secondly, there were only five main segmented tissues used to build a heterogeneous model while others simply were regarded as the muscle, which also led to errors. Finally, the accuracy of the photon propagation model was very important for source reconstruction. The diffusion approximation was used due to its moderate computational efficiency and explicit physical meaning. However, it has some limitations in certain regions, such as void or more absorptive regions. Therefore, the error brought on by the DA model is inevitable. As discussed above, our future work will focus on studying more accurate forward models to describe photon propagation in biological tissues and improving the experimental procedures and imaging system to further promote the reconstruction quality.

Acknowledgments

This paper is supported by the National Basic Research Program of China (973 Program) under Grants no. 2011CB707700, the Knowledge Innovation Project of the Chinese Academy of Sciences under Grant no. KGCX2-YW-907, the National Natural Science Foundation of China under Grant Nos. 30970780, 81000624, 81027002, 81071205, 30970778, 81071129, and 81101095, the Beijing Natural Science Foundation No. 4111004, the Science and Technology Key Project of Beijing Municipal Education Commission under Grant No. KZ200910005005, the Doctoral Fund of the Ministry of Education of China under Grant No. 20091103110005, the Fellowship for Young International Scientists of the Chinese Academy of Sciences under Grant No. 2010Y2GA03, and the Chinese Academy of Sciences Visiting Professorship for Senior International Scientists under Grant No. 2010T2G36.

References

V. Ntziachristos, J. Ripoll, L. V. Wang, and R. Weissleder, “Looking and listening to light: the evolution of whole-body photonic imaging,” Nature Biotechnology, vol. 23, no. 3, pp. 313–320, 2005.

J. K. Willmann, N. van Bruggen, L. M. Dinkelborg, and S. S. Gambhir, “Molecular imaging in drug development,” Nature Reviews Drug Discovery, vol. 7, no. 7, pp. 591–607, 2008.

A. M. Loening, A. M. Wu, and S. S. Gambhir, “Red-shifted Renilla reniformis luciferase variants for imaging in living subjects,” Nature Methods, vol. 4, no. 8, pp. 641–643, 2007.

G. Wang, W. Cong, K. Durairaj et al., “In vivo mouse studies with bioluminescence tomography,” Optics Express, vol. 14, no. 17, pp. 7801–7809, 2006.

G. Wang, Y. Li, and M. Jiang, “Uniqueness theorems in bioluminescence tomography,” Medical Physics, vol. 31, no. 8, pp. 2289–2299, 2004.

H. Dehghani, S. C. Davis, and B. W. Pogue, “Spectrally resolved bioluminescence tomography using the reciprocity approach,” Medical Physics, vol. 35, no. 11, pp. 4863–4871, 2008.

C. Kuo, O. Coquoz, T. L. Troy, H. Xu, and B. W. Rice, “Three-dimensional reconstruction of in vivo bioluminescent sources based on multispectral imaging,” Journal of Biomedical Optics, vol. 12, no. 2, Article ID 024007, 2007.

H. Dehghani, S. C. Davis, S. Jiang, B. W. Pogue, K. D. Paulsen, and M. S. Patterson, “Spectrally resolved bioluminescence optical tomography,” Optics Letters, vol. 31, no. 3, pp. 365–367, 2006.

J. Feng, K. Jia, G. Yan et al., “An optimal permissible source region strategy for multispectral bioluminescence tomography,” Optics Express, vol. 16, no. 20, pp. 15640–15654, 2008.

M. A. Naser and M. S. Patterson, “Algorithms for bioluminescence tomography incorporating anatomical information and reconstruction of tissue optical properties,” Biomedical Optics Express, vol. 1, no. 2, pp. 512–526, 2010.

C. Qin, X. Yang, J. Feng et al., “Adaptive improved element free Galerkin method for quasi-or multi-spectral bioluminescence tomography,” Optics Express, vol. 17, no. 24, pp. 21925–21934, 2009.

Y. Lv, J. Tian, W. Cong et al., “A multilevel adaptive finite element algorithm for bioluminescence tomography,” Optics Express, vol. 14, no. 18, pp. 8211–8223, 2006.

Y. Lu, H. B. Machado, A. Douraghy, D. Stout, H. Herschman, and A. F. Chatziioannou, “Experimental bioluminescence tomography with fully parallel radiative-transfer-based reconstruction framework,” Optics Express, vol. 17, no. 19, pp. 16681–16695, 2009.

H. Gao and H. Zhao, “Multilevel bioluminescence tomography based on radiative transfer equation part 1: 11 regularization,” Optics Express, vol. 18, no. 3, pp. 1854–1871, 2010.

X. He, J. Liang, X. Wang et al., “Sparse reconstruction for quantitative bioluminescence tomography based on the incomplete variables truncated conjugate gradient method,” Optics Express, vol. 18, no. 24, pp. 24825–24841, 2010.

Y. Lu, H. B. Machado, A. Douraghy, D. Stout, H. Herschman, and A. F. Chatziioannou, “Experimental bioluminescence tomography with fully parallel radiative-transfer-based reconstruction framework,” Optics Express, vol. 17, no. 19, pp. 16681–16695, 2009.

A. Cong, W. Cong, Y. Lu, P. Santago, A. Chatziioannou, and G. Wang, “Differential evolution approach for regularized bioluminescence tomography,” IEEE Transactions on Biomedical Engineering, vol. 57, no. 9, pp. 2229–2238, 2010.

M. Freiberger, C. Clason, and H. Scharfetter, “Total variation regularization for nonlinear fluorescence tomography with an augmented Lagrangian splitting approach,” Applied Optics, vol. 49, no. 19, pp. 3741–3747, 2010.

J. Feng, C. Qin, K. Jia et al., “An adaptive regularization parameter choice strategy for multispectral bioluminescence tomography,” Medical Physics, vol. 38, no. 11, pp. 5933–5944, 2011.

A. D. Klose, V. Ntziachristos, and A. H. Hielscher, “The inverse source problem based on the radiative transfer equation in optical molecular imaging,” Journal of Computational Physics, vol. 202, no. 1, pp. 323–345, 2005.

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,” Medical Physics, vol. 22, no. 11 I, pp. 1779–1792, 1995.

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,” Physics in Medicine and Biology, vol. 53, no. 14, pp. 3921–3942, 2008.

J. A. Tropp, A. C. Gilbert, and M. J. Strauss, “Algorithms for simultaneous sparse approximation. Part I: greedy pursuit,” Signal Processing, vol. 86, no. 3, pp. 572–588, 2006.

J. A. Tropp, “Algorithms for simultaneous sparse approximation. Part II: convex relaxation,” Signal Processing, vol. 86, no. 3, pp. 589–602, 2006.

G. Davis, S. Mallat, and M. Avellaneda, “Adaptive greedy approximations,” Constructive Approximation, vol. 13, no. 1, pp. 57–98, 1997.

S. Koenig, M. Likhachev, Y. Liu, and D. Furcy, “Incremental heuristic search in AI,” Artificial Intelligence Magazine, vol. 25, no. 2, pp. 99–112, 2004.

N. B. Karahanoglua and H. Erdogan, “Compressed sensing signal recovery via A* Orthogonal Matching Pursuit,” in Proceedings of the IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP '11), pp. 3732–3735, 2011.

S. Koenig, M. Likhachev, and D. Furcy, “Lifelong Planning A*,” Artificial Intelligence, vol. 155, no. 1-2, pp. 93–146, 2004.

R. Dechter and J. Pearl, “Generalized best-first search strategies and the optimality of A*,” Journal of the ACM, vol. 32, no. 3, pp. 505–536, 1985.

R. Zhou and A. Hansen, “Beam-stack search: integrating backtracking with beam search,” in Proceedings of the The International Conference on Automated Planning & Scheduling (ICAPS '05), pp. 90–98, 2005.

D. C. Liu and J. Nocedal, “On the limited memory BFGS method for large scale optimization,” Mathematical Programming B, vol. 45, no. 3, pp. 503–528, 1989.

J. Feng, K. Jia, C. Qin et al., “Three-dimensional bioluminescence tomography based on bayesian approach,” Optics Express, vol. 17, no. 19, pp. 16834–16848, 2009.

X. Ma, J. Tian, C. Qin et al., “Early detection of liver cancer based on bioluminescence tomography,” Applied Optics, vol. 50, no. 10, pp. 1389–1395, 2011.

Y. Lv, J. Tian, W. Cong et al., “Spectrally resolved bioluminescence tomography with adaptive finite element analysis: methodology and simulation,” Physics in Medicine and Biology, vol. 52, no. 15, pp. 4497–4512, 2007.

D. Han, J. Tian, S. Zhu et al., “A fast reconstruction algorithm for fluorescence molecular tomography with sparsity regularization,” Optics Express, vol. 18, no. 8, pp. 8630–8646, 2010.

D. Han, J. Tian, C. Qin, K. Liu, and X. Ma, “A fast reconstruction method for fluorescence molecular tomography based on improved iterated shrinkage,” in Proceedings of the Medical Imaging : Biomedical Applications in Molecular, Structural, and Functional Imaging (SPIE '011), Lake Buena Vista, Fla, USA, February 2011.

B. Zhang, X. Yang, C. Qin et al., “A trust region method in adaptive finite element framework for bioluminescence tomography,” Optics Express, vol. 18, no. 7, pp. 6477–6491, 2010.

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,” Physics in Medicine and Biology, vol. 50, no. 17, pp. 4225–4241, 2005.