## Abstract

Fluorescence molecular tomography (FMT) has become an important method for in-vivo imaging of small animals. It has been widely used for tumor genesis, cancer detection, metastasis, drug discovery, and gene therapy. In this study, an algorithm for FMT is proposed to obtain accurate and fast reconstruction by combining an adaptive mesh refinement technique and an analytical solution of diffusion equation. Numerical studies have been performed on a parallel plate FMT system with matching fluid. The reconstructions obtained show that the algorithm is efficient in computation time, and they also maintain image quality.

© 2007 Optical Society of America

## Corrections

Daifa Wang, Xiaolei Song, and Jing Bai, "Adaptive-mesh-based algorithm for fluorescence molecular tomography using an analytical solution: erratum," Opt. Express**15**, 10789-10789 (2007)

https://www.osapublishing.org/oe/abstract.cfm?uri=oe-15-17-10789

## 1. Introduction

Molecular imaging, especially small-animal imaging, has become an important method for biomedical research. It has been widely used for tumor genesis, cancer detection, metastasis, drug discovery, and gene therapy [1, 2]. Fluorescence molecular tomography (FMT) has become a promising imaging modality of molecular imaging. This is because it is low-cost, sensitive, and involves fluorescence/bioluminescence resonance energy transfer processes among specific molecular interactions.

In the past years, rapid progress has been made in the reconstruction of FMT. Several reconstruction algorithms for FMT have been proposed such as Born-type approximation techniques [3–5], Newton’s or Newton-type optimization methods [6–8], Bayesian nonlinear least squares approaches [9,10], adaptive finite-element-based method [11–13], and finiteelement- based method with dual excitation mode [14].

In all the reconstruction methods proposed, the volume of investigation is discretized to small voxels for computation. Hence, the quality of FMT reconstruction is not only determined by the signal-to-noise ratio (SNR) but also by the refined mesh of the volume. To a large extent, the finer the refined mesh is, the better the resolution. Usually, the refined mesh is implemented uniformly all over the volume of investigation, but the global uniform discretization always results in an overly detailed mesh, which increases the numerical instability and requires more computational resources. Therefore, adaptive mesh refinement is very useful, providing finer meshes around target locations and coarser meshes in other regions. Two adaptive mesh algorithms based on FEM have been reported previously, one for FMT by Joshi *et al*. [11] and the other for bioluminescence tomography (BLT) by Lv *et al*. [12]. Both methods have shown the abilities of improving image quality and decreasing computational cost. But, since they are based on FEM, the computational cost is still big. On the other hand, analytical solution of diffusion equation is two orders of magnitude faster [15] than numerical methods such as FEM.

In this paper, the advantages of the adaptive mesh method and analytical solution are combined in order to reduce the computation time while maintaining reconstructed image quality. In the algorithm proposed, an analytical solution of diffusion equation for slab geometry [16] is used to calculate light propagation, and adaptive mesh refinement is performed during reconstruction. In the adaptive procedure, the optimization is implemented using an interior point conjugate gradient method with trust region (Interior/CG) [17] to determine the fluorescence map that best predicts the measured boundary fluorescence.

The outline of this article is presented as follows. In Section 2, a novel adaptive-mesh-based algorithm for FMT using an analytical solution is presented. Section 3 details the numerical experiments for investigating the computational capacity of the algorithm proposed by reconstructing single and dual fluorescent targets. In Section 4, experiment results are summarized and discussed. Finally, we discuss and conclude our study in Section 5.

## 2. Methods

#### 2.1 Forward problem

When fluorescent objects are embedded in turbid media, the fluorescence measurements *U _{m}*(

*r*,

_{s}*r*) detected at position

_{d}*r*by a steady light source at position

_{d}*r*with the emission wavelength

_{s}*λ*can be formulated as [4,18]

_{m}where ${D}^{{\lambda}_{m}}$ is the diffusion coefficient with the emission wavelength *λ _{m}*. ${k}^{{\lambda}_{x}}$ is the wave number with the excitation wavelength

*λ*, and ${k}^{{\lambda}_{m}}$ is the wave number with the emission wavelength

_{x}*λ*. ${U}_{x}({r}_{s},r,{k}^{{\lambda}_{x}})$ is the theoretically calculated average photon intensity at position

_{m}*r*, and $G({r}_{d}-r,{k}^{{\lambda}_{m}})$ is the Green’s function that describes photon propagation from a point

*r*to the detector at location

