Image reconstruction for bioluminescence tomography from partial measurement

The bioluminescence tomography is a novel molecular imaging technology for small animal studies. Known reconstruction methods require the completely measured data on the external surface, although only partially measured data is available in practice. In this work, we formulate a mathematical model for BLT from partial data and generalize our previous results on the solution uniqueness to the partial data case. Then we extend two of our reconstruction methods for BLT to this case. The first method is a variant of the well-known EM algorithm. The second one is based on the Landweber scheme. Both methods allow the incorporation of knowledgebased constraints. Two practical constraints, the source non-negativity and support constraints, are introduced to regularize the BLT problem and produce stability. The initial choice of both methods and its influence on the regularization and stability are also discussed. The proposed algorithms are evaluated and validated with intensive numerical simulation and a physical phantom experiment. Quantitative results including the location and source power accuracy are reported. Various algorithmic issues are investigated, especially how to avoid the inverse crime in numerical simulations. © 2007 Optical Society of America OCIS codes: (170.0110) Imaging systems; (170.6960) Tomography; (170.3010) Image reconstruction techniques; (170.3660) Light propagation in tissues. References and links 1. C. Contag and M. H. Bachmann, “Advances in bioluminescence imaging of gene expression,” Annu. Rev. Biomed. Eng. 4, 235 – 260 (2002). 2. V. Ntziachristos, J. Ripoll, L. H. V. Wang, and R. Weissleder, “Looking and listening to light: the evolution of whole-body photonic imaging,” Nat. Biotech. 23, 313 – 320 (2005). 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. 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). 5. Z. Paroo, R. A. Bollinger, D. A. Braasch, E. Richer, D. R. Corey, P. P. Antich, and R. P. Mason, “Validating bioluminescence imaging as a high-throughput, quantitative modality for assessing tumor burden,” Molecular Imaging 3, 117–124 (2004). #81953 $15.00 USD Received 9 Apr 2007; revised 17 Jul 2007; accepted 16 Aug 2007; published 20 Aug 2007 (C) 2007 OSA 3 September 2007 / Vol. 15, No. 18 / OPTICS EXPRESS 11095 6. A. Rehemtulla, L. D. Stegman, S. J. Cardozo, S. Gupta, D. E. Hall, C. H. Contag, and B. D. Ross, “Rapid and quantitative assessment of cancer treatment response using in vivo bioluminescence imaging,” Neoplasia 2, 491 – 495 (2002). 7. A. McCaffrey, M. A. Kay, and C. H. Contag, “Advancing molecular therapies through in vivo bioluminescent imaging,” Molecuar Imaging 2, 75 – 86 (2003). 8. A. Soling and N. G. Rainov, “Bioluminescence imaging in vivo application to cancer research,” Expert Opinion on Biological Therapy 3, 1163 – 1172 (2003). 9. J. C. Wu, I. Y. Chen, G. Sundaresan, J. J. Min, A. De, J. H. Qiao, M. C. Fishbein, and S. S. Gambhir, “Molecular imaging of cardiac cell transplantation in living animals using optical bioluminescence and positron emission tomography,” Circulation 108, 1302 – 1305 (2003). 10. C. H. Contag and B. D. Ross, “It’s not just about anatomy: in vivo bioluminescence imaging as an eyepiece into biology,” J. Magn. Reson. 16, 378 – 387 (2002). 11. G. Wang, E. A. Hoffman, and G. McLennan, “Bioluminescent CT method and apparatus,” (2003). US provisional patent application. 12. G. Wang et al, “Development of the first bioluminescent tomography system,” Radiology Suppl. (Proceedings of the RSNA) 229(P) (2003). 13. G. Wang, Y. Li, and M. Jiang, “Uniqueness theorems for bioluminescent tomography,” Med. Phys. 31, 2289 – 2299 (2004). 14. M. Jiang and G. Wang, “Image reconstruction for bioluminescence tomography,” in “Proceedings of SPIE: Developments in X-Ray Tomography IV,” , vol. 5535 (2004), vol. 5535, pp. 335 – 351. Invited talk. 15. M. Jiang and G. Wang, “Image reconstruction for bioluminescence tomography,” in “Proceedings of the RSNA,” (2004). 16. 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 the Monte Carlo method,” Academic Radiology 11, 1029 – 1038 (2004). 17. X. J. Gu, Q. H. Zhang, L. Larcom, and H. B. Jiang, “Three-dimensional bioluminescence tomography with model-based reconstruction,” Opt. Express 12, 3996–4000 (2004). 18. M. Jiang, T. Zhou, J. T. Cheng, W. Cong, K. Durairaj, and G. Wang, “Image reconstruction for bioluminescence tomography,” in “Proceedings of the RSNA,” (2005). 19. 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, 6756–6771 (2005). 20. A. Cong and G. Wang, “A finite-element-based reconstruction method for 3D fluorescence tomography,” Opt. Express 13, 9847–9857 (2005). 21. C. Kuo, O. Coquoz, T. Troy, N. Zhang, D. Zwarg, and B. Rice, “Bioluminescent tomography for in vivo localization and quantification of luminescent sources from a multiple-view imaging system,” in “SMI Fourth Conference,” (Cologne, Germany, 2005). 22. 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). 23. 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). 24. H. Dehghani, S. Davis, S. D. Jiang, B. Pogue, K. Paulsen, and M. Patterson, “Spectrally resolved bioluminescence optical tomography,” Optics Letters 31, 365 – 367 (2005). 25. S. R. Arridge, “Optical tomography in medical imaging,” Inverse Problems 15, R41 – R93 (1999). 26. F. Natterer and F. Wübbeling, Mathematical Methods in Image Reconstruction (SIAM, Philadelphia, PA, 2001). 27. A. P. Gibson, J. C. Hebden, and S. R. Arridge, “Recent advances in diffuse optical imaging,” Phys. Med. Biol. 50, R1–R43 (2005). 28. A. Cong and G. Wang, “Multispectral bioluminescence tomography: Methodology and simulation,” International Journal of Biomedical Imaging 2006 (2006). Article ID 57614. doi:10.1155/IJBI/2006/57614. 29. C. Q. Li and H. B. Jiang, “Imaging of particle size and concentration in heterogeneous turbid media with multispectral diffuse optical tomography,” Opt. Express 12, 6313–6318 (2004). 30. A. Kak and M. Slaney, Principles of Computerized Tomographic Imaging (IEEE Press, New York, 1987). 31. F. Natterer, The Mathematics of Computerized Tomography (SIAM, Philadelphia, PA, 2001). 32. A. P. Dempster, N. M. Laird, and D. B. Rubin, “Maximal likelihood form incomplete data via the EM algorithm,” Journal of the Royal Statistical Society. Series B. 39, 1 – 38 (1977). 33. L. A. Shepp and Y. Vardi, “Maximum likelihood restoration for emission tomography,” IEEE Transactions on Medical Imaging 1, 113 – 122 (1982). 34. D. L. Snyder, T. J. Schulz, and J. A. O’Sullivan, “Deblurring subject to nonnegativity constraints,” IEEE Transactions on Signal Processing 40, 1143 – 1150 (1992). 35. M. Jiang and G. Wang, “Convergence studies on iterative algorithms for image reconstruction,” IEEE Transactions on Medical Imaging 22, 569 – 579 (2003). #81953 $15.00 USD Received 9 Apr 2007; revised 17 Jul 2007; accepted 16 Aug 2007; published 20 Aug 2007 (C) 2007 OSA 3 September 2007 / Vol. 15, No. 18 / OPTICS EXPRESS 11096 36. M. Jiang and G. Wang, “Development of iterative algorithms for image reconstruction,” J. X-Ray Sci. Technol. 10, 77 – 86 (2002). Invited Review. 37. M. Piana and M. Bertero, “Projected Landweber method and preconditioning,” Inverse Problems 13, 441 – 463 (1997). 38. A. Sabharwal and L. C. Potter, “Convexly constrained linear inverse problems: iterative leat-squares and regularization,” IEEE Transactions on Signal Processing 46, 2345 – 2352 (1998). 39. A. Ishimaru, Wave Propagation and Scattering in Random Media (IEEE Press, New York, 1997). 40. A. D. Klose and A. H. Hielscher, “Quasi-Newton methods in optical tomographic image reconstruction,” Inverse Problems 19, 387–409 (2003). 41. D. S. Anikonov, A. E. Kovtanyuk, and I. V. Prokhorov, Transport equation and tomography, Inverse and Ill-posed Problems Series (VSP, Utrecht, 2002). 42. D. Gilbarg and N. S. Trudinger, Elliptic Partial Differential Equations of Second Order, vol. 224 of Grundlehren der mathematischen Wissenschaften (Springer–Verlag, Berlin–Heideberg–New York, 1983). 43. R. Dautray and J. L. Lions, Mathematical Analysis and Numerical Methods for Science and Technology, vol. I (Springer-Verlag, Berlin, 1990). 44. V. Isakov, Inverse Problems for Partial Differential Equations, vol. 127 of Applied Mathematical Series (Springer, New York–Berlin–Heideberg, 1998). 45. W. Rudin, Functional analysis, International Series in Pure and Applied Mathematics (McGraw-Hill, New York, 1991), 2nd ed. 46. M. H. Protter and H. F. Weinberger, Maximum Principles in Differential Equations (Prentice-Hall, Englewood Cliffs, N. J., 1967). 47. B. Eicke, “Konvex-resringierte schlechtgestellte Problems und ihr Regularisierung durch Iterationsverfahren,” Thesis, Technischen Universität Berlin (1991). 48. B. Eicke, “Iteration methods for convexly constrained ill-posed problems in Hilbert space,” Numerical Functional Analysis and Optimization 13, 413 – 429 (1992). 49. S. C. Brenner and L. R. Scott, The mathematical theory of finite element methods, Texts in applied mathematics; 15 (Springer-Verlag, New York, NY, 2002), 2nd ed. 50. D. L. Colton and R. Kress, Inverse acoustic and elctromagnetic scattering theory (Sp


