Computational and Mathematical Methods in Medicine

Volume 2013 (2013), Article ID 767296, 8 pages

http://dx.doi.org/10.1155/2013/767296

## Improved Reconstruction Quality of Bioluminescent Images by Combining SP_{3} Equations and Bregman Iteration Method

College of Electronic Information and Control Engineering, Beijing University of Technology, Beijing 100124, China

Received 13 September 2012; Accepted 24 December 2012

Academic Editor: Jimin Liang

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

#### Abstract

Bioluminescence tomography (BLT) has a great potential to provide a powerful tool for tumor detection, monitoring tumor therapy progress, and drug development; developing new reconstruction algorithms will advance the technique to practical applications. In the paper, we propose a BLT reconstruction algorithm by combining SP_{3} equations and Bregman iteration method to improve the quality of reconstructed sources. The numerical results for homogeneous and heterogeneous phantoms are very encouraging and give significant improvement over the algorithms without the use of SP_{3} equations and Bregman iteration method.

#### 1. Introduction

As an emerging molecular imaging technique, bioluminescence imaging (BLI) is potentially well suited for early detection, clinical drug development and monitoring, and regeneration research [1–5]. Therefore, this imaging modality has received increasingly intense research interest worldwide over the recent years.

To date, planar BLI is commonly used because of its ease of implementation and operational simplicity, but it also suffers from significant limitations, including the low resolution, the lack of quantification, and the incapacity of accurately providing depth information [6]. In contrast, bioluminescence tomography (BLT) could overcome these limitations by using accurate reconstruction algorithms coupled with theoretical models of photon propagation in biological tissues, providing higher resolution, quantification accuracy, and depth information [7]. In comparing BLT to planar BLI, planar BLI is a qualitative analysis and BLT is a quantitative analysis [8]. Therefore, scientists are now paying more attention to the advancement of BLT research.

The objective of BLT is to recover the unknown bioluminescent source distribution based on the noisy surface measurements [6, 7]. Indeed, the problem is also called the inverse problem. However, a major difficulty in recovering the bioluminescent source distribution is imposed by multiple scattering which occurs when light propagates through biological tissues. This makes the inverse problem severely ill-posed [7]. Furthermore, the number of recovered unknown source distributions is usually far more than the number of detected boundary measurements, that is, (in many cases, ). Hence, BLT is also a typically underdetermined problem. To obtain a meaningful solution, regularization techniques are usually adopted, which consist of solving the following constrained optimization problem [9]: where is a properly chosen regularization term, represents regularization parameter, and is a linear operator, typically formed by discretizing diffusion equation with finite element methods [10].

When , the above regularized problem becomes the popular Tikhonov regularization, which inherently provides smoothed solutions and therefore offers compromised accuracy in localizing bioluminescent sources [11]. Recently, -regularized problems, that is, , have received an increasing amount of attention in optical imaging, which allow high-quality images to be reconstructed from a small amount of boundary measurements [11–14]. However, -regularized problems can sparsify the bioluminescent source distribution, which affects the quality of reconstructed images [13, 15].

Furthermore, in order to obtain the matrix in (1), the diffusion approximation (DA) to radiative transfer equation (RTE) is widely used as the forward model for BLT reconstructions. Although the DA is one of the most important approximation methods in BLT [6–11], it suffers from some limitations [12–14]. Firstly, the scattering is dominated over absorption and secondly, the DA fails in modeling light propagation in the vicinity of those highly vascularized tissue parts [12–14]. Therefore, the DA will introduce significant error in some BLT cases [14]. In contrast, the RTE is widely accepted as an accurate model for light propagation in biological tissues. However, the use of the RTE as the forward model for BLT is often not feasible due to the facts that analytical solutions cannot exist for biological tissues with spatially nonuniform scattering and absorption properties and the computation of numerical approximations for the solution is extremely time consuming [16, 17]. A generalized delta-Eddington phase function was recently presented to simplify the RTE, and the more accurate solution was obtained relative to the DA [18, 19]. However, the parameter used in the model is difficult to compute [18, 19]. In addition, the system matrix for the model is also difficult to construct for complex heterogeneous geometries. These factors seriously limit the utilization of the model in BLT. The use of simplified spherical harmonics (SP_{N}) equations to approximate the RTE has been demonstrated to significantly improve the diffusion solution in domains with high absorption and small geometries [5, 12–14, 16, 20]. Meanwhile, the SP_{N} methods are computationally less expensive than the RTE ones.