*r*.

_{d}*S*

_{0}is a calibration factor that accounts for the system error such as laser power and the unknown gain and attenuation factors of the system.

*n*(

*r*) is the fluorescence distribution representing the fluorochrome concentration at location

*r*multiplied by the fluorescent yield. For a single source-detector pair (

*s*

_{1},

*d*

_{1}), Eq. (1) can be discretized in vector form as follows:

where *N* is the number of the discretized voxels. For *M* source-detector pairs, the weight matrix is generated as

where *W _{ij}* is an element of the weight matrix

*W*representing the weight of voxel

*j*to measurement

*i*.

#### 2.2 Reconstruction method

The level of discretization is important for the achievable spatial resolution of FMT. Image quality can be improved by uniformly refining the level of discretization. However, the global refinement increases the computational cost including memory and time by increasing the number of unknowns. A fine mesh is only necessary where small cells are needed to accurately resolve the fluorescent targets. Therefore, the technique of adaptive mesh refinement is adopted for multilevel meshing, which has been studied over two decades, especially in adaptive finite element analysis [19].

Let {Θ_{1},Θ_{2},…,Θ* _{k}*} represent a sequence of hexahedron mesh levels of the volume of investigation, in which the mesh changes from coarse to fine with the increase in the mesh level

*k*. Therefore, Eq. (3) on each level can be written as

where the index *k* represents the kth level. Practically, weight matrix *W _{k}* at the

*k*th level is usually ill-posed, because of the inherent property of diffusion equation and the existence of noise. Therefore, it is usually impractical to solve for fluorescence distribution

*n*from the linear system directly. In the algorithm proposed, an optimization approach is employed by minimizing objective function Ψ

_{k}*(*

_{k}*n*) based on the Tikhonov regularization method as follows:

_{k}where *U _{m}*(

*r*,

_{s}*r*) is the fluorescence measurements and

_{d}*λ*is the regularization parameter.

_{k}*n*and

^{L}_{k}*n*are the lower and upper bounds of the fluorescence distribution at the

^{U}_{k}*k*th level. Ψ

*(*

_{k}*n*) is optimized using Interior/CG [17]. Then, Eq. (5) can be written as

_{k}$$s.t.\phantom{\rule{.9em}{0ex}}{n}_{k}-{n}_{k}^{L}-s=0\phantom{\rule{.9em}{0ex}},$$

$${n}_{k}^{U}-{n}_{k}-s=0$$

where *s* is a vector of slack variables and *ν*>0. The optimization consists of finding solutions of Eq. (6) for a sequence of positive barrier parameters {*ν _{l}*} that converges to zero. When switching to a finer mesh, the final solution

*n*of the

^{r}_{k}*k*th level is used to generate the initial value

*n*

^{0}

_{k+1}of the (

*k*+1)th level as follows:

Then, the optimization on (*k*+1)th level is performed near the solution, which improves the reconstruction stability and accelerates the convergence speed.

A *posteriori* error estimates play an important role in adaptive mesh refinement. These estimates use the computed fluorescence distribution on a coarse mesh to indicate the cells where mesh refinement will be most beneficial for the reduction of measurement error. In this work, the direct maximum selection method [12] is adopted for mesh refinement. On a coarse mesh, the elements with greater distribution value most likely represent targets’ locations. They are where mesh refinement should be performed. Elements selected are then refined via dividing into eight equal son hexahedron elements that have the same fluorescence distributions with their father elements. In the implementation, 40% of the elements are selected for refinement each time.

#### 2.3 Software implementation

At the beginning of the algorithm, an initial coarse mesh is refined with mesh level *k*=1. Then, on each mesh level, optimization is performed until either (i) the number of major interior/CG iterations exceeds the maximum iteration number *L*
_{max}, or (ii) the relative change in the parameter estimates between the last two steps is less than the predefined threshold ξ. Mesh refinement is then performed. The process iterates until the mesh level exceeds the maximum mesh level *K*
_{max}. In this work, *L*
_{max}=40, ξ=10^{-7}, and *K*
_{max}=3 were used.

Most parts of the algorithm proposed were coded in C++ for computational speed and the complex data structure needed for adaptive mesh refinement, and the others were coded in Matlab for flexibility. In this paper, analytical solution of diffusion equation for slab geometry [16] was used for the generation of weight matrix *W _{k}* on each mesh level. Therefore, the algorithm is suitable for the parallel plate FMT system with matching fluid as described in Grave’s study [5]. Reconstructions were carried out on a Pentium 4, 3GHz and 2GB RAM computer.