Introduction
Gene therapy is a breakthrough in the modern medicine, which promises to cure diseases by modifying gene expressions.A key for the development of gene therapy is to monitor the in vivo gene transfer and its efficacy in the mouse model.Traditional biopsy methods are invasive, insensitive, inaccurate, inefficient, and limited in the extent.To map the distribution of the administered gene, reporter genes such as those producing luciferase are used to generate light signals within a living mouse.Bioluminescent imaging (BLI) is an optical technique for sensing gene expression, protein functions and other biological processes in mouse models by reporters such as luciferases that generate internal biological light sources [1,2].The light emitted within the mouse can be captured externally using a highly sensitive CCD camera [3]."The collection times for BLI are relatively short compared to non-optical modalities, an advantage of optical imaging methods in general" [1].
BLI has great potentials in various biomedical applications, including regenerative medicine, developmental therapeutics, treatment of residual minimal disease, and studies on cancer stem cells."Additionally, time-lapse whole body imaging of animals bearing xenografts of appropriately labeled cancer cells can provide new information on tumor growth dynamics and metastasis patterns that it is not possible to obtain by invasive experimental approaches" [4]."Bioluminescence imaging (BLI) is a highly sensitive tool for visualizing tumors, neoplastic development, metastatic spread, and response to therapy" [5].Although the spatial resolution is limited when compared with other imaging modalities, the advantages of BLI in vivo imaging are sensitivity, speed, non-invasiveness, cost, ease, and low background noise (in contrast to fluorescence, PET and other non-optical modalities) [1].BLI could be applied to study almost all diseases in small animal models [6,1,7,8,9,10,2].
However, the current BLI technology has not explored the full potential of this approach.It only works in 2D imaging mode and is incapable of 3D imaging of the light source location associated with specific organs and tissues [2].Since its first introduction in 2003 [11,12], the bioluminescence tomography (BLT) has been undergone a rapid development [11,12,13,14,15,16,17,18,19,20,21,4,22,2,23,24].BLT is to address the needs for 3D localization and quantification of an internal bioluminescent source distribution in a small animal.With BLT, optimal analyzes on a bioluminescent source distribution become feasible inside a living mouse.Although there are currently several approaches for this technique [11,12,13,14,15,16,17,18,19,20,21,4,22,2,23,24], the originally proposed approach in [11,12] is still one of the promising techniques in this field.With this technique, the complete knowledge on the optical properties of anatomical components is assumed to be available from an independent tomography scan, such as a CT/micro-CT, MRI scan and/or diffuse optical tomography (DOT), by image segmentation and optical property mapping.That is, we can segment the image volume into a number of anatomical structures, and assign optical properties to each component using a database of the optical properties, or use the DOT technique for this purpose [11,12].
Traditionally, optical tomography sends visible or near infra-red light to probe a scattering object, and reconstructs the distribution of internal optical properties, such as absorption and scattering coefficients [25,26,27].In contrast to this active imaging modality, BLT reconstructs an internal bioluminescent source distribution generated by luciferase activated with reporter genes, from BLI measurements at the object surface.Mathematically, BLT is a source inversion problem based on the boundary measurement, and hence is a highly ill-posed inverse problem per se.To obtain satisfactory results for the BLT, prior knowledge must be utilized to regularize the problem.The tomographic feasibility and the solution uniqueness have been theoretically studied in [13].It was proved that the uniqueness of the BLT does not hold in general.In our previous studies, we utilized constraints such as the non-negativity and source support [14,15,18,19].Other constraints such as the range of the source intensity may be effective as well.Another approach is to utilize spectral resolved measurement such as the hyper-spectral or multi-spectral measurement [4,22,24,28] and multi-spectral source information [29] to improve the BLT reconstruction.
In previous studies, complete data measured on the full object boundary is required to reconstruct an internal bioluminescent source distribution.In practice, the measured data for BLT is often incomplete due to physical limitations as in X-ray CT [30,31,26].The BLT in this situation is similar to CT image reconstruction from angle-limited data.In this work, we propose an approach for BLT from incomplete or partial boundary measurement and establish two iterative reconstruction algorithms based on the diffusion approximation.First we generalize the results on the solution uniqueness in [13] to the partial measurement case.It is proved that similar results still hold given practical optical signals on the mouse body surface, although the solution characterization is more complicated.A new important formulation is the definition of the Dirichlet-to-Neumann (or Steklov-Poincaré) map in (21) for the partial measurement case.Then we extend two of our reconstruction methods for BLT in [14,15,18] to the partial measurement case.The first algorithm is a variant of the well-known expectation-maximization (EM) algorithm [32,33,34].The second one is based on the projected Landweber scheme [35,36,37,38].Both methods allow the incorporation of knowledge-based constraints.Two practical constraints, the source non-negativity and support constraints, are introduced to regularize the BLT problem and produce stability.The initial choice of both methods and its influence on the regularization and stability are also discussed.Both algorithms are evaluated and validated with intensive numerical simulation and a physical phantom.Finally algorithmic issues are investigated, especially the method to avoid the inverse crime in numerical simulations.
The organization of the paper is as follows.We introduce the problem of BLT from partial measurement in § 2, reformulate it in an operator form by generalizing the Dirichlet-to-Neumann map in § 3, and extend our results on the solution uniqueness of BLT in § 4.Then, we present our iterative BLT algorithms in § 6.We report our numerical and experimental results in § 7. Finally we discuss technical problems and research topics with the current BLT techniques in § 8 and conclude the paper in § 9.

