Three-dimensional Bioluminescence Tomography based on Bayesian Approach

Bioluminescence tomography (BLT) poses a typical ill-posed inverse problem with a large number of unknowns and a relatively limited number of boundary measurements. It is indispensable to incorporate a priori information into the inverse problem formulation in order to obtain viable solutions. In the paper, Bayesian approach has been firstly suggested to incorporate multiple types of a priori information for BLT reconstruction. Meanwhile, a generalized adaptive Gaussian Markov random field (GAGMRF) prior model for unknown source density estimation is developed to further reduce the ill-posedness of BLT on the basis of finite element analysis. Then the distribution of bioluminescent source can be acquired by maximizing the log posterior probability with respect to a noise parameter and the unknown source density. Furthermore, the use of finite element method makes the algorithm appropriate for complex heterogeneous phantom. The algorithm was validated by numerical simulation of a 3-D micro-CT mouse atlas and physical phantom experiment. The reconstructed results suggest that we are able to achieve high computational efficiency and accurate localization of bioluminescent source. © 2009 Optical Society of America OCIS codes: (100.3190) Inverse problems; (110.6960) Tomography; (170.3010) Image reconstruction techniques; (170.3660) Light propagation in tissues. References and links 1. V. Ntziachristos, J. Ripoll, L. V. Wang, and R. Weissleder, “Looking and listensing to light: the evolution of whole-body photonic imaging,” Nat. Biotechnol 23, 313-320 (2005). 2. C. H. Contag and M. H. Bachmann, “Advances in bioluminescence imaging of gene expression,” Annu. Rev. Biomed. Eng. 4, 235-260 (2002). 3. B. W. Rice, M. D. Cable, and M. B. Nelson, “In vivo imaging of light-emitting probes,” J. Biomed. Opt. 6, 432-440 (2001). 4. Z. Paroo, R. A. Bollinger, D. A. Braascb, E. Ricber, and D. R. Corey, “Validating Bioluminescence Imaging as a High-Throughput, Quantitiative Modality for Assesing Tumor Burden,” Mol. Imaging 3, 117-124 (2004). 5. G. Wang, E. A. Hoffman, G. McLennan, L. V. Wang, M. Suter, and J. Meinel, “Development of the first bioluminescenct CT Scanner,” Radiology 229(P), 566 (2003). #111857 $15.00 USD Received 26 May 2009; revised 23 Aug 2009; accepted 1 Sep 2009; published 8 Sep 2009 (C) 2009 OSA 14 September 2009 / Vol. 17, No. 19 / OPTICS EXPRESS 16834 6. E. E. Graves, J. Ripoll, R. Weissleder, and V. Ntziachristos, “A submillimeter resolution fluorescence molecular imaging system for small animal imaging,” Med. Phys. 30 901-911 (2003). 7. T. F. Massoud and S. S. Gambhir, “Molecular imaging in living subjects: seeing fundamental biological processes in a new light,” Genes Dev. 17, 545-580 (2003). 8. M. Allard, D. Côté, L. Davidson, J. Dazai, and R. M. Henkelman, “Combined magnetic resonance and bioluminescence imaging of live mice,” J. Biomed. Opt. 12(3), 034018-1-11 (2007). 9. 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, 6765-6771 (2005), http://www.opticsinfobase.org/oe/abstract.cfm?URI=OPEX-13-18-6756. 10. W. Cong and G. Wang, “Iterative method for bioluminescence tomography based on the radiative transport equation,” Proc. SPIE 6318, 631826-1-7 (2006). 11. G. Wang, Y. Li, and M. Jiang, “Uniqueness theorems in bioluminescence tomography,” Med. Phys. 31, 22892299 (2004). 12. 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, 323-345 (2005). 13. A. P. Gibson, J. C. Hebden, and S. R. Arridge “Recent advances in diffuse optical imaging,” Phys. Med. Biol. 50 R1-R43 (2005). 14. X. Gu, Q. Zhang, L. Larcom, and H. Jiang, “Three-dimensional bioluminescence tomography with model-based reconstruction,” Opt. Express 12, 3996-4000 (2004), http://www.opticsinfobase.org/oe/abstract.cfm?URI=OPEX-12-17-3996. 15. 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, 4225-4241 (2005). 16. 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, 5421-5441 (2005). 17. N. V. Slavine, M. A. Lewis, E. Richer, and P. P. Antich, “Iterative reconstruction method for light emitting sources based on the diffusion equation,” Med. Phys. 33, 61-68 (2006). 18. 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, 7801-7809 (2006), http://www.opticsinfobase.org/abstract.cfm?URI=oe-14-17-7801. 19. 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, 8211-8223 (2006), http://www.opticsinfobase.org/abstract.cfm?URI=oe-14-18-8211. 20. G. Wang, H. Shen, W. Cong, S. Zhao, and G. Wei, “Temperature-modulated bioluminescence tomography,” Opt. Express 14, 7852-7871 (2006), http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-14-17-7852. 21. M. Jiang, T. Zhou, J. Cheng, W. Cong, and G. Wang, “Image reconstruction for bioluminescence tomography from partial measurement,” Opt. Express 15, 11095-11116 (2007), http://www.opticsinfobase.org/abstract.cfm?URI=oe-15-18-11095. 22. Y. Lv, J. Tian, W. Cong, G. Wang, W. Yang, C. Qin, and M. Xu, “Spectrally resolved bioluminescence tomography with adaptive finite element: methodology and simulation,” Phys. Med. Biol. 52 4497-4512 (2007). 23. Q. Zhang, L. Yin, Y. Tan, Z. Yuan, and H. Jiang, “Quantitative bioluminescence tomography guided by diffuse optical tomography,” Opt. Express 16, 1481-1486 (2008), http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-16-3-1481. 24. 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, 15640-15654 (2008), http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-16-20-15640. 25. J. Tian, J. Bai, X. Yan, S. Bao, Y. Li, W. Liang, and X. Yang, ”Multimodality molecular imaging,” IEEE Eng. Med. Biol. Mag. 27, 48-57 (2008). 26. J. C. Ye, K. J. Webb, and C. A. Bouman, “Optical diffusion tomography by iterative coordinate descent optimization in a Bayesian framework,” J. Opt. Soc. Am. A 16, 2400-2412 (1999). 27. J. C. Ye, C. A. Bouman, K. J. Webb, and R. P. Millane, “Nonlinear multigrid algorithms for bayesian optical diffusion tomography,” IEEE Trans. Image Process. 10, 909-922, (2001). 28. H. Li, J. Tian, F. Zhu, W. Cong, L. V. Wang, E. A. Hoffman, and G. Wang, “A Mouse Optical Simulation Environment (MOSE) to Investigate Bioluminescent Phenomena in the Living Mouse with Monte Carlo Method,” Acad. Radiol. 11, 1029-1038 (2004). 29. S. R. Arridge, M. Schweiger, M. Hiraoka, and D. T. Delpy, “A finite element approach for modeling photon transport in tissue,” Med. Phys. 20, 299-309 (1993). 30. S. R. Arridge, “Optical tomography in medical imaging,” Inverse Probl. 15, 41-93 (1999). 31. 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, 1779-1792 (1995). 32. C. A. Bouman and K. Sauer, “A generalized Gaussian imaging model for edge-preserving MAP estimation,” #111857 $15.00 USD Received 26 May 2009; revised 23 Aug 2009; accepted 1 Sep 2009; published 8 Sep 2009 (C) 2009 OSA 14 September 2009 / Vol. 17, No. 19 / OPTICS EXPRESS 16835 IEEE Trans. Image Process. 2, 296-310 (1993). 33. S. S. Saquib, K. M. Hanson, and G. S. Cunningham, “Model-based image reconstruction from time-resolved diffusion data,” Proc. SPIE 3034, 369-380 (1997). 34. A. Mohammad-Djafari, “Joint estimation of parameters and hyperparameters in a Bayesian approach of solving inverse problems,” In Proceedings of IEEE Internation Conference on Image Processing 2, 473-476 (1996). 35. P. E. Gill, W. Murray, and M. Wright, Practical optimization, (Academic Press, New York, 1981). 36. E. G. Birgin, J. M. Martinez, “A box-constrained optimization algorithm with negative curvature directions and spectral projected gradients,” Computing, Sup. 15, 49-60 (2001). 37. E. G. Birgin, J. M. Martinez, “Large-scale Active-Set Box-Constrained Optimization Method with Spectral Projected Gradients,” Comput. Optim. Appl. 23, 101-125, (2002). 38. N. Dehiolanis, T. Lasser, D. Hyde, A. Soubret, J. Ripoll, and V. Ntziachristos, “Free-space fluorescence molecular tomography utilizing 360◦ geometry projections,” Opt. Lett. 32, 382-384 (2007). 39. H. Meyer, A. Garofalakis, G. Zacharakis, S. Psycharakis, C. Mamalaki, D. Kioussis, E. N. Economou, V. Ntziachristos, and J. Ripoll, “Noncontact optical imaging in mice with full angular coverage and automatic surface extraction,” Appl. Opt. 46, 3617-3627 (2007). 40. T. Chen, “Digital Camera System Simulator and applications,” Ph. D. Thesis, Stanford University (2003). 41. D. Qin, H. Zhao, Y. Tanikawa, and F. Gao, “Experimental determination of optical properties in turbid medium by TCSPC technique,” Proc. SPIE 6434, 64342E (2007). 42. C. Qin, J. Tian, X. Yang, K. Liu, G. Yan, J. Feng, Y. Lv, and M. Xu, “Galerkin-based meshless methods for photon transport in the biological tissue,” Opt. Express 16, 20317-20333 (2008), http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-16-25-20317. 43. S. Zhu, J. Tian, G. Yan, C. Qin and J. Liu, “An experimental cone-beam micro-CT system for small animal imaging,” Proc. SPIE 7258, 72582S (2009). 44. G. Yan, J. Tian, S. Zhu, Y. Dai, and C. Qin, “Fast cone-beam CT image reconstruction using GPU hardware,” IJ. X-Ray Sci. Tech