Large efforts in combining multiple types of *a priori* information to develop BLT reconstruction algorithms to improve the quality of reconstructed images, particularly the permissible source region and multispectral information, have formed the grounds of BLT reconstructions [9–11, 20–26]. Despite the recent advances in BLT reconstruction algorithms and light propagation models, it is necessary to develop and refine reconstruction methods to improve image quality.

Bregman iteration method has been studied recently and is widely used in compressed sensing [27, 28]. The idea is to add the residual, that is, the error produced at the current iteration, back to the data for the next iteration to be corrected [27]. The method is particularly attractive for sparse reconstruction, but so far it has not been fully investigated and analyzed in BLT, and this is the goal of this paper.

In this paper, we propose a BLT algorithm to improve the quality of reconstructed images. In the algorithm, SP_{3} equations are adapted to model light propagation, and Bregman iteration method is used to solve the inverse problem for BLT. Numerical results demonstrate that the quality of reconstructed images is improved greatly. The rest of the paper is organized as follows. In the following section, we described SP_{3} equations as light propagation model and Bregman iteration method. Last, numerical experiments were performed to evaluate the proposed algorithm, and corresponding conclusions were made.

#### 2. Methods

##### 2.1. SP_{3} Equations

The propagation of light in biological tissues can be well modeled by SP_{3} equations. SP_{3} equations are two coupled diffusion equations for the moments and [16, 17]:
where , and and are the absorption and scattering parameters, respectively. is the anisotropy parameter.

The boundary conditions are given by

The coefficients can be found in [16]. Furthermore, the partial current can be obtained from solutions and :

The coefficients can also be found in [16]. Solving the above equations by finite element methods, a linear operator can be established [29].

##### 2.2. Bregman Iteration Method

Bregman iteration method is based on the definition of Bregman distance. The Bregman distance associated with a convex function at the point is given as [27] where is in the subgradient of at . Clearly, this is not a distance in the usual sense because it is not in general symmetric. However, it does measure closeness in the sense that and for on the line segment between and [27].

Based on Bregman iteration method, (1) can be reformulated as where and is the adjoint operator of . Since the operator is linear in BLT reconstructions, the above complicated iteration can be transformed to the following two-stage iteration procedure with [27]:

This is done by iteratively solving the optimization problem (7) and then modifying the measured value of used in the next iteration. And (8) is usually referred as “adding back the noise” [30]. In the paper, is fixed as the regularizer. The implementation of (7) was performed by a gradient projected (GP) algorithm [31]. The proposed algorithm was depicted in Algorithm 1.

#### 3. Results

To fully evaluate the performance of the proposed algorithm, homogeneous and heterogeneous experiments were performed. In the experiments, the parameters and were set to and 10, respectively. The parameters in GP algorithm set default values, except the maximum iteration number is fixed at 50000 to ensure the convergence of the algorithm unless otherwise is specified.

##### 3.1. Homogeneous Phantom Experiments

In this section, 2D numerical simulations were used to investigate the performance of the proposed algorithm since less computational time was required for 2D data. Here, two individual cases were considered. In the first case, numerical simulations were performed on a homogenous circle with 10 mm radius. Within this circle, two sources (source 1 and source 2) were placed in (−5, 0) mm and (0, 5) mm, respectively and each source had a radius of 1.0 mm. The corresponding optical parameters were listed in Table 1. The boundary data were generated for two wavelengths (600 and 620 nm) with finite element methods, and different levels of Gaussian noise (0%, 10%, and 30%) were added to the datasets. BLT reconstructions were performed without and with Bregman iteration method. Corresponding results were shown in Figure 1. In this case, the ratios of are larger than 10; therefore, the circular phantom has high-scattering characteristics. Hence, the DA is suitable for the simulation. For comparison, we carried out BLT reconstructions with the DA as the forward model; reconstructed images were also illustrated in Figure 1. From Figure 1, we can see that the results with SP_{3} equations are better than those obtained with the DA and Bregman iteration method can improve the quality of reconstructed images. The best results are obtained by combing SP_{3} equations and Bregman iteration method. In addition, quantitative results were summarized in Table 2. Data in Table 2 show that reconstructed position errors can be significantly reduced when SP_{3} equations are used together with Bregman iteration method.