Formulation of BLT
Let Ω be a bounded smooth domain in the three-dimensional Euclidean space R 3 that contains an object to be imaged.BLT is to reconstruct the source q from measurement on the boundary of Ω enclosing the source.Let u(x, θ ) be the radiance in direction θ ∈ S 2 at x ∈ Ω, where S 2 is the unit sphere.A general model for light migration in a random medium is the radiative transfer equation (RTE) [39,25,26] (1) for t > 0, and x ∈ Ω, where c denotes the particle speed, μ = μ a + μ s with μ a and μ s being the absorption and scattering coefficients respectively, the scattering kernel η is normalized such that S 2 η(θ • θ ) dθ = 1, and q is the internal light source.In (1), the radiance u(x, θ ,t) is in W cm −2 sr −1 , the source term q(x, θ ,t) in W cm −3 sr −1 , the scattering coefficient μ s and the absorption coefficient μ a both are given in cm −1 , and the scattering phase function η is in sr −1 [40].To find a solution for (1), we need the initial condition and the boundary condition for u where ν is the exterior normal on the boundary Γ of Ω.The forward problem (1), ( 2) and (3) admits a unique solution under natural assumptions on μ, μ a and η [41].The homogeneous condition g − (x, θ ,t) = 0 specifies that no photons travel in an inward direction at the boundary, except for source terms [25, p. R50], which is the case for BLT.
In terms of the radiance u(x, θ ,t), the measurement at a boundary points is given by With the RTE, it is quite complex to reconstruct the light source q from the measurement g and with the above initial-boundary problem being as the forward process.This is largely due to the difficulty in computing the solution u for the forward problem (1), ( 2) and (3).
Typical values of the optical parameters for biological tissues are μ a = 0.1 − 1.0mm −1 , μ s = 100 − 200 mm −1 , respectively.This means that the mean free path of the particles is between 0.005 and 0.01 mm, which is very small compared to a typical object.Thus the predominant phenomenon in optical tomography is scattering rather than transport.Therefore, the diffusion approximation has been widely applied to simplify the RTE (1) in optical tomography [25,26,27].The diffusion approximation assumes that the radiance u(x, θ ,t) is isotropic and takes the average radiance u 0 (x,t) to approximate the radiance u(x, θ ,t) The forward process in the diffusion regime is then given by the following initial-boundary problem [39,25,26] 1 c For diffusion approximation based BLT, the measured quantity at x ∈ Γ for t ≥ 0 is given by Because the internal bioluminescence distribution is relatively stable, we can use the stationary version of ( 9) and (11) as the forward model for BLT.By discarding all the time dependent terms in (9) and (11), the stationary forward model is given by the following boundary value problem (BVP) In practice, it is difficult to obtain all the measurement on the boundary Γ.We consider the case in which the measurement is conducted on some disjoint patches Γ j ⊂ Γ for j = 1, ••• , J, each of which is smooth, closed and connected.Let Γ P = J j=1 Γ j .The BLT problem from partial measurement can be stated as follows: Given the incoming light g − on Γ and outgoing radiance g on Γ P , find a source q 0 with the corresponding diffusion approximation u 0 such that In a typical BLT setting, g − = 0, because there is no incoming light.The optical parameters D and μ a in the above formulations can be established point-wise as mentioned above [11,12].