## 3. Computational experiments

A series of computational experiments were designed to evaluate the efficiency of the proposed algorithm on the parallel plate FMT system with matching fluid, which was modeled after an existing system reported in Graves *et al*. [5]. As shown in Fig. 1(a), a laser source was used as the light source, which was decoupled into optical fibers by an optical switch. The imaging chamber filled with optical matching fluid was illuminated sequentially by the optical fibers. Fluorescence was detected using a highly sensitive charge-coupled device (CCD) after propagating through a band-pass optic filter.

The synthetic measurements were generated on a compressed parallel-plate imaging chamber of which the width was 1.5 cm, as shown in Fig. 1(a). Thirty-two optic-fiber light sources were arranged in a space-filling grid pattern with a tip-tip distance of 0.3 cm at illumination plane *Z*=0, as shown in Fig. 1(b). Detection points were arranged in a space-filling grid pattern with a tip-tip distance of 0.15 cm at the detection plane *Z*=1.5 cm, as shown in Fig. 1(c). The discretized voxel size was selected as 0.01×0.01×0.01 cm for synthetic measurements. Two percent of random Gaussian noise is then added to these synthetic measurements.

#### 3.1 Single fluorescent target

In the experiment, fluorescent target of Cy5.5 fluorochrome was investigated. A laser source with light wavelength 675 nm was chosen, and the corresponding emission wavelength was 694 nm. The matching fluid had a reduced scattering coefficient *µ*′* _{s}* of 6 cm

^{-1}and an absorption coefficient

*µ*of 0.3 cm

_{a}^{-1}for both the excitation and emission wavelengths. The optical properties used were consistent with estimates of the bulk optical properties of mouse tissues [20] for different organs or body regions, and they had been demonstrated successfully in phantom experiments [5]. A single small cubic fluorescent target was embedded at the position (

*x*=0 cm,

*y*=0 cm,

*z*=0.75 cm), which was at a depth of 0.75 cm from the illumination surface. The size of the cubic fluorescent target was 0.01×0.01×0.01 cm. The fluorescence distribution

*n*was set to unit 1 in the target and unit 0 in the background.

#### 3.2 Dual fluorescent targets

All the settings were the same as those in the single-target experiment except two targets were investigated. Two small cubic targets with volume size 0.01×0.01×0.01 cm were placed at an edge-to-edge distance varying from 0.3 cm to 0.1 cm. Reconstructions were performed to evaluate the spatial resolution of the algorithm proposed.

## 4. Results

#### 4.1 Single fluorescent target

A centrally located volume of size 1.8×1.2×1.0 cm was selected as the reconstruction region, as shown by the gray region of the imaging chamber in Fig. 1(a). In the reconstruction the lower bound was set to *n ^{L}_{k}*=0 and the upper bound was set to

*n*=100. Tomographic reconstructions were obtained as illustrated in Fig. 2(d) with zero initial guess

^{U}_{k}*n*

^{0}

_{1}=

*n*

^{L}_{1}and regulation parameter

*λ*=10

_{k}^{-9}.

Figures 2(a), 2(b), and 2(c) show the evolution of the mesh refinement during the reconstruction procedure. The algorithm began with a coarse mesh of 8×5×5 hexahedral elements. Two mesh refinements were performed. The final mesh had 2888 elements that were mostly located around the target position. A total of 12,800 unknowns would be needed for a global uniformly refined mesh with the same mesh resolution around the target. A FMT problem with 12,800 unknowns would spend memory of 0.351GB on the storage of the weigh matrix *W* with double precision. But memory of 0.079 GB would be needed in a problem with 2888 unknowns. Obviously, adaptive refinement adopted largely decreases the number of unknowns, which thus reduces computational cost and enhances computational robustness.

In Fig. 2(d), the top 70% of the contour levels of reconstructed fluorescence distribution are shown, and the red cube represents the real target. The target center was reconstructed at the position (*x*=-0.02 cm, *y*=0.01 cm, *z*=0.75 cm), which is close to the real one. The reconstructed target fluorescence distribution *n _{recons}* was 1.1, with relative error (

*RE*=|

*n*-

_{real}*n*|/|

_{recons}*n*|) 10% to the real target fluorescence distribution

_{real}*n*. The relative error is reasonable, because fluorescent targets with different fluorescence distribution and different size may have similar measurements, especially in the presence of noise. It could be illustrated by the 2.0% measurement error

