Open Access
Add to BibSonomyAbstract
Bioluminescence tomography (BLT) is an effective molecular imaging (MI) modality. Because of the ill-posedness, the inverse problem of BLT is still open. We present a trust region method (TRM) for BLT source reconstruction. The TRM is applied in the source reconstruction procedure of BLT for the first time. The results of both numerical simulations and the experiments of cube phantom and nude mouse draw us to the conclusion that based on the adaptive finite element (AFE) framework, the TRM works in the source reconstruction procedure of BLT. To make our conclusion more reliable, we also compare the performance of the TRM and that of the famous Tikhonov regularization method after only one step of mesh refinement of the AFE framework. The conclusion is that the TRM can get faster and better results after only one mesh refinement step of AFE framework than the Tikhonov regularization method when handling large scale data. In the TRM, all the parameters are fixed, while in the Tikhonov method the regularization parameter needs to be well selected.
©2010 Optical Society of America
I. Introduction
Among many molecular imaging modalities, optical imaging, especially bioluminescence imaging, has attracted remarkable attention for its unique advantages in probing capabilities, sensitivity, specificity, and cost-effectiveness [1, 2, 3] in cancer research [4] and drug development [5]. By utilizing the surface light distribution of an object, BLT can reconstruct the bioluminescent light source distribution inside, which is called the inverse source problem of partial differential equations (PDE) [6]. The inverse problem of BLT is an ill-posed problem.
Till now, the inverse problem of BLT is still open. A common way to overcome the ill-posed property is the regularization. The Tikhonov regularization strategy is usually used [7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20]. Given a suitable regularization parameter, the Tikhonov method can get excellent results. But it is difficult to choose a proper regularization parameter. Here we want to introduce a new kind of method into the reconstruction procedure of BLT problem, the trust region method with all involved parameters fixed.
As of the origin of the TRM, the work by Levenberg [21] in 1944 are the mostly cited among the papers related to TRM. And a more direct link to TRM is the work by Marquardt in 1963 [22]. Based on the work of Levenberg and Marquardt, Powell derives the first trust region algorithm in 1970s [23, 24]. After Powell’s work, TRM becomes prosperous. The convergence and regularity of the standard trust region method when applying it to ill-posed problems has also been studied by Yanfei Wang and Yaxiang Yuan [25]. Comparing with the Tikhonov method, the advantages of TRM are that it can get better results after only one step of mesh refinement process in the AFE framework, it is faster when handling large scale data and the parameters in TRM are all fixed while in the Tikhonov method the regularization parameter has to be properly selected. So we introduce the TRM to solve the BLT inverse problem.
The paper is organized in the following sequence. In section 2, we firstly formulate the BLT forward problem. Then after a brief introduction of the AFE framework that the TRM works in for the inverse problem, we will focus on the TRM for the regularization and optimization procedure. In section 3, the results of both numerical simulations and physical experiments on cube phantom and nude mouse can draw us to the conclusion that TRM can work in solving BLT inverse problem in the AFE framework comparing with the Tikhonov regularization method. Finally, we will give our comments and conclude the paper in section 4.
2. Method
2.1. Diffusion model for forward problem
In bioluminescence imaging, biological entities, such as tumor cells, genes and compounds of drug, are tagged with luciferase enzymes. When the luciferase is combined together with the substrate luciferin, oxygen and ATP, a biochemical reaction occurs that transforms part of the chemical energy into the bioluminescent photons with a wavelength of about 600nm [26]. The radiative transfer equation (RTE) [27, 28, 29] can describe the photon propagation. However, the directly solving of RTE is very costly. In near infrared light spectrum, photon scattering predominates over absorption in the biological tissue, the photon propagation can be described by the following steady-state domain (SSD) diffusion equation (DE) depending on the light wavelength λ:
where Ω is the problem domain; Φ(x, λ) is the photon flux density [Watts/mm 2]; S(x, λ) is the bioluminescent source density [Watts/mm 3]; μa(x,λ) is the absorption coefficient [mm -1]; D(x,λ) = 1/(3(μa(x,λ)+μs ′(x,λ))) is the optical diffusion coefficient [mm];μs ′(x,λ) = (1 - g)μs(x,λ) is the reduced scattering coefficient; μs(x,λ) is the scattering coefficient [mm -1] and g is the anisotropy parameter.
When the bioluminescence imaging experiment is performed in a totally dark environment, that is to say no external photon travels into Ω through the boundary ∂Ω, we can get the boundary condition (Robin boundary condition) for DE [30, 31].
where ∂Ω is the boundary of the problem domain and v is the unit outer normal on ∂Ω. As the experiment is usually carried out with the surrounding medium Ω being air, n ′ is approximately 1. So the index mismatch parameter A(x;n,n ′) can be approximated by A(x;n,n ′) ≈ (1 + R(x))/(1 - R(x)) to incorporate diffuse boundary reflection arising from a refractive index mismatch between the problem domain Ω and the surrounding medium, where R(x) is a parameter governing the internal reflection at the boundary ∂Ω and can be approximated with R(x) ≈ 1.4399n -2 + 0.7099n -1 + 0.6681 + 0.0636n [30]. The measured quantity is the outgoing photon density on ∂Ω [31]:
2.2. Adaptive finite element framework for inverse problem
According to the finite element theory [32], in the Sobolev space H 1(Ω), we can get the weak solution of the flux density Q(x, λ) through Eqs. (1) and (2):
According to the adaptive finite element framework introduced by Lv et al. [7], we can get the matrix form of Eqs. (3) for the l-th level of the mesh refinement process: ([Kl]+[Cl]+[Bl])Φl = Ml Φl = Fl Sl where the components of the matrices Kl,Cl,Bl are obtained by
kij (l),cij (l),bij (l) and sij (l) are the elements of K (l),C (l),B (l) and S (l) with the row number i and column number j, respectively. φi (l)(x) and φj (l)(x) are the interpolation basis functions. si (l) is the source density at i. Ml is a symmetric positive-definite matrix. After applying the permissible source region method [33], the linear relationship between the boundary measured photon flux density Φlmeas and the unknown source density in the permissible source region SlP can be obtained:
where Al can be got by retaining those columns of Ml -1 Fl corresponding to SlP and those rows corresponding to Φlmeas. Then, the objective function fl(x) of the l-th level can be obtained as
The inverse problem of BLT is to reconstruct a 3D bioluminescent source distribution inside an object, such as a mouse or phantom, given the surface bioluminescence distribution of the object. Mahtematically, BLT is to reconstruct the source distribution SlP in Eq. 5 according to the measured surface distribution Φlmeas.
2.3. Trust region method
The optimization procedure is the minimization problem
where f(x), A , x and b stand for fl(SlP), Al,SlP and Φlmeas respectively in Eq. (5) for short. Since the inverse problem of BLT is ill-posed, hence regularization is necessary. The trust region method is usually posed for well-conditioned problems. The reason of applying trust region method to ill-posed inverse problems is that TRM was proved to be a regularization method [25]. Now, we will give the detailed description of TRM introduced by Wenyu Sun and Yaxiang Yuan [34]. The basic idea of TRM is to approximate the objective function f(x) around xk by choosing a quadratic model of the form
where k denotes the k-th iteration of the calculation in TRM, gk = AT Axk - ATb,Gk = AT A and use the minimizer sk of qk(s) to modify xk,
Then, we define a region around the current iterate
where Δk is the radius of Ωk, where the model is trusted to be adequate to the objective function. Then we choose a step to be the approximate minimizer of the quadratic model in the trust-region, such that xk + sk is the approximately best point on the generalized sphere
with center xk and radius Δk. If the step is not acceptable, the size of the trust-region will be reduced and we will continue to find a new minimizer. Since the step is restricted by the trust region, it is also called the restricted step method.
The model subproblem of the trust-region method can be formatted as
where Δk > 0 is the trust-region radius, Bk is symmetric and approximate to the Hessian Gk. In Eq. (6), we set Bk = Gk and the method is called a Newton-type trust-region method. In general, when there is good agreement between the model qk(s) and the objective function value f(xk + s), one should select Δk as large as possible. Let
which is called the actual reduction, and let
which is called the predicted reduction. Define the ratio
which measures the agreement between the model function q (k) and the objective function f. The ratio rk plays an important role in selecting new iterate x k+1 and updating the trust-region radius Δk. If rk is close to 1, it means there is good agreement, and we can expand the trust-region for the next iteration; if rk is close to zero or negative, we shrink the trust-region; otherwise, we do not alter the trust-region. The following is the trust-region algorithm:
- Given initial point x 0, , Δ0 ∈ (0, ), ε ≥ 0, 0 < η 1 ≤ η 2 < 1 and 0 < γ 1 < 1 < γ 2, k := 0.
- If ∥gk∥ ≤ η, stop.
- Approximately solve the subproblem Eq. (6) for sk.
- Compute f(xk + sk) and rk. Set
- If rk < η 1, then Δk+1 ∈ (0, γ 1Δk];
If rk ∈ [η 1, η 2), then Δk+1 ∈ [γ 1 Δk, Δk];
If rk ≥ η 2 and ∥sk∥ = Δk, then Δk+1 = min(γ 2Δk, );
else Δk+1 = θ∥sk∥.
- Generate B k+1, update qk, set k := k + 1, go to Step 2.
In the above algorithm, is an overall bound for all Δk. The iterations for which rk ≥ η 2 and thus for which Δk+1 ≥ Δk, are said to be very successful iterations; the iterations for which rk ≥ η 1 and thus for which x k+1 = xk + sk, are said to be successful iterations; otherwise the iterations for which rk < η 1 and thus for which x k+1 = xk, are said to be unsuccessful iterations. Sometimes, the iterations in the first two cases are said to be successful iterations.
2.4. Relevant application issues
All the AFE framework and the Tikhonov regularization method are implemented in C++ code and the TRM is implemented in Matlab. The reconstruction procedure is performed on an Intel(R) Core(TM)2 Duo E4600 CPU(2.4GHz) platform with 2GB memory.
The processing speed of AFE framework is increased by using the multi-thread ZSM technology [35] when generating Al in Eq. (4).
All the parameters of the AFE framework and Tikhonov regularization method are the same as those used in literature [7] except the refinement threshold that is responsible for controlling the mesh quality during the mesh refinement procedure. The smaller the refinement threshold is set, the finer and more dense mesh is generated. In another way, bigger refinement threshold will result in more coarse mesh. A modified Newton method is used to solve the minimization problem in the Tikhonov regularization procedure.
As of the parameters of TRM, we set = 0.25, η 1 = 0.01, η 2 = 0.75, γ 1 = 0.5, γ 1 = 2, and η = 1 × 10-50. We use
to update Δk+1 in Step 5. We use the “TRUST” function that is supplied by Matlab tool box to solve the trust region subproblem of Eq. (6).
In order to analyze the algorithms more reasonably, we define the Distance Error of the distance between the actual source position and the reconstructed source position and Relative Error of the density between the actual source position and the reconstructed source position as:
where (x, y, z) is the center coordinate of the reconstruction source with the maximum density and (x 0, y 0, z 0) is that of the actual source center, Sreconstruct and Sreal are density of reconstruction source and actual source, respectively.
3. Experiments
3.1. Numerical simulation
We designed a heterogeneous cylindrical phantom of 30mm in height and 10mm in radius. The phantom consisted of four ellipsoids and one cylinder to represent muscle, lungs, heart, bone and liver, as shown in Fig. 1(a). The phantom was dis The optical parameters were all obtained from the literature [36] and listed in Table 1.
Table 1. Optical parameters of different tissues of the heterogeneous cylindrical phantom
Since Monte Carlo (MC) method remained a gold standard for photon transportation simulation because of its accuracy and flexibility [37, 38], when carrying out numerical experiments, we used the MC method based molecular optical simulation environment (MOSE) [39] that could take 2D/3D analytical models, micro-CT and micro-MRI slices to define the object geometry to get the surface photon distribution data. Besides, MC method could avoid the inverse crime problem. In the MOSE simulation, the source was sampled by 106 photons and was assumed to obey the uniform distribution. The aforementioned heterogeneous phantom was discretized into 34072 triangles and 11499 surface measurement points with an average element diameter of about 0.5mm for MOSE simulation, while in the reconstruction procedure, the aforementioned heterogeneous phantom was discretized into 1537 points and 6878 tetrahedrons with an average element diameter of about 2mm.
Fig. 1. Heterogeneous cylindrical numerical phantom with single source (a), consisted of muscle (white), bone (black), heart (pink), lungs (green), liver (yellow) and a ball source (blue) in the right lung. Homogeneous cube phantom with single source (b) and double sources (c). Those blue cylinders in sub figures (b) and (c) denote the light source.
A solid sphere source of 1mm in radius and 0.238nano – Watts/mm 3 in power density was centered at (3,5,15) inside the right lung as shown in Fig. 1(a). To reduce the ill-posedness of the BLT inverse problem, we incorporated the permissible source region of
as a priori information, according to the surface light distribution. The refinement threshold was set to 7 × 10-3. The matrix A used in the last reconstruction procedure was 677 × 487. After one step of mesh refinement procedure of the AFE framework, we got the reconstruction results as shown in Fig. 2 and Table 2. Sub figures (a) to (d) in Fig. 2 were the views of Tikhonove method, while (e) to (l) were the views of TRM. The reconstructed power density of Tikhonov method was 0.088 with a relative error of 0.630. The reconstructed power density of TRM method was 0.271 with a relative error of 0.133. The reconstruction time of both the methods in the last mesh refinement procedure were shown in Table 2.
Table 2. Reconstruction results comparison between Tikhonov method and TRM in single source heterogeneous cylindrical numerical phantom case.
Fig. 2. Reconstruction results comparison between Tikhonov method (sub figures (a) to (d)) and TRM (sub figures (e) to (l)) in single source heterogeneous cylindrical numerical phantom case. Sub figures (a), (e) and (i) are 3D views; (b), (f) and (j) are front views; (c), (g) and (k) are side views; (d), (h) and (l) are top views. Sub figures (i) to (l) are zoom in images of sub figures (e) to (h) around the real source, respectively. The blue ball in each sub figure denotes the real source and the red tetrahedron denotes the reconstructed source with the maximum density. For concision, only the real source and the reconstructed source are displayed.
3.2. Physical experiments
3.2.1. Experimental setup
In our imaging system [8] as shown in Fig. 3, a liquid-nitrogen-cooled back-illuminated charge coupled device (CCD) camera (Princeton Instruments VersArray 1300B, Roper Scientific, Trenton, NJ) was adopted to collect the transmitted and scattered near-infrared photons on the surface of the phantom, and the CCD array could generate 1340 × 1300 pixels and 16 bits dynamic range images with 20mm × 20mm sized pixel. The dark current of the camera could be reduced largely through cooling the CCD chip down to – 110°C using liquid nitrogen for long exposures. Furthermore, quantum efficiency (QE) of the CCD camera was greater than 80% for the wavelength range between 500nm and 700nm. 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[40]. The rotation stage was designed to carry the imaging object, such as a phantom or a mouse, and rotate to different angles for the CCD camera to take different views. The translation stage was designed to control the distance between the imaging object and the CCD camera. The camera holder was designed for modulating the height of the CCD camera.
Fig. 3. Overview of our imaging system that consists of a CCD camera, a camera holder, a translation stage and a rotation stage [8].
3.2.2. Homogeneous cube phantom
In the phantom experiments, two cube resinous phantoms with 20mm in length were designed. The two phantoms were both made from polyoxymethylene, and one or two small holes of 1.25mm in radius were drilled in the phantoms to emplace the light source respectively, as shown in Fig. 1(b) and (c). Peroxide, ester compound and fluorescent dye were injected into the holes in the phantoms after thorough mix, and then the red light whose central wavelength located about 650nm was emitted due to the chemical reaction of the mixed resolutions. The aforementioned red light served as the internal light source in the physical phantom experiments. The phantoms were put in the imaging system as show in Fig. 3.
The final measured optical parameters of the phantom at the wavelength around 660nm were listed as follows: the absorption coefficient ua ≈ 0.0002mm -1 and the reduced scattering coefficient us ′ = (1 - g)us ≈ 1.0693mm -1. All the physical phantom experiments were performed in a light-tight imaging chamber to avoid external disturbance and light leaking. Under the computer control, the motorized rotation stage was used to rotate the phantom for acquisition of the photon flux density on the four sides of the phantom. The photo data was then mapped to the surface of the phantom for reconstruction. The aforementioned homogeneous cube phantom was discretized into 1457 points and 6444 tetrahedrons with an average element diameter of about 2mm for the reconstruction procedure.
3.2.3. Homogeneous cube phantom experiment with single source
A cylindrical source of 2.5mm in height and 1.25mm in radius was centered at (12.5,12.5,11.25) as shown in Fig. 1(b). According to the surface light distribution, the permissible source region was set to
The refinement threshold was set to 1 × 10-3. The matrix A used in the last reconstruction procedure was 681 × 3925. After one step of mesh refinement procedure of the AFE framework, we got the reconstruction results as shown in Fig. 4 and Table 3. Sub figures (a) to (d) in Fig. 4 were the views of Tikhonov method, while (e) to (l) were the views of TRM. The reconstruction time of both the methods in the last mesh refinement procedure were shown in Table 3.
Fig. 4. Reconstruction results comparison between Tikhonov method (sub figures (a) to (d)) and TRM (sub figures (e) to (l)) in single source homogeneous cube phantom case. Sub figures (a), (e) and (i) are 3D views; (b), (f) and (j) are front views; (c), (g) and (k) are side views; (d), (h) and (l) are top views. Sub figures (i) to (l) are zoom in images of sub figures (e) to (h) around the real source, respectively. The blue cylinder in each sub figure denotes the real source and the red tetrahedron denotes the reconstructed source with the maximum density.
Table 3. Reconstruction results comparison between Tikhonov method and TRM in single source homogeneous cube phantom case
3.2.4. Homogeneous cube phantom experiment with double sources
Two cylindrical sources of 2.5mm in height and 1.25mm in radius were centered at (6.25,6.25,11.25) and (13.75,13.75,11.25) respectively, as shown in Fig. 1(c). According to the surface light distribution, the permissible source region was set to
The refinement threshold was set to 6.5 × 10-2. The matrix A used in the last reconstruction procedure was 680 × 3082. After one step of mesh refinement procedure of the AFE framework, we got the reconstruction results as shown in Fig. 5 and Table 4. Sub figures (a) to (d) in Fig. 5 were the views of Tikhonov method, while (e) to (l) were the views of TRM. The reconstruction time of both the methods in the last mesh refinement procedure were shown in Table 4.
Fig. 5. Reconstruction results comparison between Tikhonov method (sub figures (a) to (d)) and TRM (sub figures (e) to (l)) in double sources homogeneous cube phantom case. Sub figures (a), (e) and (i) are 3D views; (b), (f) and (j) are front views; (c), (g) and (k) are side views; (d), (h) and (l) are top views. Sub figures (i) to (l) are zoom in images of sub figures (e) to (h) around the real source, respectively. The blue cylinder in each sub figure denotes the real source and the red tetrahedron denotes the reconstructed source with the maximum density.
Table 4. Reconstruction results comparison between Tikhonov method and TRM in double sources homogeneous cube phantom case
3.2.5. Mouse experiment
Besides the numerical simulation experiment and the real phantom experiment, we’d successfully carried out an experiment on a nude mouse. For getting better anatomical structure information in Micro-CT scanning, 0.3ml Fenestra VC was injected intravenously into the lateral tail vein of the nude mouse. A cylindrical light source with 3mm in length and 2mm in diameter that was made of hollow plastic catheter filled with lightening material was sewed into the epigastrix torso of the mouse. The mouse was put in a mouse holder that could keep the mouse from moving during the data acquisition procedure. The mouse holder including the mouse was then put in the rotation stage in Fig. 3. The rotation stage was then set to rotate to 0°, 90°, 180° and 270° for taking photos at those positions. At each position, one photo of white light and one photo of bioluminescence light were taken. After acquiring the optical data, the mouse was scanned by the Micro-CT system to incorporate the multi-modality information [41] and the platform independent parameters of power, voltage, and exposure time were set to 50W, 50Kvp, and 0.467s, respectively. The number of views was 360 with 1120 × 2344 pixel per view. The raw CT data was then reconstructed and the center of the real source was (25.0,20.0,7.8). The dimension of the reconstructed CT data was 512 × 512 × 720. The voxel of reconstructed CT data was 0.1mm × 0.1mm × 0.1mm. The reconstructed CT data was segmented into 5 major tissues, including muscle, lungs, heart, liver and bone. The segmented CT data was discretized into 4361 points and 22614 tetrahedrons for the reconstruction procedure as shown in Fig. 6(a). Different from the mesh used in the numerical and cube phantom experiments, the mesh used in the nude mouse experiment had to be generated from segmented CT slices, as those mesh used in numerical and cube phantom experiments had regular shape, which made it possible to be generated on computer, given the geometric information of the object. The 2D photos were then mapped onto the surface of the mesh as shown in Fig. 6(b). The optical parameters of the mouse were shown in Table 5. According to the surface light distribution, the permissible source region was set to
Fig. 6. Sub figure (a) is the mesh used in the reconstruction procedure. The mesh consists 5 tissues, including the heart in blue, the bone in red the lung in yellow, the liver in green and the muscle in gray. Sub figure (b) is the 3D bioluminescence mapping result from 2D bioluminescence photos.
Table 5. Optical parameters of the nude mouse
The refinement threshold was set to 1 × 10-3. The matrix A used in the last reconstruction procedure was 1457 × 4235. After one step of mesh refinement procedure of the AFE framework, we got the reconstruction results as shown in Table 6 and Fig. 7. Sub figures (a) to (d) in Fig. 7 were the views of Tikhonov method, while (e) to (l) were the views of TRM. The reconstruction time of both the methods in the last mesh refinement procedure were shown in Table 6.
Fig. 7. Reconstruction results comparison between Tikhonov method (sub figures (a) to (d)) and TRM (sub figures (e) to (l)) in single source heterogeneous nude mouse case. Sub figures (a), (e) and (i) are 3D views; (b), (f) and (j) are front views; (c), (g) and (k) are side views; (d), (h) and (l) are top views. Sub figures (i) to (l) are zoom in images of sub figures (e) to (h) around the real source, respectively. The blue cylinder in each sub figure denotes the real source and the red tetrahedron denotes the reconstructed source with the maximum density. For concision, only the real source and the reconstructed source are displayed.
Table 6. Reconstruction results comparison between Tikhonov method and TRM in single source heterogeneous nude mouse case
4. Discussions and conclusions
Based on the adaptive finite element framework, the trust region method has been proposed the first time for BLT. In order to compare with the famous Tikhonov regularization method, we’ve carried out experiments both numerically and physically. The results could give us the conclusion that TRM can work in solving BLT inverse problem.
Compared with the Tikhonov method, the TRM has the following advantages:
Firstly, the TRM can get better results after only one step of mesh refinement procedure in the AFE framework by modifying the refinement threshold. In order to save more processing time, we only perform one mesh refinement step in the AFE framework while usually more steps are taken. Theoretically, the more the refinement steps are carried out, the better the reconstruction results. However, the more processing time will surely be cost if more mesh refinement steps are taken. In this way, we can take the advantages of the AFE framework and at the same time save processing time.
Secondly, the TRM is more time efficient than the Tikhonov method when dealing with large scale data. In the numerical experiment above, the permissible source region is restricted in the right lung in order to compare with previous work. As a result, the dimensions of matrix A used in the reconstruction procedure are relatively small. So the processing time of the TRM is longer than that of the Tikhonov method. In physical cube phantom and nude mouse experiment cases, the permissible source regions are all big compared with the one used in numerical experiment. And the processing time of the TRM in the physical experiment cases is shorter than that of the Tikhonov method when the dimensions of matrix A becomes large.
Thirdly, the TRM can eliminate the need of choosing regularization parameter as all the used parameters are fixed in TRM while in the Tikhonov method the regularization parameter has an important influence. And it is difficult to choose an appropriate regularization parameter.
We also discover that by modifying the refinement threshold of the AFE frame work, we can get acceptable reconstruction results only after one mesh refinement procedure to save more processing time. The influence of the refinement threshold is beyond the scope of this paper and will be reported later.
The numerical experiment is designed to compare with previous work. The physical cube phantom and real nude mouse experiments are designed to demonstrate that TRM can handle large scale data acquired in real experiments. In the future, our work will be focused on the mouse experiments of tumor model to show the biological application of the proposed TRM.
Acknowledgments
This paper is supported by the Project for the National Basic Research Program of China (973) under Grant No.2006CB705700, Changjiang Scholars and Innovative Research Team in University (PCSIRT) under Grant No.IRT0645, CAS Hundred Talents Program, CAS scientific research equipment develop program under Grant No. YZ200766, the Knowledge Innovation Project of the Chinese Academy of Sciences under Grant No. KGCX2-YW-129, KSCX2-YW-R-262, the National Natural Science Foundation of China under Grant No. 30672690, 30600151, 60532050, 60621001, 30873462, 60910006, 30970769, 30970771, Beijing Natural Science Fund under Grant No.4071003, Science and Technology Key Project of Beijing Municipal Education Commission under Grant No.KZ200910005005.
References and links
1. R. Weissleder and V. Ntziachristos, “Shedding light onto live molecular targets,” Nature Medicine 9, 123–128 (2003). [CrossRef] [PubMed]
2. V. Ntziachristos, J. Ripoll, L. H. V. Wang, and R. Weissleder, “Looking and listening to light: the evolution of whole-body photonic imaging,” Nature Biotechnol. 23, 313–320 (2005). [CrossRef]
3. M. K. So, C. J. Xu, A. M. Loening, S. S. Gambhir, and J. H. Rao, “Self-illuminating quantum dot conjugates for in vivo imaging,” Nature Biotechnol. 24, 339–343 (2006). [CrossRef]
4. R. Weissleder and M. J. Pittet, “Imaging in the era of molecular oncology,” Nature 452, 580–589 (2008). [CrossRef] [PubMed]
5. J. K. Willmann, N. van Bruggen, L. M. Dinkelborg, and S. S. Gambhir, “Molecular imaging in drug development,” Nat. Rev. Drug Discov. 7, 591–607 (2008). [CrossRef] [PubMed]
6. V. Isakov, Inverse Problems for Partial Differential Equations (Springer-Verlag, New York, 1998).
7. 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. [CrossRef] [PubMed]
8. C. Qin, J. Tian, X. Yang, J. Feng, K. Liu, J. Liu, G. Yan, S. Zhu, and M. Xu, “Adaptive improved element free Galerkin method for quasi or multi spectral bioluminescence tomography,” Opt. Express 17, 21925–21934 (2009), http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-17-24-21925. [CrossRef] [PubMed]
9. J. Feng, K. Jia, C. Qin, G. Yan, S. Zhu, X. Zhang, J. Liu, and J. Tian, “Three-dimensional Bioluminescence Tomography based on Bayesian Approach,” Opt. Express 17, 16834–16848 (2009), http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-17-19-16834. [CrossRef] [PubMed]
10. Y. Lu, X. Zhang, A. Douraghy, D. Stout, J. Tian, T. F. Chan, and A. F. Chatziioannou, “Source Reconstruction for spectrally-resolved bioluminescence tomography with sparse a priori information,” Opt. Express 17, 8062–8080 (2009), http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-17-10-8062. [CrossRef] [PubMed]
11. Y. Lu, A. Douraghy, H. B. Machado, D. Stout, J. Tian, H. Herschman, and A. F. Chatziioannou, “Spectrally-resolved bioluminescence tomography with the third-order simplified spherical harmonics approximation. Physics in Medicine and Biology,” 59, 6477–6493 (2009).
12. 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. [CrossRef] [PubMed]
13. Y. Lv, J. Tian, H. Li, 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, 1–16 (2007). [CrossRef]
14. W. M. Han, W. X. Cong, and G. Wang, “Mathematical theory and numerical analysis of bioluminescence tomography,” Inverse Problems 22, 1659–1675 (2006). [CrossRef]
15. V. Soloviev, “Tomographic bioluminescence imaging with varying boundary conditions,” Applied Optics 46, 2778–2784 (2007). [CrossRef] [PubMed]
16. H. Dehghani, S. C. Davis, and B. W. Pogue, “Spectrally resolved bioluminescence tomography using the reciprocity approach,” Medical Physics 35, 4863–4871 (2008). [CrossRef] [PubMed]
17. W. Gong, R. Li, N. N. Yan, and W.B. Zhao, “An improved error analysis for finite element approximation of bioluminescence tomography,” Journal of Computational Mathematics 26, 297–309 (2008).
18. M. B. Unlu and G. Gulsen, “Effects of the time dependence of a bioluminescent source on the tomographic reconstruction,” Appl. Opt. 47, 799–806 (2008). [CrossRef]
19. X. L. Cheng, R. F. Gong, and W. M. Han, “Numerical approximation of bioluminescence tomography based on a new formulation,” Journal of Engineering Mathematics 63, 121–133 (2009). [CrossRef]
20. M. Chua and H. Dehghani, “Image reconstruction in diffuse optical tomography based on simplified spherical harmonics approximation,” Opt. Express 17, 24208–24223, (2009), http://www.opticsinfobase.org/abstract.cfm?URI=oe-17-26-24208. [CrossRef]
21. K. Levenberg, “A method for the solution of certain nonlinear problems,” Quart. Appl. Math. 2, 164–168 (1944).
22. D. W. Marquardt, “An algorithm for least-squares estimation of nonlinear parameters,” SIAM J. Appl. Math. 11, 431–441 (1963). [CrossRef]
23. M. J. D. Powell, “A new algorithm for unconstrained optimization,” in Nonlinear Programming, J. B. Rosen, O. L. Mangasarian, and K. Ritter, eds. (Academic Press, New York, 1970), 31–65.
24. M. J. D. Powell, “Convergence properties of a class of minimization algorithms,” in Nonlinear Programming, O. L. Mangasarian, R. R. Meyer, and S. M. Robinson, eds. (Academic Press, New York, 1975), 1–27.
25. Y. Wang and Y. Yuan, “Convergence and regularity of trust region methods for nonlinear ill-posed inverse problems,” Inverse Problems 21, 821–838, (2005). [CrossRef]
26. C. Contag and M. H. Bachmann, “Advances in Bioluminescence imaging of gene expression,” Annu. Rev. Biomed. Eng. 4, 235–260 (2002). [CrossRef] [PubMed]
27. V. Ntziachristos, C. Tung, C. Bremer, and R. Weissleder, “Fluorescence molecular tomography resolves protease activity in vivo,” Nat. Med. 8, 757–760 (2002). [CrossRef] [PubMed]
28. R. Schultz, J. Ripoll, and V. Ntziachristos, “Experimental fluorescence tomography of tissues with noncontact measurements,” IEEE Trans. Med. Imag. 23, 492–500 (2004). [CrossRef]
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). [CrossRef] [PubMed]
30. 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). [CrossRef] [PubMed]
31. J. J. Duderstadt and L. J. Hamilton, Nuclear Reactor analysis (Wiley, New York, 1976).
32. S. S. Rao, The finite element method in engineering (Butterworth-Heinemann, Boston, 1999).
33. W. Cong, D. Kumar, Y. Liu, A. Cong, and G. Wang, “A practical method to determine the light source distribution in bioluminescent imaging,” Proc. SPIE 5535, 679–686 (2004). [CrossRef]
34. W. Sun and Y.x. Yuan, “Chapter 6 Trust-Region Methods and Conic Model Methods” in Optimization Theory and Methods: Nonlinear Programming (Springer US, 2006).
35. B. Zhang, J. Tian, D. Liu, L. Sun, X. Yang, and D. Han, “A multithread based new sparse matrix method in bioluminescence tomography”, presented at Conference 7626 of SPIE on Medical Imaging, San Diego, USA , 13–18 February 2010.
36. G. Alexandrakis, F. R. Rannou, and A. F. Chatziioannou, “Tomographic bioluminescence imaging by use of a combined optical-PET (OPET) system: a computer simulation feasibility study,” Phys. Med. Biol. 50, 4225–4241 (2005). [CrossRef] [PubMed]
37. L. H. Wang, S. L. Jacques, and L. Q. Zheng, “MCML-Monte Carlo modeling of photon transport in multi-layered tissues,” Comput. Meth. Prog. Biomed. 47, 131–146 (1995). [CrossRef]
38. D. Boas, J. Culver, J. Stott, and A. Dunn, “Three dimensional Monte Carlo code for photon migration through complex heterogeneous media including the adult human head,” Opt. Express 10, 159–169 (2002), http://www.opticsinfobase.org/abstract.cfm?URI=OPEX-10-3-159. [PubMed]
39. H. Li, J. Tian, F. Zhu, W. Cong, L. V. Wang, E. A. Hoffman, and G. Wang, “A mouse optical simulation enviroment (MOSE) to investigate bioluminescent phenomena in the living mouse with the Monte Carlo Method,” Acad. Radiol. 11, 1029–1038 (2004). [CrossRef] [PubMed]
40. 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). [CrossRef]
41. J. Tian, J. Bai, X.-P. Yan, S. Bao, Y. Li, W. Liang, and X. Yang, “Multimodality molecular imaging,” IEEE Eng. Med. Bio. Mag. 27, 48–57 (2008). [CrossRef]