Reformulation of BLT(P)
The uniqueness property of the BLT problem was already studied [13].For the BLT(P) problem (15), we analyze its solution uniqueness in this section.The reader is reminded that the presentation in its mathematically accurate form will require rather technical and tedious assumptions on the domain and the coefficients.Here we will make the mathematical presentations as precise as possible while keeping a reasonable readability.For details, please refer to [13,14] and the references therein.
In additions to the assumptions for the the BLT(P) problem (15), we assume that the parameters D > D 0 > 0 for some positive constant D 0 and that μ a ≥ 0 are bounded functions.We further assume that D is sufficiently regular near Γ.For example, D is equal to a constant near Γ.Let γ 0 and γ 1 be the boundary value maps and L be the differential operator be the solution of the following mixed boundary value problem (MBVP) [42,43] We define a linear operator On the other hand, for q 0 ∈ L 2 (Ω), we consider the following MBVP and define another linear operator Both operators N Γ P and Λ Γ P are extensions of the well-known Dirichlet-to-Neumann (or Steklov-Poincaré) map ( [44]) and an relevant operator Λ in the complete measurement case ( [14,15,18]) to the partial measurement case, respectively.Assume that q 0 is a solution of the BLT(P) problem (15) with one corresponding radiance diffusion approximation solution u, given the observed g on Γ P and assumed g − on Γ.Then u is the unique solution of the following MBVP On the other hand, let w 1 be defined as in ( 18) -( 20) with f = g − + 2g, w 2 be defined as in ( 22) -( 24), and v = w 1 + w 2 .Then we have By the solution uniqueness of this MBVP [42], it follows that v satisfies ( 26) - (28).Hence, u = v is the required radiance u that generates the measurement on Γ P .The measurement equation ( 12) implies that Hence, q 0 satisfies the following equation Conversely, if there exists a q 0 satisfying (33), we can construct u = v = w 1 + w 2 as above.It follows that u satisfies the forward model and the measurement equation.
In summary we have Proposition 3.1.q 0 is a solution to the (BLT(P)) problem (15) if and only if it is a solution to the equation (33).
This is an extension of the result for the BLT problem in our previous studies.[13,14,15,18]

Solution structure of BLT(P)
In practice the BLT(P) problem is expected to have at least one solution.Therefore, we will not discuss the existence of the BLT(P) problem but focus on the uniqueness of the BLT(P) solution.The uniqueness of the BLT solution was discussed in our previous study [13].In what follows we will demonstrate how to extend the results to the BLT(P) problem.We need the following notations from functional analysis [45].Let A be a linear operator from a Banach space X to a Banach space Y .The kernel or null space of A is defined as for some x ∈ X}.For a subspace M of a Hilbert space H, M ⊥ is the set of all y ∈ H, such that y, x = 0 for all x ∈ M.
By Theorem 3.1, the uniqueness of BLT(P) solution can be characterized by determining the kernel We need the Green's formula [42,43] For ψ ∈ H 1 2 (Γ P ), let T Γ P be defined by as the unique solution in H 1 (Ω) ⊂ L 2 (Ω) of the boundary problem Then, we can prove with the Green's formula that for the operators i.e., they are adjoint to each other, Hence, the kernel of To obtain an explicit characterization of N [Λ Γ P ], let