_{real}*RE*=‖

_{meas}*U*

_{m,meas}-

*U*

_{m,cal}‖/‖

*U*

_{m,meas}‖ between the measurements

*U*

_{m,meas}and the theoretically calculated data

*U*

_{m.cal}, which is equal to the Gaussian noise percentage added to the measurement data.

All the results demonstrated were obtained with an approximate time of 140 s. The time is shorter when compared with an approximate 12 min computation time of the algebraic reconstruction technique (ART) [5, 21] with the same mesh resolution and 100 iterations, in which the analytical solution of diffusion equation is used to generate the weight matrix. The fast computation speed obtained can be attributed to two reasons: First, it is the result of the analytical solution of diffusion equation used, which is two orders of magnitude faster than numerical methods such as FEM [15]. Second, it is the result of the adaptive refinement technique used, which largely reduces the number of unknowns [11, 19].

#### 4.2 Dual fluorescent targets

All the reconstruction configurations were the same as those in the single-target reconstruction. Figure 3 shows the tomographic reconstructions for dual targets separated by 0.3 cm, 0.15 cm and 0.1 cm. The centers of the targets were reconstructed accurately with distance further than 0.15 cm, and only one target was reconstructed when the distance was 0.1 cm. Therefore, the maximum spatial resolution is 0.15 cm with the given experimental configurations in the presence of measurement noise. All the results for dual targets reconstruction are also summarized in Table 1.

## 5. Discussion and conclusion

In this study, a novel adaptive mesh based algorithm for fluorescence molecular tomography (FMT) using analytical solution was investigated. The algorithm was demonstrated using synthetic measurements on a parallel-plate FMT system with matching fluid modeled after the system reported in Graves *et al*. [5].

Computational experiment results show that reconstructions of targets embedded in turbid media are accurate and fast using the algorithm proposed. The reconstruction spends an average computation time of 140 s, which is 20% of the time spent using ART. The fast reconstruction speed is the result of the combination of the adaptive mesh refinement technique and the analytical solution of diffusion equation. The improvement in reconstruction speed is very important for FMT applications in biological experiments and clinical research studies. We believe that the combination of these two techniques used in the algorithm proposed will be a basic strategy for the development of real time FMT.

Although the algorithm is only demonstrated on a FMT system with simple boundary geometry, the idea behind it can be used for a FMT system in the presence of complex boundaries because of the advances in analytical solutions such as Kirchhoff approximation [15, 22–25]. Future work will be focusing on the improvements of the algorithm for a FMT system in the presence of complex boundaries and will target reconstruction in physical phantom using the algorithm proposed.

## Acknowledgments

This work is partially supported by the National Nature Science Foundation of China, the Tsinghua-Yue-Yuen Medical Science Foundation, the National Basic Research Program of China, the Hi-Tech Research and Development Program of China, and the Special Research Fund for the Doctoral Program of Higher Education of China.

## References and links

**1. **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]

**2. **E. E. Graves, R. Weissleder, and V. Ntziachristos, “Fluorescence molecular imaging of small animal tumor models,” Curr. Mol. Med. **4**, 419–430 (2004). [CrossRef] [PubMed]

**3. **M. A. O’Leary, D. A. Boas, B. Chance, and A. G. Yodh, “Fluorescence lifetime imaging in turbid media,” Opt. Lett. **21**, 158–160 (1996). [CrossRef] [PubMed]

**4. **V. Ntziachristos and R. Weissleder, “Experimental three-dimensional fluorescence reconstruction of diffuse media by use of a normalized Born approximation,” Opt. Lett. **26**, 893–895 (2001). [CrossRef]

**5. **E. E. Graves, J. Ripoll, R. Weissleder, and V. Ntziachristos, “A submillimeter resolution fluorescence molecular imaging system for small animal imaging,” Med. Phys. **30**, 901–911 (2003). [CrossRef] [PubMed]

**6. **D. Y. Paithankar, A. U. Chen, B. W. Pogue, M. S. Patterson, and E. M. Sevick-Muraca, “Imaging of fluorescent yield and lifetime from multiply scattered light reemitted from random media,” Appl. Opt. **36**, 2260–2272 (1997). [CrossRef] [PubMed]

**7. **R. Roy and E. M. Sevick-Muraca, “Truncated Newton’s optimization scheme for absorption and fluorescence optical tomography: Part I theory and formulation,” Opt. Express **4**, 353–371 (1999). [CrossRef] [PubMed]