Introduction
Molecular imaging, especially small-animal molecular imaging, is rapidly becoming an active research filed [1].Among all the imaging modalities that are used for molecular imaging, optical techniques stand out because of their ability to reveal molecular and cellular activities in a small animal in vivo [1,2,3,4], including bioluminescence tomography (BLT) [5] and fluorescence tomography (FMT) [6], and so on.Their mechanism is to identity light source from the collected emitted photons from multiple 3-D directions with respect to a living mouse marked by reporter luciferases, which has been applied to cancer research including studies of tumor burden, response to therapy, assessment of gene expression, and development of metastatic lesions [7].In this paper, we focus exclusively on bioluminescence tomography.
There has been a great effort lately devoted to transforming bioluminescence imaging from a 2-D, planar bioluminescent imaging technique into a truly 3-D tomographic imaging modality applicable to small animals, because planar bioluminescent imaging cannot generate depth information [3,8].Therefore, bioluminescence tomography (BLT) has become a research hot spot.Compared with other imaging modalities, the advantages of BLT are sensitivity, speed, non-invasiveness, low cost and background noise [9,10].The aim of BLT is to reconstruct an internal bioluminescent source distribution based on both the outgoing bioluminescent photons captured using a highly sensitive charge-coupled device (CCD) camera and a prescanned tomographic volume, such as a CT/Micro-CT or MRI volume, of the same mouse [11].
Bioluminescent photons propagation in biological tissue is governed by the radiative transfer equation (RTE) [12].However, to solve for the distribution of bioluminescent source in a large 3-D volume, as is required in BLT, the RTE is computationally expensive [13].Given the dominance of scattering over absorption in bioluminescent imaging, RTE can be replaced by diffusion equation and Robin boundary condition, which provide a quite accurate description of propagation of the photon transport.However, a major difficulty in determining the bioluminescent source distribution is imposed by multiple scattering of photons that propagate through heterogeneous biological tissues.In addition, the absence of external illumination accords a high sensitive signal, it complicates the tomographic problem.
Wang et al. [11] has theoretically proved the solution uniqueness for BLT under practical constrains and presented that the BLT problem can be reduced to an inverse source problem by incorporating sufficient a priori information.The commonly used a priori information is the permissible source region strategy and its importance in BLT reconstruction has been well recognized [9].Now, various reconstruction methods and real experiments have been developed and reported based on the strategy [10,18,19,21,22,24] and promising results have demonstrated the merits and potential of the strategy.Aside from the permissible source region strategy, the anatomic structure information of the small animal which can be imaged by Xray CT and MRI techniques can be considered as a priori information.The geometrical shape of major organ regions can be acquired by 3-D computer graphic techniques.Then, the organ region can be associated with tissue optical parameters which can be determined by diffusion optical tomography (DOT).Based on the permissible source region strategy, anatomic structure information and tissue optical parameters, the ill-posedness of BLT will be further reduced and the viable solutions can be solved [25].
Bayesian approach provides a nature framework to incorporate multiple types of a priori information.Recently, Bayesian approachs have been applied to nonlinear inverse problems such as diffusion optical tomography (DOT) [26,27].But surprisingly, Bayesian approach has not been applied to the BLT problem.In addition, Bayesian reconstruction methods for DOT can not be directly applied to BLT problem because of the differences between DOT problem and BLT problem.Furthermore, the known Bayesian methods used in DOT problem are not appropriate to complex phantoms such as real mouse [27].Therefore, there is a critical need to derive and construct a method based on Bayesian approach that can accurately reconstruct bioluminescent sources in heterogeneous and complex tissues.
In the work, various a priori information is incorporated into a Bayesian network to reduce the ill-posedness of BLT problem and an adaptive Gaussian Markov random field prior model for unknown source density estimation is developed on the basis of finite element theory.Under the framework of Bayesian approach, we develop a reconstruction algorithm to identify the bioluminescent source distribution in complex phantoms.The paper is organized as follows.The next section describes the details of the proposed reconstruction algorithm.In the third section, numerical experiments are presented to demonstrate the feasibility and merits of the algorithm with a 3-D micro-CT mouse atlas.With the model, the reconstructions were carried out from the data on the surface of the phantom, which were synthesized by molecular optical simulation environment (MOSE) developed by our group [28].Furthermore, real physical phantom experiment was performed to valuate the algorithm.Finally, Section 4 discusses the relevant issues and makes concluding remarks.

