Open Access
Add to BibSonomyAbstract
Bioluminescence tomography (BLT) is a promising optical molecular imaging technique on the frontier of biomedical optics. In this paper, a generalized hybrid algorithm has been proposed based on the graph cuts algorithm and gradient-based algorithms. The graph cuts algorithm is adopted to estimate a reliable source support without prior knowledge, and different gradient-based algorithms are sequentially used to acquire an accurate and fine source distribution according to the reconstruction status. Furthermore, multilevel meshes for the internal sources are used to speed up the computation and improve the accuracy of reconstruction. Numerical simulations have been performed to validate this proposed algorithm and demonstrate its high performance in the multi-source situation even if the detection noises, optical property errors and phantom structure errors are involved in the forward imaging.
© 2013 Optical Society of America
1. Introduction
Bioluminescence tomography (BLT) is an established optical molecular imaging technique for quantitatively monitoring the in vivo biological process of small animals [1–3]. It involves a typical inverse source problem and aims to localize the bioluminescent sources inside the tissues and quantify their distributions from light detection on the surface [4]. Compared to other tomographic imaging techniques, such as X-ray computed tomography (CT), magnetic resonance imaging (MRI), positron emission tomography (PET) and single photon emission computed tomography (SPECT), etc., BLT is a nonradiative imaging technique and has low background noise, high sensitivity to the lower bioluminescent light over the visible spectral band [1, 3]. That guarantees it being competent for many different tasks of biomedical research, including gene therapy, metastasis detection, cancer diagnosis and drug discovery and development [1, 2, 5, 6].
In 2003, Wang et al. first introduced a bioluminescent tomography assembly into a multimodality tomography system to depict the internal bioluminescent source distribution upon the structures of object [4]. Hereafter, hundreds of research literatures have been reported focusing on the system development [4, 7–9], experiment scheme [10, 11] and reconstruction method [12–19]. For BLT reconstruction, the major challenge is to overcome its ill-posedness which is derived from the non-uniqueness of source distribution by nature [20]. So far, dozens of algorithms have been proposed based on optimization with different regularization terms [16, 17, 19], source support constraint [13, 15, 21–24], multi-spectral image acquisition [7, 14, 16, 17, 25], adaptive discretization [22, 26], and those mixed.
Generally, reconstruction algorithms with respect to BLT problem are mainly based on the iterative-type. In analytic case, Jiang et al. successfully employed the expectation maximization (EM) algorithm and Landweber algorithm both in complete and partial measurement situations [12, 15]. Alternatively in discrete case, Cong et al. first established a linear relationship between the internal source distribution and the surface measurement, and adopted a modified Newton algorithm with the active set strategy to solve the problem [13]. Then, Lv et al. introduced a multilevel adaptive finite element method into the reconstruction algorithm [26]. Some other gradient based algorithms, such as preconditioned conjugate gradient (PCG) algorithm [14, 16], gradient projection (GP) algorithm [16], quasi-Newton algorithm [17], trust region algorithm [18], etc., have also been applied to the BLT reconstruction. Here, it should be noted that, in order to improve the feasibility of above algorithms, a prior source support must be used.
However, a reliable source support is sometimes hardly available in practice. Lv, Feng and Naser et al. came up with different strategies for posterior estimation of the permissible source support, successively [21–23]. In these algorithms, the permissible source support in next iteration could be determined from the reconstructed source distribution in current iteration basically by regional erosion of source intensity. Nevertheless, the magnitude of reconstructed intensity at each point does not absolutely consist with its actual contribution to the source distribution. This inconsistency would probably make the source support inappropriate and cause the iteration dropping into the ambiguity.
Recently, Liu et al. developed a new BLT localization algorithm [19]. Based on the global optimization, the regularized BLT objective function can be represented by a directed graph and then be resolved by graph cuts algorithm in the whole region without source support constraint. Their single-source BLT experiment has demonstrated that the localization of reconstructed source is accurate at sub-millimeter level. However, it was shown that when we used this algorithm in the multi-source situation, the reconstructed source regions would be interconnected where sources were close to each other and this may ultimately invalidate the source localization.
In this paper, considering the complicated source situation we propose a generalized reconstruction method under the hybrid algorithm scheme. Therein, graph cuts algorithm is only adopted to estimate the source support and different gradient-based algorithms are used to reconstruct the source distribution sequentially with different performances. In this way, a more reliable source support is acquired without prior knowledge, and a more accurate and finer source distribution is finally obtained. For the further enhancement of reconstruction performance, we keep descending the source scale according to the iteration status. Numerical simulations have been performed to validate this proposed algorithm and demonstrate its high performance in the multi-source situation even if the detection noises, optical property errors and phantom structure errors are involved in the forward imaging.
This paper is organized as follows. We formulate the BLT problem based on the diffusion approximation of the radiative transfer equation in §2, and analyze the performances of some popular previously proposed algorithms in §3. In §4, a generalized hybrid algorithm is proposed and its numerical results are given in §5. Finally we conclude the paper in §6.
2. Formulation of BLT
The light propagation in a random medium can be accurately described by the radiative transfer equation (RTE) [27, 28]. However, it is quite complex and computionally expensive to reconstruct the internal source from the boundary measurements with RTE. In practice, light scattering predominates over absorption in the biological tissues. Therefore the diffusion approximation (DA) is widely used to simplify RTE in medical imaging field [27, 28].
Let Ω ∈ R3 be a bounded domain with its boundary denoted by Γ, u0(x) be the radiance at x ∈ Ω. BLT problem can be stated as follows: Given the measured outgoing radiance g(x) on Γ, find a source distribution q0(x) with the corresponding diffusion approximation u0(x) such that[15, 20]
Where μa(x) and μs(x) are absorption and scattering coefficients, respectively. D(x) = 1/[3(μa(x) + μ′s(x))] is diffusion coefficient with μ′s = (1 − η̄)μs being the reduce scattering coefficients and η̄(0 ≤ η̄ ≤ 1)) the anisotropy parameter.According to Eq. (1), there is a direct linear relationship between the measured outgoing radiance and the unknown source distribution. Assume that there are N source notes in the domain Ω, and M detectors on the boundary Γ. Then, the discrete model for BLT can be established as [13]
Where ℳ = [m1, m2, ···, mN] ∈ RM×N denotes the system matrix with mi being the ith column of it, which can be calculated based on the nodal basis functions ϕi(x) for the source q0(x) [25, 29]. q = (q1, q2, ···, qN)T ∈ RN is the discrete distribution of the source with qi = q0(xi) being the nodal value of q0(x), and b ∈ RM being the measurement on the boundary.Regularization is crucial for BLT due to its ill-posed nature. Generally, the BLT problem is converted into a regularized least-square problem
Where λ is the regularization parameter, G(q) is the regularization term. In some of BLT algorithms, the source support constraint has been well used as a substitute for the regularization term.3. Reconstruction of BLT
In this section we briefly introduce some previously proposed reconstruction algorithms for BLT which will be referred to in our proposed algorithm.
3.1. EM, Landweber and modified Newton algorithm
EM, Landweber and modified Newton are three popular iterative algorithms in BLT reconstruction. Here, we only give their iterative schemes and corresponding references.
Where 1 ∈ RM is a vector with only 1’s components.Landweber algorithm[12, 15, 30–32]
Modified Newton algorithm[13, 26, 33]
3.2. Graph cuts algorithm
Recently, a graph cuts (GC) algorithm has been presented to localize the bioluminescent source in BLT reconstruction which has been widely used in computer vision [19, 34, 35].
In the graph framework, the energy function Eq. (3) can be reformulated as follows
Where, θconst is the constant term of the energy, θi(·) are the unary terms which denote the weight of the edges between the source nodes and the terminals, and θij(·, ·) are the pairwise terms which denote the weight of the edges between the internal source nodes. The concrete forms of θconst, θi and θij depend on the regularization term G(q). When , i.e. the L2 regularization is chosen, the forms of θconst, θi and θij are given in [19, Eqs. 8(a)–8(c)]. When G(q) = ||q||L1, i.e. the L1 regularization is chosen, θconst, θi and θij are as follows Other kinds of regularizations, such as total variation (TV) regularization, can also be chosen for BLT.Based on the quadratic pseudo-boolean optimization (QPBO), E(q) can be coverted to a graph-representable function and a directed graph can be constructed for it [19, 35, 36]. Then a graph cuts algorithm can be adopted to minimize the objective function and get the source distribution in polynomial time. In this paper we choose the min-cut/max-flow algorithm proposed by Boykov and Kolmogorov to compute the min-cut whose efficiency has been validated in computer vision [34].
3.3. Performance of above algorithms
To fully investigate the performances of EM, Landweber and modified Newton algorithms, numerical simulations have been carried out whilean appropriate source support is given. The results of those algorithms presented similar source distribution estimations but different speed of convergence. Therein, EM algorithm converges slightly faster than Landweber algorithm, and modified Newton converges remarkably faster than both EM and Landweber. However, the modified Newton algorithm needs to compute the Hessian inverse which is very time consuming for large-scale problems.
Numerical simulations have also been carried out for investigating the performance of graph cuts algorithm. The results demonstrated that when the bioluminescent source is single or the sources are of enough distance from each other, graph cuts algorithm could accurately locate the source (sources). Otherwise, if the sources are close to each other, the reconstructed regions would be interconnected and ultimately invalidate the source location. Unexpectedly, we found that the reconstructed results is an ideal estimation of the source support which can be used for other BLT algorithms.
4. Generalized hybrid algorithm
In this section, we establish an optimized solution for the BLT problem based on the hybrid algorithm scheme. As introduced in former context, most of BLT reconstruction algorithms are gradient-based and could be mainly classified into two groups according to the order of derivative those algorithms concerning with. For our statement convenience, one group is entitled first-order algorithm, and the other group is the second-order algorithm.
4.1. Algorithm description
As shown in Fig. 1, this hybrid algorithm involves two steps, including source support computation and source distribution computation. Each step has its own generalized iteration scheme to optimize the objective estimation. Support computation mainly depends on the graph cuts algorithm to achieve a more reliable and finer support. And source distribution computation considerately combines the first-order algorithms with second-order algorithms to accelerate the calculation.
It is known that, using graph cuts algorithm to estimate a binary source distribution in BLT problem is namely to label source nodes essentially independent on the magnitude of source intensity. Meanwhile, to guarantee the validity of the reconstruction via graph cuts algorithm, one first-order algorithm with low computational cost needs to be pre-executed to provide the source number and total power for it. And being well adjusting the iteration factors, the iteration scheme involving these two algorithms can be built up in source support computation step (indicated as the right dashed box in Fig. 1) to make the support shrink to some stable regions. Then to the step of source distribution computation, all the gradient-based algorithms are appropriate as the source support is given. Considering the speed of algorithm convergence, one second-order algorithm is first executed to obtain the primary source distribution, and one first-order algorithm is executed to update the source distribution successively. Furthermore, based on the multilevel meshes concept, we also build up an iteration scheme here (indicated as the left dashed box in Fig. 1) to descend the scale of source discretization step by step to improve the performance of BLT reconstruction.
The generalized hybrid algorithm is briefly described as follows. Here, we assume there is a Cartesian mesh level sequence L1, L2,... , LK0 for the internal sources, with K0 the maximum number of mesh refinement.
Algorithm 1. (Generalized hybrid algorithm for BLT)
4.2. Relevant issues
4.2.1. Estimation of source number and source power
In the support computation step, after using a first-order algorithm to preliminarily obtain a source distribution, the source total power could be easily estimated by integrating the source distribution. Then using the region growing algorithm, the source number could also be estimated [37]. At first the local extreme points of the source distribution were taken as the seed points, and after region growing we took the number of disjoint regions as the source number.
4.2.2. Source deblurring
In the distribution computation step, due to the multilevel meshes iteration scheme, the source may be blurred due to the distribution prolongation from the coarser mesh to the finer mesh. Therefore, an efficient image deblurring algorithm needs to be adopted here. From the experience of our numerical simulations, the Wiener filtration and a pan-like point spread function is suggested [37].
4.2.3. Related criteria in iterations
In the support computation step, the criterion aims to indicate whether the estimations of source number and total power are getting stable. In distribution computation step, the criteria involve the maximum number of source scale descending and residual error of the distribution in current iteration which should ultimately decrease to a neglectable level.
5. Numerical simulations
Numerical simulations have been performed to validate the proposed algorithm and three influence factors in the forward imaging model are taken into account. For better evaluating the reconstructed results quantitatively, the location and power of each source are respectively calculated which are mostly independent on the non-uniqueness of the source distribution.
5.1. Simulation settings
5.1.1. Phantom and optical properties
In our numerical simulations, we use the same phantom and optical property parameters as reported in [13, 15]. The cylindrical heterogeneous mouse chest phantom of 30mm height and 30mm diameter contains four kinds of materials to represent bone (B), heart (H), lungs (L), and muscle (M). The structure of the phantom is shown in Fig. 2(a) and Fig. 2(b), and the optical properties of the four components are given in Table 1.
Table 1. Optical Properties for the Numerical Phantom
5.1.2. Source setting
Four spherical sources are set up in the simulations, which are shown in Fig. 2(c). The detailed information about the source radius, intensity and power are shown in Table 2.
Table 2. Setting up of the Internal Sources
5.1.3. Algorithm setting
In our hybrid algorithm, L1 regularization is employed in the graph cuts algorithm. In fact, numerical simulations have demonstrated that the differences among the reconstructed results respectively based on L1, L2 and TV regularization are in a negligible level. That because the graph cuts algorithm is only adopted in the support computation step, and the differences among these regularizations are dramatically suppressed by the coarse mesh used in that step.
In the numerical simulations, we select EM algorithm as the first-order algorithm to compute source distribution and provide source number and total power for graph cuts algorithm in the support computation step of the hybrid algorithm. In the source distribution computation step, the modified Newton algorithm is first applied to obtain the primary source distribution, and then the Landweber algorithm is applied to update the result.
Numerical simulations indicate that the estimations of source number and source power are getting stable after two iterations, and two levels of meshes are enough for the algorithm to acquire a fairly good source distribution. Hence, we terminate the source support computation step after two iterations and the maximum number of mesh refinement K0 is set to 2 manually in our simulations.
5.1.4. Source discretization and forward process
The nodal basis functions are chosen as piecewise constants for internal source notes in our numerical simulations. That is, for source node xi, the nodal basis function is
Where Δxi is the spatial resolution. There are two levels of meshes for the internal sources in our hybrid algorithm, for the coarse mesh Δxi = 3mm, and for the fine mesh Δxi = 1mm. The system matrices corresponding to the meshes are computed in advance and stored before source reconstruction procedure.Both the external measurements and the system matrices are obtained based on finite element method (FEM). In this work, all the algorithms are implemented using Matlab™, and performed in a desktop computer with Intel Core 2 Duo 2.93 GHz CPU and 4GB RAM. To avoid the inverse crime, we use different finite elements for the computation of the external measurements and the system matrices, respectively. The information of the finite elements is given in Table 3.
Table 3. Finite Element Information for the Simulations
Furthermore, for the forward process of BLT we have carried out numerical simulations to compare the external measurements obtained by Monte Carlo (MC) method and finite element method (FEM) respectively. To guarantee the accuracy of the Monte Carlo simulations, 108 rays were traced for each source. The results show that the measurements obtained by FEM are in good agreement with that obtained by MC method, and compared with MC, the FEM had an excellent computation efficiency. The results also demonstrate the diffusion approximation equation is a proper approximation to RTE in weakly absorbing and highly scattering media, and FEM is capable of solving the forward problem of BLT efficiently with a high accuracy in this case.
5.2. Algorithm comparison
In this part, we compare our hybrid algorithm with EM, Landweber, modified Newton algorithm and the algorithm that proposed in [22] in the ideal case. The ideal case means that the phantom structure and the optical properties are accurately known, and the boundary measurements are free of noise.
Firstly, comparison between the hybrid algorithm and EM, Landweber and modified Newton algorithm have been carried out according to the simulation settings in §5.1. And the prior source support used for those single algorithms is set as
Figure 3 illustrates the reconstructed results obtained by the hybrid algorithm, EM, Landweber and modified Newton algorithm, respectively. As shown in Figs. 3(a) and 3(b), the proposed algorithm is capable of reconstructing the source distribution and separating the sources close to each other. In contrast, as shown in Figs. 3(c)–3(h), although the EM, Landweber and modified Newton algorithm can also recover the source distribution, they are incapable of well separating the sources that close to each other.
Fig. 3 Reconstructed results of the hybrid algorithm, EM, Landweber and modified Newton algorithm. (a), (c), (e) and (g) are the cross-section at z = 0mm of the reconstructed results corresponding to the four algorithms, and the white circles with actual source radius (1mm) are marked to show the true source location. (b), (d), (f) and (h) are the 3D view of the reconstructed results.
Table 4 shows the quantitative comparison of the reconstructed results between the hybrid algorithm, EM, Landweber and modified Newton algorithm. In Table 4 LE (Location Error) denotes the error of the reconstructed source center, the components of LE denotes the errors in x, y and z directions respectively; RMS (Root Mean Square) denotes the root mean square of LE, i.e. the distance between the centers of the actual and reconstructed sources. PE (Power Error) denotes the error of the reconstructed source power; RPE (Relative Power Error) denotes the relative error of the reconstructed source power; Source Size denotes the radius of reconstructed source; Computation Time denotes the time cost of the reconstruction except the time of computing the system matrix.
Table 4. Quantitative Comparison of the Reconstructed Results Between the Hybrid Algorithm and the Single EM, Landweber, and Modified Newton Algorithm
The quantitative analysis indicates that, for the hybrid algorithm, the location error from the actual source is within 0.02mm, and the relative error of the reconstructed source power is within 0.3%. In contrast, for the EM algorithm, the location error is greater than 0.05mm, and the power error is greater than 18%; for the Landweber algorithm, the location error is greater than 0.31mm, and the power error is greater than 22%; for the modified Newton algorithm, the location error is greater than 0.11mm, and the power error is greater than 14%.
As it shown in Table 4, all the algorithms overestimated the source size more or less. However, the reconstructed source radius by our hybrid algorithm is much closer to the true source radius (1mm) compared to the other algorithms. It is found that the maximum intensity of each reconstructed source obtained by the hybrid algorithm is greater than the true value, that because the power of the reconstructed source mainly concentrates in a region which is smaller than the actual source support.
It can be seen that the proposed hybrid algorithm can avoid the prior knowledge of the source support and obtain more accurate reconstructed results than the single EM, Landweber and modified Newton algorithm. In fact, the hybrid algorithm is also more efficient due to the use of multilevel meshes for the internal sources, the whole reconstruction process of the hybrid algorithm costs less than 20s.
Secondly, we compare our algorithm with the algorithm proposed in literature [22] which can self-provide the source support during iteration. Here, we implement the same phantom, optical property parameters and source setting as that in [22] to promote the comparison. And the difference is we only used the monochromatic measured data in [500nm, 550nm] rather than multispectral measured data in the reconstruction. The reconstructed source center obtained by our algorithm is (−3.01, 5.05, 15.00) with a 0.05mm location error to the actual center (−3, 5, 15), which is much smaller than the 1.36mm location error that reported in [22]. Moreover, the reconstructed power obtained by our algorithm is 0.9844nw, and the power error is only 1.56%. Quantitative comparison shows that our algorithm is superior to the algorithm proposed in [22] and has a better performance though only a single band measurement is used.
5.3. Consideration of influence factors involved
Influence factors taking into account here are the detection noise, optical property error, phantom structure error and those mixed.
5.3.1. Detection noise
The external measurements are usually captured by the high sensitive CCD camera in practice. To explore the influence of the detection noise on the hybrid algorithm, different levels of Gaussian noises are added to the measured data. In the simulations, we discretized the measured data to a gray image with a maximum value 104, and add a zero-mean Gaussian noise with standard deviation σ, Gaussian(0, σ), to the image.
Simulations have been carried out with different levels of Gaussian noise (σ = 10, 20, 50, 100), respectively. The numerical results demonstrate that the reconstructed results are slightly influenced by Gaussian noise even if the standard deviation is up to 100. The location and power errors of the reconstructed source slightly increase with the increasing of the Gaussian noise. Table 5 shows the quantitative results obtained by the hybrid algorithm with the measurement being corrupted by Gaussian(0,100). In this case the location error is within 0.07mm, and the power error is within 3% which is only slightly different from the quantitative results in the ideal case. The reconstructed result is not given here as it is very similar to the result illustrated in Figs. 3(a) and 3(b).
Table 5. Quantitative Analysis of the Reconstructed Results Obtained by the Hybrid Algorithm with Some Specific Influence Factors*
5.3.2. Optical property error
In BLT problem, the optical properties of the tissues are assumed to be determined by other techniques, such as diffuse optical tomography (DOT). To explore the influence of the optical property errors on our algorithm, different magnitude perturbations are added to the optical properties. As the errors of absorption and scattering coefficients in opposite directions may cancel out each other [8], we only take the optical properties all with positive or negative errors into account. In this part, optical properties in Table 1 are taken as baseline values for each component of the numerical phantom. The system matrix ℳ is computed once using these baseline values. External measurements are computed using optical properties that different from the baseline values by various amount. Positive/negetive errors mean the optical properties that be used to compute the measurements are greater/lower than the baseline values.
We have carried out simulations with +25%, +50%, +100%, −10%, −25%, −50% optical property errors, respectively. The numerical results indicate that the reconstructed source centers are off their actual locations to the phantom edge, and the reconstructed source powers are lower than the actual values in the positive optical error cases. In contrast, the reconstructed source centers are off their actual locations to the phantom center, and the reconstructed source powers are greater than the actual values in the negative optical error cases. The location and power errors increase with the increasing of the optical property errors. The reconstructed source powers are influenced more greatly than the reconstructed source centers by optical property errors. The power errors are even over 100% with large optical property errors, however the location errors are within 1mm in all the cases. Figure 4 shows the reconstructed results obtained by the hybrid algorithm with +50% optical property error. In this case the location error is within 0.80mm, and the power error is about 65% as it shown in Table 5.
5.3.3. Phantom structure error
The anatomical structure of the object is assumed to be provided by other tomographic technique, such as CT/micro-CT or MRI. To explore the influence of the phantom structure errors on the hybrid algorithm, we make some shifts to different components of the numerical phantom. In this part, the system matrix ℳ is computed once using the phantom in Fig. 2. External measurements are computed using the phantom with shifted components.
In our simulations both the heart (H) and the lungs (L) have been shifted with different distances (0.5mm, 1.0mm, 1.5mm), respectively. The numerical results indicate that when the lungs are shifted outward, the reconstructed source centers are off their actual locations to the phantom center, and the reconstructed source powers are lower than the actual values. In contrast, when the lungs are shifted inward, the reconstructed source centers are off their actual locations to the phantom edge, and the reconstructed source powers are greater than the actual values. The influence on the reconstructed results increases with the increasing of the phantom structure errors. The influence of the heart shift on the reconstructed results is comparatively less than that of the lungs shift, that because the sources are set in the lungs. Figure 5 shows the reconstructed results obtained by the hybrid algorithm with the lungs being shifted 1mm outward to the phantom edge. In this case the location error is within 0.39mm, and the power error is within 7% as it shown in Table 5.
Fig. 5 Reconstructed result obtained by the hybrid algorithm with the lungs being shifted 1mm outward to the phantom edge.
5.3.4. Mixed influence
We have explored the influence of the detection noises, optical property errors and phantom structure errors on the hybrid algorithm respectively. In this part, we explore the performance of the hybrid algorithm with the influence of the mixed factors. We have carried out a simulation with +50% optical property error, the lungs being shifted 1mm outward to the phantom edge, and the measurement being corrupted by Gaussian(0,100).
Figure 6 shows the reconstructed results obtained by the hybrid algorithm with the influence of the mixed factors. In this case the location error is within 0.37mm, the reconstructed source centers are off their actual locations to the phantom edge, and the power error is about 67%. Compared to the reconstructed results with the separate influence factors, the reconstructed result with the mixed factors is the most similar to that with +50% optical property error. As shown in Table 5, the location error is greater than 0.70mm with +50% optical property error, however that is only about 0.37mm with the influence of the mixed factors. That because the reconstructed source centers are off their actual locations to the phantom edge with the positive optical property error, and the influence is opposite with the lungs being shifted outward to the phantom edge, they partially canceled out each other which leads to the improved source localization. Consequently, the reconstructed results with the mixed factors is the comprehensive influence of all the factors. The numerical result indicates that the optical property error influences the reconstructed results most, the phantom structure error less, and the detection noise the least.
Fig. 6 Reconstructed result obtained by the hybrid algorithm with the influence of the mixed factors.
6. Conclusions and discussions
A generalized hybrid algorithm has been proposed based on the graph cuts algorithm and gradient-based algorithms. Here, graph cuts algorithm is adopted to estimate the source support and different gradient-based algorithms are sequentially used to reconstruct the source distribution according to the iteration status. Therefore, a more reliable source support can be obtained without prior knowledge, and a more accurate and finer source distribution can be finally acquired. Furthermore, multilevel meshes for the internal sources are used to enhance the performance of reconstruction. Three algorithmic issues, including estimation of source number and total power, deblurring of source distribution and stopping criteria, have also been investigated. Numerical simulations have been performed to validate this proposed algorithm and demonstrate its high performance in the multi-source situation even if the detection noises, optical property errors and phantom structure errors are involved in the forward imaging.
Generally, the system matrices corresponding to different source scales have been all calculated in advance. However, if a compact source support has been determined, only a few columns of the system matrices related to the nodes in the support need to be calculated, and this may reduce the computational cost to a great extent. Therefore, it implies that the system matrices in the step of source distribution computation can be calculated off-the-cuff.
Nevertheless with a good estimation of the source location and power, the proposed algorithm does not work well for the estimation of source size and intensity distribution. This is because that the sources with different sizes and different intensity distributions can generate the same boundary measurement. To improve the uniqueness of the solution, the multi-spectral technique and some other effective approaches, such as the eigenvectors expansion approach [24], could be introduced into our generalized hybrid algorithm scheme.
Besides, due to its tolerance of phantom structure errors and optical property errors, our algorithm can be probably used for the bioluminescent source reconstruction in the multi-modality imaging system.
Acknowledgments
This work is supported in part by the National Basic Research Program of China (973 Program) (2011CB809105), the National Science Foundation of China (61101156, 61121002, 10990013), Key Laboratory of Machine Perception (Ministry of Education) of Peking University, and Microsoft Research of Asia. The authors would like to thank Yuanzheng Si for valuable discussions.
References and links
1. C. H. Contag and M. H. Bachmann, “Advances in in vivo bioluminescence imaging of gene expression,” Annu. Rev. Biomed. Eng. 4, 235–260 (2002) [CrossRef] [PubMed] .
2. 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. Imaging 16, 378–387 (2002) [CrossRef] [PubMed] .
3. V. Ntziachristos, J. Ripoll, L. V. Wang, and R. Weissleder, “Looking and listening to light: the evolution of whole-body photonic imaging,” Nat. Biotechnol. 23, 313–320 (2005) [CrossRef] [PubMed] .
4. G. Wang, E. A. Hoffman, G. McLennan, L. V. Wang, M. Suter, and J. Meinel, “Development of the first bioluminescent CT scanner,” Radiology 229(P), 566 (2003).
5. 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 (2000) [CrossRef] .
6. M. Rudin and R. Weissleder, “Molecular imaging in drug discovery and development,” Nat. Rev. Drug Discovery 2, 123–131 (2003) [CrossRef] .
7. 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] .
8. G. Alexandrakis, F. R. Rannou, and A. F. Chatziioannou, “Effect of optical property estimation accuracy on tomographic bioluminescence imaging: simulation of a combined optical-PET (OPET) system,” Phys. Med. Biol. 51, 2045–2053 (2006) [CrossRef] [PubMed] .
9. G. Wang, H. Shen, Y. Liu, A. Cong, W. Cong, Y. Wang, and P. Dubey, “Digital spectral separation methods and systems for bioluminescence imaging,” Opt. Express 16, 1719–1732 (2008) [CrossRef] [PubMed] .
10. 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) [CrossRef] [PubMed] .
11. X. Ma, J. Tian, C. Qin, X. Yang, B. Zhang, Z. Xue, X. Zhang, D. Han, D. Dong, and X. Liu, “Early detection of liver cancer based on bioluminescence tomography,” Appl. Opt. 50, 1389–1395 (2011) [CrossRef] [PubMed] .
12. M. Jiang and G. Wang, “Image reconstruction for bioluminescence tomography,” Proc. SPIE 5535, 335–351 (2004) [CrossRef] .
13. 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) [CrossRef] [PubMed] .
14. 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) [CrossRef] [PubMed] .
15. 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) [CrossRef] [PubMed] .
16. S. Ahn, A. J. Chaudhari, F. Darvas, C. A. Bouman, and R. M. Leahy, “Fast iterative image reconstruction methods for fully 3D multispectral bioluminescence tomography,” Phys. Med. Biol. 53, 3921–3942 (2008) [CrossRef] [PubMed] .
17. Y. Lu, X. Zhang, A. Douraghy, D. Stout, J. Tian, T. F. Chan, and A. F. Chatziioannou, “Source reconstruction for spectrally-resolved bioluminescence tomography with sparse a priori information,” Opt. Express 17, 8062–8080 (2009) [CrossRef] [PubMed] .
18. B. Zhang, X. Yang, C. Qin, D. Liu, S. Zhu, J. Feng, L. Sun, K. Liu, D. Han, X. Ma, X. Zhang, J. Zhong, X. Li, X. Yang, and J. Tian, “A trust region method in adaptive finite element framework for bioluminescence tomography,” Opt. Express 18, 6477–6491 (2010) [CrossRef] [PubMed] .
19. K. Liu, J. Tian, Y. Lu, C. Qin, X. Yang, S. Zhu, and X. Zhang, “A fast bioluminescent source localization method based on generalized graph cuts with mouse model validations,” Opt. Express 18, 3732–3745 (2010) [CrossRef] [PubMed] .
20. G. Wang, Y. Li, and M. Jiang, “Uniqueness theorems in bioluminescence tomography,” Med. Phys. 31, 2289–2299 (2004) [CrossRef] [PubMed] .
21. 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) [CrossRef] [PubMed] .
22. 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) [CrossRef] [PubMed] .
23. M. A. Naser and M. S. Patterson, “Algorithms for bioluminescence tomography incorporating anatomical information and reconstruction of tissue optical properties,” Biomed. Opt. Express 1, 512–526 (2010) [CrossRef] .
24. M. A. Naser and M. S. Patterson, “Bioluminescence tomography using eigenvectors expansion and iterative solution for the optimized permissible source region,” Biomed. Opt. Express 2, 3179–3193 (2011) [CrossRef] [PubMed] .
25. H. Dehghani, S. C. Davis, and B. W. Pogue, “Spectrally resolved bioluminescence tomography using the reciprocity approach,” Med. Phys. 35, 4863–4871 (2008) [CrossRef] [PubMed] .
26. 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) [CrossRef] [PubMed] .
27. S. R. Arridge, “Optical tomography in medical imaging,” Inverse Prob. 15, R41–R93 (1999) [CrossRef] .
28. F. Natterer and F. Wübbeling, Mathematical Methods in Image Reconstruction (SIAM, 2001) [CrossRef] .
29. S. R. Arridge and M. Schweiger, “Photon-measurement density functions. Part 2: Finite-element-method calculations,” Appl. Opt. 34, 8026–8037 (1995) [CrossRef] [PubMed] .
30. M. Piana and M. Bertero, “Projected landweber method and preconditioning,” Inverse Prob. 13, 441–463 (1997) [CrossRef] .
31. M. Jiang and G. Wang, “Convergence studies on iterative algorithms for image reconstruction,” IEEE Trans. Med. Imaging 22, 569–579 (2003) [CrossRef] [PubMed] .
32. C. Byrne, “Iterative oblique projection onto convex sets and the split feasibility problem,” Inverse Prob. 18, 441–453 (2002) [CrossRef] .
33. P. E. Gill, W. Murray, and M. H. Wright, Practical Optimization (Academic Press, 1981).
34. Y. Boykov and V. Kolmogorov, “An experimental comparison of min-cut/max-flow algorithms for energy minimization in vision,” IEEE Trans. Pattern Anal. Mach. Intell. 26, 1124–1137 (2004) [CrossRef] .
35. V. Kolmogorov and R. Zabin, “What energy functions can be minimized via graph cuts?” IEEE Trans. Pattern Anal. Mach. Intell. 26, 147–159 (2004) [CrossRef] [PubMed] .
36. P. L. Hammer, P. Hansen, and B. Simeone, “Roof duality, complementation and persistency in quadratic 0–1 optimization,” Math. Program. 28, 121–155 (1984) [CrossRef] .
37. M. Sonka, V. Hlavac, and R. Boyle, Image Processing, Analysis, and Machine Vision (Thompson Learning, 2008).