**8. **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]

**9. **M. J. Eppstein, D. E. Dougherty, T. L. Troy, and E. M. Sevick-Muraca, “Biomedical optical tomography using dynamic parameterization and Bayesian conditioning on photon migration measurements,” Appl. Opt. **38**, 2138–2150 (1998). [CrossRef]

**10. **M. J. Eppstein, D. J. Hawrysz, A. Godavarty, and E. M. Sevick-Muraca, “Three dimensional near infrared fluorescence tomography with Bayesian methodologies for image reconstruction from sparse and noisy data sets,” Proc. Nat. Acad. Sci. **99**, 9619–9624 (2002). [CrossRef] [PubMed]

**11. **A. Joshi, W. Bangerth, and E. M. Sevick-Muraca, “Adaptive finite element based tomography for fluorescence optical imaging in tissue,” Opt. Express **12**, 5402–5417 (2004). [CrossRef] [PubMed]

**12. **Y. Lv, J. Tian, W. Cong, G. Wang, J. Luo, W. Yang, J. Shi, and H. Li, “A multilevel adaptive finite element algorithm for bioluminescence tomography,” Opt. Express **14**, 8211–8223 (2006). [CrossRef] [PubMed]

**13. **A. Joshi, W. Bangerth, K. Hwang, J. Rasmussen, and E. M. Sevick-Muraca, “Plane wave fluorescence tomography with adaptive finite elements,” Opt. Lett. **31**, 193–195 (2006). [CrossRef] [PubMed]

**14. **A. Cong and G. Wang, “A finite-element-based reconstruction method for 3D fluorescence tomography,” Opt. Express **13**, 9847–9857(2005). [CrossRef] [PubMed]

**15. **J. Ripoll, V. Ntziachrisos, R. Carminati, and M. Nieto-Vesperinas, “The Kirchhoff approximation for diffusive waves,” Phys. Rev. E **64**, 051917 (2001). [CrossRef]

**16. **S. R. Arridge, “Photon-measurement density functions. Part I: Analytical forms,” Appl. Opt. **34**, 7395–7409 (1995). [CrossRef] [PubMed]

**17. **R. H. Byrd, M. E. Hribar, and J. Nocedal, “An interior point algorithm for large scale nonlinear programming,” SIAM J. Optimization **9**, 877–900 (1999). [CrossRef]

**18. **E. E. Graves, J. P. Culver, J. Ripoll, R. Weissleder, and V. Ntziachristos, “Singular-value analysis and optimization of experimental parameters in fluorescence molecular tomography,” J. Opt. Soc. Am. A **21**, 231–241 (2004). [CrossRef]

**19. **W. Bangerth, “Adaptive finite element methods for the identification of distributed parameters in partial differential equations,” Ph.D. thesis (University of Heidelberg, 2002).

**20. **E. E. Graves, A. Petrovsky, R. Weissleder, and V. Ntziachristos, “In vivo time-resolved optical spectroscopy of mice,” presented at the Optical Society of America Biomedical Optical Spectroscopy and Diagnostics Meeting, Miami, Fla. , April 7–10, 2002. [PubMed]

**21. **X. Intes, V. Ntziachristos, J. P. Culver, A. Yodh, and B. Chance, “Projection access order in algebraic reconstruction technique for diffuse optical tomography,” Phys. Med. Bio. **47**, N1–N10 (2002). [CrossRef]

**22. **R. B. Schulz, J. Ripoll, and V. Ntziachristos, “Experimental fluorescence tomography of tissues with noncontact measurements,” IEEE Trans. Med. Imaging **23**, 492–500 (2004). [CrossRef] [PubMed]

**23. **J. Ripoll, M. Nieto-Vesperinas, R. Weissleder, and V. Ntziachristos, “Fast analytical approximation for arbitrary geometries in diffuse optical tomography,” Opt. Lett. **27**, 527–529 (2002). [CrossRef]

**24. **J. Ripoll and V. Ntziachristos, “Iterative boundary method for diffuse optical tomography,” J. Opt. Soc. Am. A **20**, 1103–1110 (2003). [CrossRef]

**25. **J. Ripoll, R. B. Schulz, and V. Ntziachristos, “Free-space propagation of diffuse light: theory and experiments,” Phys Rev Lett. **91**, 103901 (2003). [CrossRef] [PubMed]