Forward model for BLT
In bioluminescence imaging experiment, it is necessary to depict the propagation of the photon transport accurately.Generally, the bioluminescent photon is subject to both scattering and absorption in biomedical tissues.However, in the wavelength range of bioluminescent light, the biomedical tissues are highly scattering, photon scattering predominates over photon absorption.Assuming the bioluminescent source density is stable when photons are collected, the photon propagation in biological tissues can be modeled by steady-state diffusion equation and Robin boundary condition [9,29,30]: where r ∈ R 3 represents the location vector, Ω and ∂ Ω are the domain and its boundary respectively; Φ(r) denotes the photon flux density [Watts/mm 2 ]; x(r) is the source energy density [Watts/mm 3 ]; μ a (r) is the absorption coefficient is the optical diffusion coefficient [mm], μ s (r) the scattering coefficient [mm −1 ] and g the anisotropy parameter; ν(r) is the unit outer normal on ∂ Ω.Given the mismatch between the refractive indices n for Ω and n for the external medium, A(r; n, n ) can be approximately represented: where n is close to 1.0 when the mouse is in air; R(r) can be approximated by: R(r) ≈ −1.4399n −2 + 0.7099n −1 + 0.6681 + 0.0636n [31].The measured quantity is the outgoing flux density Q(r) on ∂ Ω, that is:

Formulation of the Bayesian Problem
We use the vector x to denote the set of unknown source density.Φ meas k is assumed as the measured photon flux density at the kth detector position (k = 1, 2,...,M) and then organize the measurements as a single column vector y: y = [Φ meas

The data likelihood model
In a Bayesian framework, the data likelihood p(y|x) is required.The bioluminescence experiment is generally operated at low temperature, photon detection can be modeled using shot noise statistics, then the data likelihood can be given by [26]: Where α is the parameter related to the noise variance, Λ is the diagonal covariance matrix, ω 2 Λ = ω T Λω, and the vector value function f (x) represents the exact value of the outgoing flux for the assumed value of the source density x.For BLT, we can assume that the measurements are statistically independent with the variance of each measurement equal to its mean [27], so Λ is diagonal.In our simulations, we approximately assume Λ as: To compute the likelihood of the data, we need to solve the diffusion equation.The finite element method (FEM) has been used for the purpose and the detailed formulation of f (x) can be determined.The FEM method used in the paper is proposed in [19].In finite element analysis, the given domain Ω can be discretized into N T tetrahedron elements and N v vertex nodes.
Taking the permissible source region information PS into account, there are N p tetrahedron elements in the permissible source region which represent the possible unknown source distribution.Then the column vector x is reduced to the set of N p elements' source density distribution.

Source reconstruction
Bayes requires that we must assign a prior to unknown variable x.When the generalized adaptive Gaussian Markov random field (GAGMRF) prior model is used [32,33], Where σ is a normalization hyperparameter and 1 ≤ p ≤ 2 with p = 2 corresponding to the Gaussian case.Given the ith and jth tetrahedron elements, their four local vertexes are E ie and E je (e = 1, 2, 3, 4), respectively.If one of vertex of E ie is identified to any vertex of jth element, we assume that jth element is the adjacent element.b i j denotes the weighting assigned to be inversely proportional to the pair {i, j}, for each i, the weights b i j sum to 1. N is adaptively determined in the algorithm, which consists of all pairs of adjacent element.z(p) is partition function.If α is unknown, take the Eqs.( 9) and ( 11), the source reconstruction problem can be stated as following optimization problem: In the optimization process, α is adaptively estimated, which can be solved by viewing the Eq. ( 12) as a cost function of α, and setting the derivative with respect to α equal to zero.Then we can obtain: By substituting Eq. ( 13) into Eq.( 12), the optimization problem is converted into where x is an estimate of the unknown source density x.After neglecting constant terms, the Eq. ( 14) can simplified as : In practical calculation procedure, for maximizing l(x) by maximizing with respect to α and x using the following equations: Eq. ( 16) is a straight-forward computation, and Eq. ( 17) is a computationally expensive optimization problem.In order to calculate Eq. ( 17), the spectral projected gradient-based largescale optimization algorithm is employed [26,27,34,35,36,37].In the optimization procedure, each tetrahedron element of the phantom is updated sequentially.After every tetrahedron element has been updated, the procedure is repeated, starting from the first tetrahedron again.
We define a single update of every tetrahedron as a scan.There are a number of scans in the optimization procedure until the convergence criterion is satisfied [26].
As far as convergence criterion is concerned, we use the discrepancy between the measured and computational boundary nodal flux data, maximum number of scan, or the discrepancy between the current scan and last scan log posterior probability to evaluate if the procedure should be terminated, that is,

Experiments and results
A series of experiments were designed to evaluate our proposed algorithm, including numerical and physical phantom experiments.Reconstructions were carried out on a personal computer with 2.8 GHz Pentium 4 processor and 1.75 GB RAM.

Phantom and Synthetic data
In order to optimally establish and evaluate the proposed algorithm, a 3-D mouse atlas of micro-CT was employed to provide 3-D anatomical information.For this, we first prepared a 20 g male BALB/c mouse because of its broad applications in biological and medical research and convenient tail vein injection in albino mouse.All animal studies were performed according to procedures approved by the Animal Care and Use Committee.After fasting for 24 hours, the mouse tail was immersed in warm water for 30-60 seconds to increase blood flow to the tail and dilate the vessels prior to injection of the contrast agent.Then the Fenestra LC at a dose of 15 mL/kg was slowly injected over a period of 30-60 seconds via the lateral tail vein.Finally the mouse was anesthetized with intraperitoneal injection of a 13% aqueous urethane (0.15 mL/20 g body weight).
The anesthetized mouse was scanned 40 minutes after injection of the contrast agent using a micro-CT system.During scanning, the X-ray source was operated in a continuous mode with a 1 mm aluminum filter, and the typical tube voltage was 50 kVp with tube current 1.4 mA [43].The source to detector distance (SDD) and source to object distance (SOD) were set to 541.36 mm and 448.50 mm, respectively.500 projection views with image size of 2240 × 2344 pixels were obtained by 360 • full scan.Finally, the FDK filtered backprojection method using GPU hardware acceleration was performed for 3-D reconstruction [44].4.0 23.0 16.0 20.0 6.0 g 0.9 0.94 0.85 0.9 0.9 After micro-CT scanning and 3-D reconstruction, the real mouse was manually segmented into different tissue organs, shown in Fig. 1.In the following bioluminescence experiments, only the thorax micro-CT images were used, including lung, bone, heart, liver and muscle as shown in Fig. 2(a).Optical parameters from the literature [15,19] and the references therein were assigned to each of the five components, as listed in Table 1.For finite-element-based BLT reconstruction, the mouse phantom was discretized into a tetrahedral-element mesh, illustrated in Fig. 2(b).For inverse problems, numerical simulations of reconstruction methods usually make use of the synthetic data from the forward problem.Therefore, the inverse crime is likely to occur when closely correlated computational ingredients are used in the forward solver and the inverse scheme [19].In order to avoid the inverse crime problem, the synthetic data was produced by Molecular Optical Simulation Environment (MOSE) which is developed based on Monte Carlo method [28].In the single source simulation, the solid spherical source with 1mm radius and source density of 0.238nano − Watts/mm 3 was centered at (22.8, 28.6, 12.5) inside the lung.The bioluminescent source was sampled by 1.0 × 10 6 photons and was assumed to obey the uniform distribution.In MOSE, the 3-D atlas was discritized by triangular elements with 49992 triangles and 24998 surface measurement points, shown in Fig. 2(c).In the following experiments, the default value σ and p used in prior model were set 0.1 and 1.1, respectively.The stopping thresholds ε Φ , K max and ε log were 1.0 × 10 −8 , 10 and 1.0, respectively.In the experiment of single bioluminescent source, a priori permissible source region can be inferred from the surface light power distribution.Figure 3 shows four views of the surface light power distribution with an angular increment of 90 degrees, and red lines denotes the isoline of light power distribution on the phantom surface.From the Fig. 3, the permissible source region can be inferred as: PS = {(x, y, z)|22 < y < 30, 7 < z < 15, (x, y, z) ∈ lung} when default value σ and p were used, the BLT reconstruction was performed and the procedure was terminated after 3 scans, the computational time was about 4 minutes.The corresponding results are shown in Fig. 4 and the reconstruction result indicates that the location of light source is accurately identified.The reconstruction central position and maximum source density were (22.88, 28.93, 12.87) and 0.234nano − Watts/mm 3 , respectively.The relative error (RE) between the real source and reconstruction result was 1.7%, which was calculated according to RE = S recons − S real /S real .S recons and S real are reconstructed source density and real source density, respectively.
In Table 2, the actual bioluminescent source and the bioluminescent source reconstructed with and without Bayesian approach are compiled together for easy comparison.From the Table 2, it can be seen that when the Bayesian approach is used, the bioluminescent source is more accurately reconstructed in terms of the position and source density.
In order to evaluate the robustness of the algorithm, different values of σ and p were set to perform BLT reconstruction when Gaussian noise with different levels was added to the  synthetic data.The representative results were demonstrated in Fig. 5 and corresponding quantitative results are compiled in Table 3.The results show that the method is fairly robust with respect to both σ and p.Furthermore, the proposed algorithm was evaluated in the case of multiple bioluminescent sources.In the experiments, two light sources were placed in the lung with different edge-toedge distance and 20dB Gaussian noise was added to the synthetic data.The reconstruction results were demonstrated in Fig. 6 and quantitative results were compiled in Table 4.In both cases, the maximum relative error between the actual and reconstructed source densities were 42%, but the reconstructed center positions were around the centers of the actual light sources.The physical phantom measurements were acquired based on an imaging configuration of the free-space BLT with 360 • geometry projections [38,39].The sketch of the system structure is shown in Fig. 7.In the experiment, a cube resinous phantom with 20 mm side length was designed and manufactured.The phantom is made from polyoxymethylene, and one hole of 1.25 mm radius was drilled in the phantom to emplace the light source, as shown in Fig. 7(b).
According to luminescence principle of luminescent light stick, peroxide, ester compound solutions and fluorescent dye were injected into the hole in the phantom after thorough mix, and then the red light whose central wavelength located about 660 nm was emitted due to the chemical reaction of the mixed resolutions.In bioluminescent imaging, a liquid-nitrogen-cooled backilluminated CCD camera (Princeton Instruments VersArray 1300B) and a Micro-Nikkor 55 mm f/2.8 lens are used for data acquisition by multi-view noncontact detection mode in a dark environment.The sensitive CCD can generate 1340 × 1300 pixels and 16 bits dynamic range images with 20μm × 20μm sized pixel even if the detectable photon flux density of bioluminescent source on the phantom surface was very weak, reaching pico −Watts/mm 2 .The collected bioluminescent views need to be transformed from grey-scale pixel values into corresponding numbers in physical units.Therefore, camera calibration is a prerequisite for BLT [40].In order to do this, we used an absolutely calibrated integrating sphere of 8-inches in diameter.A filter and variable attenuator help to select a particular wavelength and to control the light level entering the sphere.For a given wavelength, gray levels are associated with varying intensity values.In our experiment, the wavelength range is about 660 nm, and the calibration formula for the CCD camera is given by the formula where ϕ represents phantom density (nano − Watts/mm 2 ) and pixel the pixel value.In addition, the optical parameters of the phantom were determined by a time-correlated single photon counting (TCSPC) system specifically constructed for the optical properties of the turbid medium [41].The final measured optical parameters of the phantom at the wavelength around 660 nm were given as follows: the absorption coefficient μ a ≈ 0.0002mm −1 and the reduced scattering coefficient μ s = (1 − g)μ s ≈ 1.0693mm −1 [42].

Experimental data acquisition and analysis
The physical phantom containing one source centered at (8.75, 8.75, 10) was mounted on a sample stage in front of the CCD camera and the experimental setup is placed in a totally dark imaging chamber to avoid external disturbance and light leaking.Under the computer control, a motorized rotation stage was used to rotate the phantom for recording the flux density with the CCD camera on the four surfaces of the phantom, as schematically shown in Fig. 7(c).
During each data acquisition, one luminescent view was taken by exposing the camera for 10 seconds.The collected pixel gray levels of each luminescent view were transformed into corresponding light units according to the aforementioned calibration formula Eq. ( 18), and the light distribution of each view is demonstrated in Fig. 8.For BLT reconstruction, a finiteelement mesh was built consisting of 1934 tetrahedron elements and 495 nodes with 225 datum nodes on the phantom surface, as shown in Fig. 9(a).Then, the proposed algorithm was applied to reconstruct the light source distribution in the phantom, the result is presented in Fig. 9(b).The reconstructed source located at (8.03, 8.40, 9.31), the computational time was about 37 seconds.The reconstruction results indicates that the position of the light source is accurately identified.

Discussion and conclusions
We have developed a novel BLT reconstruction algorithm to identify the bioluminescent source distribution based on Bayesian approach.The proposed algorithm was not only evaluated with a 3-D micro-CT atlas of inhomogeneous optical properties but also tested through a physical phantom experiment.All results demonstrate preferred localization and quantification both in source center and power.Our work shows the merits and potential of the Bayesian approach for optical molecular tomography.Taking into account the whole computational framework, there are several novel features which have been applied to the proposed algorithm.First, the Bayesian approach for BLT is firstly proposed.BLT attempts to reconstruct the distribution of the bioluminescent source from boundary measurements.While Bayesian approaches are well suited to the ill-posed problem and it provides a framework to incorporate multiple types of a priori information to reduce the ill-posedness of BLT such as anatomical tissue information, tissue optical parameter information and permissible source information.The algorithm maximizes the log posterior probability with respect to a noise parameter and the unknown source.In each scan, the spectral projected gradient-based optimization algorithm is used.Compared with known EM algorithm for BLT, the reconstruction time reported for EM is about several hours [21], but our algorithm needs much less time than that, about tens or hundreds of seconds.Thus, the proposed algorithm is relative computationally inexpensive.
Our study shows that the Bayesian approach provides an effective framework to reconstruct the distribution of bioluminescent source.Furthermore, the proposed Bayesian formulation can be extended to incorporate spectral character of bioluminescent source.Recently, the fusion of multi-modality is becoming an active research field.The Bayesian framework will present a new avenue to deal with the multi-modality information.
Second, in the Bayesian framework, a generalized adaptive Gaussian Markov random field (GAGMRF) prior model for unknown source density estimation is developed.The prior model introduce a prior information for unknown source, therefore, the ill-posedness of BLT can be further reduced.In addition, the prior model is adaptively determined on the basis of finite element analysis, which makes the proposed algorithm possible to handle a more complex geometrical model such as the real mouse phantom.
Another advantage of the proposed algorithm is its regularization property.For most algorithms existed for BLT, Tikhonov approaches are adopted and the regularization parameter has an important role in reconstruction result.Although regularization parameter can be determined with the well-known L-curve methods, the computational burden is expensive.Furthermore, the regularization parameter is not often selected accurately because of the ill-posedness of BLT.However, our algorithm is iterative and imposes regularization to the BLT problem by setting iteration numbers.
With the development of in vivo molecular imaging for small animals, it is necessary for reconstruction using a real animal-shape model.Furthermore, there is a need to take into consideration the heterogeneous optical properties of biomedical tissues.In our algorithm, a geometrical model is generated from a 3-D micro-CT mouse atlas, which is more accurate to evaluate the proposed algorithm than previous regular models.In addition, in order to avoid inverse crime, a Monte Carlo method is adopted to synthesize the measured data, which further ensures the presented algorithm evaluated properly.
In summary, we have presented a novel BLT reconstruction algorithm based on Bayesian approach and demonstrated its feasibility and potential with a micro-CT atlas.It can provide a preferred performance in view of reconstruction quality.Future work will focus on extending our algorithm to real mouse experiments to evaluate its performance as well as developing adaptive finite element meshes to enable increased resolution without overly increasing the computation burden.Relevant results will be reported later.

1 , Φ meas 2 , 5 ) 6 ) 7 )
...,Φ meas M ] T .Bayesian methods provide a natural framework for incorporating prior information about unknown source density x.The MAP estimate of x given the measurement vector y can be represented as follows: xMAP = arg max x≥0 log p(x|y) = arg max x≥0 {log p(y|x) + log p(x)} (Considering the real physical meaning, x ≥ 0 nonnegative constrain is adopted.Taking into the anatomical tissue information C which is derived from a priori anatomical image, the MAP estimate can be modified as xMAP = arg max x≥0 log p(x|y, C) = arg max x≥0 {log p(y|x, C) + log p(x|C)} (When the permissible source region information PS is used as a priori information to reduce the ill-posedness of BLT, the MAP estimate can be further modified as xMAP = arg max x≥0 log p(x|y, C, PS) = arg max x≥0 {log p(y|x, C, PS) + log p(x|C, PS)} (Where p(y|x, C, PS) is the data likelihood and p(x|C, PS) is the conditional probability density function of x given the tissue information C and permissible source region information PS.Given the forward model Eqs.(1) and (2), the data likelihood is governed mainly by the noise statistics, therefore Eq. (7) reduces to xMAP = arg max x≥0 {log p(y|x) + log p(x|C, PS)} (8)

Fig. 1 .
Fig. 1.The view of mouse phantom.(a) Coronal and transverse image of micro-CT data; (b) Labelling of organs and tissues in the coronal and transverse planes; (c) Dorsal-ventral view of surface rendering of the separated organs and tissues: skin, skeleton, heart, lung and liver; (d) Left lateral view of surface rending of above five tissues.

Fig. 2 .
Fig. 2. (a) Heterogeneous phantom with muscle, bone, heart, liver and lung; (b) The mesh used in the tomographic algorithm; (c) The discretized mesh of the phantom used in MOSE.

Fig. 3 .
Fig. 3. Four views of the phantom surface with an angular of 90 degrees.Red lines denote the isoline of the surface light power.(a) Front view; (b) Right view; (c) Back view and (d) Left view.

Fig. 4 .
Fig. 4. BLT Reconstruction with the proposed algorithm.(a) The three-dimensional rendering of the reconstructed result.(b), (c) and (d) are three different slices of the reconstruction which are selected to illustrate the result in more detail.(c) is through the actual source's center; (b) and (d) are perpendicular to the z-axis direction off the actual source's center about ∓0.5 mm.The black circle denotes the actual source.

Fig. 5 .
Fig. 5. BLT reconstruction with different σ and p.(a) and (b) Noise level was 30dB, and for (c) and (d), (e) and (f) was 20dB and 10dB, respectively.From (a) to (f), p gradually reduced and its value was 2, 1.8, 1.6, 1.4, 1.2, 1.1, respectively.The cross-section is through the actual source's center.The black circle denotes the actual source.

Fig. 6 .
Fig. 6.Reconstruction results with two sources.The green mesh denotes the surface mesh of the 3D mouse phantom and the red sphere is the actual source.(a) BLT reconstruction in case of I, and (b) The BLT reconstruction in case of II.

Fig. 7 .
Fig. 7. Experimental setup and phantom.(a) The sketch of the free-space BLT system; (b) The phantom with one light source and (c) The middle cross-section of the phantom.The four degrees show the direction for data acquisition with the CCD camera.The black circle denotes the actual source.

Table 2 .
Quantitative reconstructed results with and without Bayesian approach.

Table 3 .
Quantitative results with different σ and p.

Table 4 .
Quantitative BLT reconstruction results in case of two sources