Furthermore, we tested the proposed algorithm by using experiments with multiple bioluminescent sources. The optical properties of a real mouse muscle for different wavelengths (580 and 620 nm) were assigned as listed in Table 3 [29]. Four identical sources with 1 mm radii were placed different positions. First, the sources were placed near the surfaces, and the distance to the center of the circle was 7.07 mm. The boundary measurements were also produced by finite element methods, and 20% Gaussian noise was added into the simulated data. Note that in the test, for two wavelengths are less than 10; therefore, the condition does not hold and the DA is less valid. Hence, BLT reconstructions with the DA were not implemented. The results with SP_{3} equations are shown in Figures 2(a) and 2(b). Next, the sources were placed at 5 mm positions off the center. Then BLT reconstructions were performed, as shown in Figures 2(c) and 2(d). Furthermore, quantitative results were shown in Table 4. It is worthy of mentioning that BLT reconstructions without and with Bregman iteration method use the same regularization parameter (i.e., ), but the reconstructed results are different. From Figure 2 and Table 4, it is easily concluded that better images can be obtained by combining SP_{3} equations and Bregman iteration method.

##### 3.2. Heterogeneous Phantom

In the subsection, a micro-MRI-based heterogeneous mouse model (MOBY) was used to validate the proposed algorithm [32]. About 2/3 of the entire phantom was used for mesh generation, and a volumetric mesh with 17661 nodes and 93312 tetrahedron elements was obtained by iso2mesh [33], as shown in Figure 3. The optical properties of different tissues were assigned according to Table 5, reproduced from Alexandrakis et al. [21]. The forward simulation data was produced by finite element methods, and 10% Gaussian noise was added. Then BLT reconstructions were performed without and with Bregman iteration method. The regularization parameters used in the two methods were the same, and the value was 0.1. The maximum iteration number in the GP algorithm was set to 5000, and other parameters remained unchanged. The reconstructed results without and with Bregman iteration method were shown in Figure 4. From the images, we can see that the quality of reconstructed images can be improved with the use of Bregman iteration method. Furthermore, the reconstructed central positions for the two algorithms are (22.77, 14.95, 13.33 mm) and (22.24, 13.95, 14.49 mm), respectively. The real source position is (22.07, 14.43, 13.06 mm). The absolute distances between reconstructed sources and the real source are 0.91 mm and 1.52 mm, respectively. The quantitative results also demonstrate that Bregman iteration method can improve the quality of reconstructed images.

#### 4. Conclusion

We have presented a BLT reconstruction algorithm by combing SP_{3} equations and Bregman iteration method as a competitive method for reconstructing bioluminescent sources and validated the proposed algorithm using homogeneous and heterogeneous experiments. It has been demonstrated that the proposed algorithm can enhance the recovery of bioluminescent sources in terms of the quality of reconstructed images and localization error.

The use of SP_{3} equations is a helpful technique to improve BLT reconstructions. Our experiments have illustrated that the appearance of artifacts can be reduced when SP_{3 }equations are used as the forward model. However, the computation of the system matrix by solving SP_{3} equations is very expensive, especially when the imaged objects are very complex, irregular, and heterogeneous. Fortunately, with the fast development of graphics processing unit (GPU), the computation of can be significantly accelerated.

One merit of the proposed algorithm is that the improved results are obtained by making use of the available boundary measurements and thus do not require increased number of boundary measurements and do not bring more hardware requirements. Meanwhile, the proposed algorithm is relatively easy to implement. Therefore, the algorithm is suitable for *in vivo* applications. As a sacrifice, the computational burden for the proposed algorithm is greatly increased, especially for the heterogeneous mouse experiment, since solving (1) brings extra cost through Bregman iteration method, and each iteration of which is equivalent of solving a standard “L1” problem. To increase computational efficiency for mouse experiments, developing fast large-scale optimization algorithms is essential.

In conclusion, we have developed a BLT reconstruction algorithm by combing SP_{3} equations and Bregman iteration method and indicated its feasibility and merits. In the near future, we expect to accelerate the proposed algorithm based on GPU and extend it to *in vivo* mouse experiments.

#### Acknowledgments

This paper is supported by the National Natural Science Foundation of China under Grant nos. 30970780 and 81000624, the Science and Technology Key Project of Beijing Municipal Education Commission under Grant no. KZ200910005005, the Doctoral Fund of the Ministry of Education of China under Grant no. 20091103110005, the Rixin Fund of Beijing University of Technology under Grant no. 2013-RX-L04", and the Scientific Research Fund of Beijing University of Technology under Grant no. X0002012201101.

#### References