R[T Γ
, by the Green's formula (34), We have, by ( 22) -( 25), there exists w 2 such that We have w 2 ∈ H 2 (Ω) by the regularity theory for second order elliptic partial differential equations [42,43].The above boundary conditions imply that w 2 ∈ H 2 0,Γ P (Ω).Hence, The conclusion follows immediately.By Proposition 3.1, the BLT(P) problem is equivalent to the linear equation (33) with q 0 as the unknown to be found.All the solutions q 0 to (33) form a convex set in L 2 (Ω).Hence, there exists one unique solution of the minimal L 2 -norm among those solutions [45].Let this minimal norm solution be denoted as q H .Then, all the solutions can be expressed as We summarize the above results into the following theorem.
Theorem 4.2.Assume that the BLT(P) problem is solvable.For any couple (g − , g) such that there is one special solution q H for the BLT(P) problem (15), which is of the minimal L 2 -norm among all the solutions.Then, any solution can be expressed as q 0 = q H + L[p], for some p ∈ H 2 0,Γ P (Ω).

Uniqueness of the BLT(P) solution in RBF
It is remarkable that the uniqueness results in [13] for sources consisting of radial base functions (RBF) still hold for the BLT(P) problem after examining the assumptions and proofs therein.
We assume the following conditions on Ω, D, μ a and q 0 .C1: Ω is a bounded C 2 domain of R 3 and partitioned into non-overlapping sub-domains Ω i , i = 1, 2, ..., I; C2: Each Ω i is connected with piecewise C 2 boundary Γ i ; C3: D and μ a are C 2 near the boundary of each sub-domain.
C4: D > D 0 > 0 for some positive constant D 0 is Lipschitz on each sub-domain; μ a ≥ 0 and μ a ∈ L p (Ω) for some p > N/2; The first uniqueness result is for sources consisting of bioluminescent point impulses where each a s is a constant, and y s the location of a point source inside some Theorem 5.1.( [13]) Assume the conditions C1 -C4 hold.If q 0 (y) = ∑ S s=1 a s δ (y − y s ) and q 0 (y) = ∑ S s=1 a s δ (y − y s ) are two solutions to the BLT(P) problem (15) with the same measurement on Γ P , then S = S and there is a permutation τ of {1, ••• , S} such that a s = a τ(s) and y s = y τ(s) .
The second uniqueness result is for sources as a linear combination of RBFs, which includes the previous type of sources in (45) as the limiting case, where each g s ∈ L 2 (B r s 0 ,r s 1 (x s )) is continuous, the source centers {x s } are distinct, and each source support B r s 0 ,r s 1 (x s ) ⊂⊂ Ω i for some 1 ≤ i ≤ I.For each 0 ≤ r 0 < r 1 < ∞, x 0 ∈ R 3 , B r 0 ,r 1 (x 0 ) denotes a hollow ball specified by r 0 < ||x − x 0 || < r 1 for r 0 > 0 or a solid ball specified by ||x − x 0 || < r 1 for r 0 = 0. Let ϕ(r) be defined as follows.For μ a = 0, and for μ a > 0, For the second result we need to replace the condition C4 with the following condition C4* and a new condition C5: Note that the condition C4* is a special case of the condition C4.The second uniqueness result can be stated as follows.
Theorem 5.2.( [13]) Assume the conditions C1 -C3, C4* and C5 hold.If q 1 (y) = ∑ S s=1 g s (||y− y s ||)χ B r s (X s ) are two solutions to the BLT(P) problem (15) with the same measurement on Γ P , then S = S and there exist a permutation τ of {1, ••• , S} and a map C : where ϕ j is given by ( 48) or (47) with D = D j and μ a = μ j .

Reconstruction methods
In our previous studies, the BLT problem with complete measurement is reformulated as a linear equation by the Dirichlet-to-Neumann map and solved with the proposed EM-like and constrained Landweber (CL) iterative algorithms therein.In this section we follow the same approach to establish reconstruction methods for the BLT(P) problem by solving the linear equation (33).
. By (33), the BLT solution q 0 satisfies the following equation The extended EM and CL iterative algorithms for solving the BLT(P) problem ( 50) are established in the following.

Method I: EM method
Based on the formulation in [26,14], let which is the log likelihood function when the measured data b is subject to Poissonian distribution.By the maximum principle of elliptic partial differential equations [42], it follows that Λ Γ P [q 0 ] ≥ 0 and b ≥ 0 if q 0 ≥ 0, q 0 = 0, g ≥ 0, g = 0 and g − ≥ 0. We try to find a solution for the BLT(P) problem by performing the following optimization arg max We first assume that q 0 > 0 is a minimizer of F. The case of q 0 ≥ 0 can be handled similarly as the limiting case.We need to find the Frechet derivative of F. Let where v is an arbitrary bounded function of L 2 (Ω), and compute Hence, the Frechet derivative of F is If q 0 > 0 is a solution of ( 52), it follows that F [q 0 ] = 0.The general case of q 0 ≥ 0 is given by the following Kuhn-Tucker condition [26] Let φ 1 = Λ * Γ P [1] = T Γ P [1], i.e., the solution of the following MBVP, by ( 40) and ( 36) - (38), It follows from the maximum principle of elliptic partial differential equations that 0 < φ 1 ≤ 1 [46].Because Λ * Γ P = T Γ P by ( 40), the Kuhn-Tucker condition (55) can be rewritten as Then we obtain the following EM formula for BLT from partial measurement . (60)

Method II: CL method
Assume that we have some prior knowledge about the source represented as a convex set C = {q 0 : q 0 satisfies some convex constraints.}, which is a closed convex subset of L 2 (Ω).Let P C be the orthogonal projection operator from L 2 (Ω) to C .Then the CL scheme, or projected Landweber scheme, for solving (50), is given as follows q (n+1) 0 where λ n is a relaxation parameter and Λ * Γ P = T Γ P by (40).The convergence property of the CL scheme was studied in [47,48] and improved in [37,38].Conditions for the relaxation parameter depend on the operator norm ||Λ * Γ P Λ Γ P || = ||Λ Γ P || 2 .For example, 0 < λ n < 2 ||Λ Γ P || 2 [37,38].The limit of q (n) 0 is a solution to the constrained least-squares problem arg min In this work, the non-negativity and source support constraints are applied.The non-negativity constraint implies that the source is non-negative and is utilized by setting the negative parts of iterated sources to be zero.The source support constraint assumes that the source is non-zero only within some sub-region Ω 0 of the region Ω, and is utilized by setting the iterated sources to be zero outside Ω 0 .The choice of Ω 0 is discussed in § 6.In all iterative image restoration methods, we have to initialize the process.We propose the following method to choose an initial guess q (0) 0 .By Green's formula, let v = u and p = 1 in (34), we have, If we replace q 0 with q (0) 0 , we obtain Ω q (0) Because u ≥ w 1 by the the maximum principle of elliptic partial differential equations [42], we have Ω q (0) In practice, the source q 0 is compactly supported on a subset Ω 0 inside Ω.We choose q 0 such that it is equal to a positive constant in its support Ω 0 and zero otherwise.Hence where Q 0 is a constant, and χ Ω 0 (x) = 1 on Ω 0 and is zero otherwise.The support Ω 0 of the source q 0 is part of the prior knowledge, which was termed the permissible region in and could be inferred from the measured data [19].Please refer to § 8 for further discussions.Because only the data g on Γ P is available, the final estimate is where |Ω 0 | is the volume of Ω 0 .

Convergence criteria
The convergence criteria for both algorithms may include (1) when the iteration number n reaches an assumed maximum number; (2) when the successive incremental |q 0 | is smaller than an assumed error level.In this work, the convergence criterion is by manually setting the iteration numbers to fixed numbers for both methods, respectively.

Numerical experiments
For inverse problems, numerical tests of reconstruction methods usually make use of simulated data from the numerical solution of the forward problem.One typical issue is coined as the inverse crimes in [50].This happens especially when insufficient rough discretization or the same discretization are used for the forward and inverse process, because "it is possible that the essential ill-posedness of the inverse problem may not be evident" [50, p. 304].Hence, the results could be overly optimistic and unreliable."Unfortunately, not all of the numerical reconstructions which have appeared in the literature meet with this obvious requirement" [50, p. 133].As suggested, "it is crucial that the synthetic data be obtained by a forward solver which has no connection to the inverse solver under consideration" [50, p. 133].
One approach to avoid the inverse crime is to use different discretizations in the forward and the inverse process [50, p. 304].Because our formulation of the reconstruction methods is analytical and independent of any discretization, we can use different finite elements (shape functions and meshes) in Comsol TM for the forward process and reconstruction algorithms, respectively.Moreover, we can change the mesh sizes with the adaptive mesh technique at each iteration step for the reconstruction algorithms to solve the MBVPs ( 22) -( 24) and ( 36) - (38).During the iteration intermediate results at different meshes are interpolated to the required nodes with the built-in bilinear interpolation method in Comsol TM , when values at the nodes are required.
We have carries out intensive numerical simulation for various numerical phantoms.Various algorithmic settings such as the initial choices, iteration numbers, shape functions and meshes of the FEM have also been evaluated.Because the algorithms depend on multiple factors, the optimal setting is still under investigation.Nevertheless both methods can reliably reconstruct the sources in most settings.Due to the limited space, we report one representative result with the CL method in the following.
In this experiment, a numerical phantom with the same geometrical structure as the physical phantom in § 7.2 is used.This is a cylindrical heterogeneous mouse chest phantom (Fig. 1(a)) of 30mm height and 30mm diameter.Its structure is shown in Fig. 1(b).Three sources of the form for i = 1, 2, 3, are set up in this simulation, where Ω i is a ball centered at x i : The radius r i are all set to 1mm.The sources are centered at x 1 = (-0.9cm,0.25cm, 0cm), x 2 = (-0.9cm,-0.25cm, 0cm) and x 3 = (0.9cm, 0.25cm, 0cm), with the intensity values A 1 = 25.1μW/cm 3 , A 2 = 23.3μW/cm 3 and A 3 = 25.1μW/cm 3 , respectively.These intensity values are set according to the total source power of the physical phantom in § 7.2 so that the total source power of each source is equal to one of the physical sources in § 7.2.The optical coefficients of the phantom are set as in Table 1.As discussed above, we use different finite elements (shape functions and meshes) for the forward process and inverse process, respectively.The information of the finite elements is  Figure 2 shows the results obtained with the CL method.The initial support or the permissible region is set to Q 0 = {(x, y, z) : 0.8 < (x 2 + y 2 ) 1/2 < 1.2, −0.15 < z < 0.15}.The relaxation coefficient λ is manually set to λ = 20.The computational overhead for the cases of complete measurement and partial measurement is about the same.It takes about 6 hours for 70 iterations to get the results in Fig. 2. The figures are generated with the commands postplot, geomplot and meshplot of Comsol TM with manually adjusted parameters at different views, respectively.
In the case of the partial measurement from the front view, the source in the right lung close to the back view is not reconstructed, see Fig. 2 (c) and (d).This is reasonable because that source is far from the measurement surface in terms of the mean-free path.When complete measured data is used for reconstruction, this source could be reliably reconstructed and is shown in Fig. 2 (a) and (b).Quantitative results about the location accuracy are compiled into Table 3.The absolute error is defined by is the reconstructed center of each source and is estimated interactively from the reconstructed source distribution from different views.The relative error is defined by the absolute error divided by The results demonstrate that the same location accuracy for the left two sources can be obtained with only the partial measurement in the front view.Table 3. Quantitative results for the reconstructed locations of the three sources at S1 = (-0.90,0.25, 0), S2 = (-0.90,-0.25, 0) and S3 = (0.90, 0.25, 0), respectively.The unit is cm.

Measurement
Reconstructed As reported in the literature, the source power was estimated as the source integral q(x) dx of the source intensity over its support [19,51].Another approach for the source power estimation for a RBF source is based on the result of Theorem 5.2 by computing the source moment ϕ s (||x The reconstructed source intensities and sizes are estimated interactively from orthogonal views of the reconstructed source distributions.Table 4 presents the results of the reconstructed source integrals and its absolute and relative errors with the true value, respectively, for each source.Table 5 shows the corresponding results for the source moments.The differences of the results in both tables are discussed in § 8.

Physical experiment
We use the same physical phantom in our previous study [19].As described in § 7.1, a cylindrical heterogeneous mouse chest phantom (Fig. 1 A luminescent light stick (Glowproducts, Canada) was selected as the testing source.The stick consisted of a glass vial containing one chemical solution and a larger plastic vial containing another solution with the former being embedded in the latter.By bending the plastic vial, the glass vial can be broken to mix the two solutions after which luminescent light was emitted.The particular dye in the chemical solution was for red light, and it could last for approximately 4 hours at an emission wavelength range between 650nm and 700nm, being close to that of the red spectral region of the luciferase.Two small holes of diameter 0.6mm and height 3mm were drilled in the phantom with their centers at (-0.9cm, 0.15cm, 0.0cm) and (-0.9cm, -0.15cm, 0.0cm) in the left lung region of the phantom, respectively.Two red luminescent liquid filled catheter tubes of 1.9mm height and 0.56mm diameter were placed inside the two small holes, respectively.We measured the total power of the red luminescent liquid filled polythene tubes with the CCD camera.They were 105.1 nW and 97.4 nW, respectively.In our bioluminescent imaging, a CCD camera (Princeton Instruments VA 1300B, Roper Scientific, Trenton, NJ) was used for data acquisition on the phantom surface.The collected bioluminescent views were transformed from grey-scale pixel values into equivalent numbers in physical units.We used the same calibration approach described in the previous study [19].The measured data on the CCD were transformed to the form of the measurement equation ( 12) through a geometric optics mapping of light beams.The required optical parameters of the phantom in Table 1 were the same as in the previous study [19].
Due to the limited space, representative results from the physical phantom are given in Fig. 4(a) and Fig. 4(b) using the EM algorithm from only the data measured in the four views of the side surface of the phantom.We conducted experiments for reconstruction from the data measured in one or several of the four views.Figure 4(c) and 4(d) are the results reconstructed using the EM algorithm from only the data measured in the front view.The results from the measured data in other single views or combinations of these three views were not encouraging, because the signal-to-noise ratios were too low.For the results in Fig. 4 (a) -(d), the iteration number for the EM algorithm was set to 50 along with an initial source support region, Q 0 = {(x, y, z) : 0.8 < (x 2 + y 2 ) 1/2 < 1.2, −0.15 < z < 0.15}.It takes about 5 hours for the 50 iterations of the EM algorithm.For the CL algorithm, similar results were obtained but the separation of the two sources was not as good as that obtained using the EM algorithm.Quantitative results about the location accuracy and source power computed by the source integral are summarized in Table 6 and Table 7.The source powers estimated by source moments are not available because the original source is not of RBF sources, please refer to discussions in § 8.The reconstructed source intensities and sizes are estimated interactively from orthogonal views of the reconstructed source distributions.The information of the finite element in this

Discussions
Because there is no unique solution to the BLT(P) problem in the general case by Theorem 4.2, one may consider to utilize the minimal norm solution q H as the solution of the BLT(P) problem.The minimal norm solution q H is unique and also called the minimal energy solution, and advocated in other applications [38,52].However, it can be similarly demonstrated as in our previous study [14] that the minimal energy source solution is not favorable for the BLT(P) problem in general.Because sources of compact supports are commonly encountered in practice, it can be proved in the same way that such sources cannot be found as the minimal norm solution for the BLT(P) problem [14].It has been reported that adequate prior knowledge must be utilized to obtain a physically favorable BLT solution [11,12,13,14,15,16,18,19,20].
From the physical perspective, the multi-spectral technique is also a promising approach to resolve this issue, although it needs extra hardware and exposure time at a compromised signal-to-noise ratio [4,22,24,28,29].The methods for BLT from partial measured data in this work can be combined with the multi-spectral technique as well.
The iterative approach provides a mechanism for incorporating prior knowledge based constraints and has been widely used in practice.In this paper we have established two iterative algorithms for BLT from partial measurement.In this work, two constraints, the non-negativity and source support constraints, are applied.The EM algorithm induces the non-negativity of the source because of the maximum principle of elliptic partial differential equations, if the initial choice is non-negative.The source support with the EM algorithm is automatically implied because the support will not increase during the iteration due to its multiplication operation in (60).Hence, it seems that the EM algorithm is more preferable than the CL algorithm with those two constraints.However, the CL algorithm is more flexible to add other constraints and has the established convergence property.Other constraints highlighted by Theorem 5.1 and 5.2 are to restrict the solution space to a sub-space of bioluminescent source distributions so that the solution uniqueness holds to a practical extent.Nevertheless, the source support constraint helps to resolve the non-uniqueness because of the property of both the EM and CL algorithms.
Another advantage of the proposed iterative methods for BLT is its regularization property.For inverse problems, popular approaches include the Tikhonov [53] and the iterative regularization methods [54].For the latter the regularization is achieved by setting the iteration number.Both of our proposed methods are iterative and impose regularization to the BLT problem by setting iteration numbers.It was proved in [55] that the CL and Tikhonov regularization methods are equivalent.As discussed in the previous paragraph, the source support of the initial choice also helps to regularize the problem and produce stability.
There are remaining mathematical issues with the proposed algorithms.For the EM algorithm, its convergence has not been established, though it converges in our experiments.For the CL method, there is no guide in choosing the relaxation coefficient λ n .This parameter depends on the operator norm ||Λ * Γ P Λ Γ P || = ||Λ Γ P || 2 , which is equivalent to find the minimal eigenvalue of Λ * Γ P Λ Γ P or the minimal singular value of Λ Γ P .This can be reduced to a boundary eigenvalue problem of partial differential equations.Both problems are left for future investigation.
There are also remaining physical issues with the proposed approach.We have used a geometric optics mapping of light beams to transform the measured data by a CCD detector to the surface measurement equation (12).More advanced techniques should be used to improve this process, such as the recently proposed non-contact measurement technique for fluorescence tomography [56].Furthermore, the diffusion approximation can be improved by the radiative transfer equation [51].
There are two methods in quantifying the source power, namely, by the source integral and source moment.By (49), the source integral is equal to the source moment when the absorption coefficient μ a = 0.As presented in Table 4 and 5, the error in source power computed by the source moment is smaller than that by the source integral in numerical simulations.This is in accordance with the theoretical result in Theorem 5.2.The results of source moments presented in this work are based on estimated source sizes interactively from orthogonal views of the reconstructed source distributions.More sophisticated techniques of RBF approximations should be applied to get more accurate estimations [57].Although estimating source powers by the source moment is more accurate than by the source integral, it is unclear how the source moment is related to the total biological activities and how to extend it to a general source without the RBF structure.Both methods of computing source powers will be investigated in real biological experiments in our future work.
The computational overhead reported in this work is about 5 -6 hours.Most of the computation is spent on solving MBVPs, for which there are sophisticated parallel algorithms.We are working to improve the computing performance with the Linux cluster technology, the most

Conclusions
To overcome the limitations in the previous BLT studies, we have formulated a mathematical framework for BLT from partial boundary measurement, extended our results on the solution uniqueness, and proposed the two iterative reconstruction algorithms based on the diffusion approximation.Either of the methods is suitable for incorporating practical constraints.Two practical constraints, the source non-negativity and support constraints, have been introduced to regularize the BLT problem and produce stability.The initial choice of both methods and its influence on the regularization and stability have also been discussed.Numerical and physical phantoms have been used to evaluate and validate the proposed algorithms with quantitative results.Technical problems and research topics with the BLT technique have also been discussed.

Fig. 1 .
Fig. 1.(a) A heterogeneous mouse phantom consisting of bone (B), heart (H), lungs (L), and muscle (M).(b) A cross-section through two luminescent sources in the left lung and another source in the right lung.The four arrows show the directions of the CCD camera for data measurement.

Fig. 2 .
Fig. 2. Reconstructed results by the CL method and a cross-section at z = 0cm.(a) and (b) are results from data measured at the four views.(c) and (d) are from data measured at the front view only.

Fig. 3 .
Fig. 3. (a) A cross-section through two hollow cylinders for hosting luminescent sources in the left lung.The four arrows show the direction of the CCD camera during data acquisition.(b) The measurements at the four views combined along the phantom side surface with unit μW/cm 2

Fig. 4 .
Fig. 4. Representative results reconstructed by the EM algorithm.(a) and (c) are the sources reconstructed by the EM algorithm from the data measured in the four views and in the front view only, respectively.(b) and (d) are cross-sections at z = 0 cm of the sources in (a) and (c), respectively.
(20)on feature of the proposed EM and CL methods is that they are of an iterative nature.At each step, they require the evaluation of the operators Λ Γ P by(25)based on the MBVP (22) -(24) and T Γ P by(35)based on the MBVP (36) -(38).Both MBVPs are of the same type can be solved with the finite element method (FEM)[49].In this work, both MBVPs are solved with the FEM software Comsol TM and Matlab TM .Note b in (50) for both methods can be computed by(21)after solving the MBVP (18) -(20).The computer is a Dell TM workstation, Precision 670, with dual Intel TM Xeon CPUs of main frequency 2.8GHz and 6GB memory.The operating system is Microsoft TM Windows XP Professional X64 Edition.Other details are reported in § 7. A

Table 1 .
Optical parameters for the phantom

Table 2 .
Finite element information for the simulation

Table 4 .
Quantitative results for the reconstructed source integrals of the sources.The sources are listed in the order as in Table3.Their true values are 105.1,97.4 and 105.1, respectively.The unit is nW.

Table 5 .
Quantitative results for the reconstructed source moments of the sources.The sources are listed in the order as in Table3.Their true values are 125.5, 116.5 and 125.5, respectively.The unit is nW.

Table 7 .
Quantitative results for the reconstructed source integrals of the sources.The sources are listed in the order as in Table6.Their true values are 105.1 and 97.4,respectively.The unit is nW.

Table 8 .
Finite element information for the physical phantom experiment