- V. Ntziachristos, J. Ripoll, L. V. Wang, and R. Weissleder, “Looking and listening to light: the evolution of whole-body photonic imaging,”
*Nature Biotechnology*, vol. 23, no. 3, pp. 313–320, 2005. View at Publisher · View at Google Scholar · View at Scopus - R. Weissleder and M. J. Pittet, “Imaging in the era of molecular oncology,”
*Nature*, vol. 452, no. 7187, pp. 580–589, 2008. View at Publisher · View at Google Scholar · View at Scopus - J. K. Willmann, N. van Bruggen, L. M. Dinkelborg, and S. S. Gambhir, “Molecular imaging in drug development,”
*Nature Reviews Drug Discovery*, vol. 7, no. 7, pp. 591–607, 2008. View at Publisher · View at Google Scholar · View at Scopus - T. Schroeder, “Imaging stem-cell-driven regeneration in mammals,”
*Nature*, vol. 453, no. 7193, pp. 345–351, 2008. View at Publisher · View at Google Scholar · View at Scopus - J. Feng, C. Qin, K. Jia, S. Zhu, X. Yang, and J. Tian, “Bioluminescence tomography imaging
*in vivo*: recent advances,”*IEEE Journal of Selected Topics in Quantum Electronics*, vol. 18, no. 4, pp. 1394–1402, 2012. View at Google Scholar - X. Gu, Q. Zhang, L. Larcom, and H. Jiang, “Three-dimensional bioluminescence tomography with model-based reconstruction,”
*Optics Express*, vol. 12, no. 17, pp. 3996–4000, 2004. View at Publisher · View at Google Scholar · View at Scopus - G. Wang, Y. Li, and M. Jiang, “Uniqueness theorems in bioluminescence tomography,”
*Medical Physics*, vol. 31, no. 8, pp. 2289–2299, 2004. View at Publisher · View at Google Scholar · View at Scopus - G. Wang, H. Shen, W. Cong, S. Zhao, and G. W. Wei, “Temperature-modulated bioluminescence tomography,”
*Optics Express*, vol. 14, no. 17, pp. 7852–7871, 2006. View at Publisher · View at Google Scholar · View at Scopus - S. Ahn, A. J. Chaudhari, F. Darvas, C. A. Bouman, and R. M. Leahy, “Fast iterative image reconstruction methods for fully 3D multispectral bioluminescence tomography,”
*Physics in Medicine and Biology*, vol. 53, no. 14, pp. 3921–3942, 2008. View at Publisher · View at Google Scholar · View at Scopus - W. Cong, G. Wang, D. Kumar et al., “Practical reconstruction method for bioluminescence tomography,”
*Optics Express*, vol. 13, no. 18, pp. 6756–6771, 2005. View at Publisher · View at Google Scholar · View at Scopus - Y. Lu, X. Zhang, A. Douraghy et al., “Source reconstruction for spectrally-resolved bioluminescence tomography with sparse a priori information,”
*Optics Express*, vol. 17, no. 10, pp. 8062–8080, 2009. View at Publisher · View at Google Scholar · View at Scopus - H. Gao and H. Zhao, “Multilevel bioluminescence tomography based on radiative transfer equation part 1: l1 regularization,”
*Optics Express*, vol. 18, no. 3, pp. 1854–1871, 2010. View at Publisher · View at Google Scholar · View at Scopus - H. Gao and H. Zhao, “Multilevel bioluminescence tomography based on radiative transfer equation part 2: total variation and l1 data fidelity,”
*Optics Express*, vol. 18, no. 3, pp. 2894–2912, 2010. View at Publisher · View at Google Scholar · View at Scopus - K. Liu, J. Tian, C. Qin et al., “Tomographic bioluminescence imaging reconstruction via a dynamically sparse regularized global method in mouse models,”
*Journal of Biomedical Optics*, vol. 16, no. 4, Article ID 046016, 2011. View at Google Scholar - J. Feng, C. Qin, K. Jia et al., “Total variation regularization for bioluminescence tomography with the split Bregman method,”
*Applied Optics*, vol. 51, no. 19, pp. 4501–4512, 2012. View at Google Scholar - A. D. Klose and E. W. Larsen, “Light transport in biological tissue based on the simplified spherical harmonics equations,”
*Journal of Computational Physics*, vol. 220, no. 1, pp. 441–470, 2006. View at Publisher · View at Google Scholar · View at Scopus - Z. Yuan, X. H. Hu, and H. Jiang, “A higher order diffusion model for three-dimensional photon migration and image reconstruction in optical tomography,”
*Physics in Medicine and Biology*, vol. 54, no. 1, pp. 65–88, 2009. View at Publisher · View at Google Scholar · View at Scopus - W. Cong, A. Cong, H. Shen, Y. Liu, and G. Wang, “Flux vector formulation for photon propagation in the biological tissue,”
*Optics Letters*, vol. 32, no. 19, pp. 2837–2839, 2007. View at Publisher · View at Google Scholar · View at Scopus - W. Cong, H. Shen, A. Cong, Y. Wang, and G. Wang, “Modeling photon propagation in biological tissues using a generalized Delta-Eddington phase function,”
*Physical Review E*, vol. 76, no. 5, Article ID 051913, 2007. View at Publisher · View at Google Scholar · View at Scopus - A. D. Klose, B. J. Beattie, H. Dehghani et al., “
*In vivo*bioluminescence tomography with a blocking-off finite-difference SP_{3}method and MRI/CT coregistration,”*Medical Physics*, vol. 37, no. 1, pp. 329–338, 2010. View at Publisher · View at Google Scholar · View at Scopus - G. Alexandrakis, F. R. Rannou, and A. F. Chatziioannou, “Tomographic bioluminescence imaging by use of a combined optical-PET (OPET) system: a computer simulation feasibility study,”
*Physics in Medicine and Biology*, vol. 50, no. 17, pp. 4225–4241, 2005. View at Publisher · View at Google Scholar · View at Scopus - H. Dehghani, S. C. Davis, S. Jiang, B. W. Pogue, K. D. Paulsen, and M. S. Patterson, “Spectrally resolved bioluminescence optical tomography,”
*Optics Letters*, vol. 31, no. 3, pp. 365–367, 2006. View at Publisher · View at Google Scholar · View at Scopus - Y. Lv, J. Tian, W. Cong et al., “A multilevel adaptive finite element algorithm for bioluminescence tomography,”
*Optics Express*, vol. 14, no. 18, pp. 8211–8223, 2006. View at Publisher · View at Google Scholar · View at Scopus - J. Feng, K. Jia, G. Yan et al., “An optimal permissible source region strategy for multispectral bioluminescence tomography,”
*Optics Express*, vol. 16, no. 20, pp. 15640–15654, 2008. View at Publisher · View at Google Scholar · View at Scopus - J. Feng, K. Jia, C. Qin et al., “Three-dimensional bioluminescence tomography based on Bayesian approach,”
*Optics Express*, vol. 17, no. 19, pp. 16834–16848, 2009. View at Publisher · View at Google Scholar · View at Scopus - J. Feng, C. Qin, K. Jia et al., “An adaptive regularization parameter choice strategy for multi-spectral bioluminescence tomography,”
*Medical Physics*, vol. 38, no. 11, pp. 5933–5944, 2011. View at Google Scholar - T. Goldstein and S. Osher, “The split Bregman method for l1-regularized problems,”
*SIAM Journal on Imaging Sciences*, vol. 2, no. 2, pp. 323–343, 2009. View at Google Scholar - S. Osher, Y. Mao, B. Dong, and W. Yin, “Fast linearized Bregman iteration for compressive sensing and sparse denoising,”
*UCLA CAM Reports*08-37, 2008. View at Google Scholar - Y. Lu, A. Douraghy, H. B. Machado et al., “Spectrally resolved bioluminescence tomography with the third-order simplified spherical harmonics approximation,”
*Physics in Medicine and Biology*, vol. 54, no. 21, pp. 6477–6493, 2009. View at Publisher · View at Google Scholar · View at Scopus - S. Osher, M. Burger, D. Goldfarb, J. Xu, and W. Yin, “An iterative regularization method for total variation-based image restoration,”
*Multiscale Modeling and Simulation*, vol. 4, no. 2, pp. 460–489, 2005. View at Publisher · View at Google Scholar · View at Scopus - M. A. T. Figueiredo, R. D. Nowak, and S. J. Wright, “Gradient projection for sparse reconstruction: application to compressed sensing and other inverse problems,”
*IEEE Journal on Selected Topics in Signal Processing*, vol. 1, no. 4, pp. 586–597, 2007. View at Publisher · View at Google Scholar · View at Scopus - W. P. Segars, B. M. W. Tsui, E. C. Frey, G. A. Johnson, and S. S. Berr, “Development of a 4-D digital mouse phantom for molecular imaging research,”
*Molecular Imaging and Biology*, vol. 6, no. 3, pp. 149–159, 2004. View at Publisher · View at Google Scholar · View at Scopus - Q. Fang and D. A. Boas, “Tetrahedral mesh generation from volumetric binary and grayscale images,” in
*Proceedings of the IEEE International Symposium on Biomedical Imaging: From Nano to Macro (ISBI '09)*, pp. 1142–1145, Boston, Mass, USA, July 2009. View at Publisher · View at Google Scholar · View at Scopus