## Abstract

With the development of *in-vivo* free-space fluorescence molecular imaging and multi-modality imaging for small animals, there is a need for new reconstruction methods for real animal-shape models with a large dataset. In this paper we are reporting a novel hybrid adaptive finite element algorithm for fluorescence tomography reconstruction, based on a linear scheme. Two different inversion strategies (Conjugate Gradient and Landweber iterations) are separately applied to the first mesh level and the succeeding levels. The new algorithm was validated by numerical simulations of a 3-D mouse atlas, based on the latest free-space setup of fluorescence tomography with 360° geometry projections. The reconstructed results suggest that we are able to achieve high computational efficiency and spatial resolution for models with irregular shape and inhomogeneous optical properties.

©2007 Optical Society of America

## 1. Introduction

In recent years, fluorescence molecular tomography (FMT) has emerged as a promising tool for *in-vivo* small animal imaging because of its ability to resolve 3D spatial distributions of fluorescence probes associated with molecular and cellular functions. Many efforts have been made to develop new probes, photon migration models, imaging systems, and the corresponding reconstruction methods [1, 2].

During the last decade, imaging systems for fluorescence tomography have evolved over several generations from the early fiber-based systems [3, 4] to the non-contact detection system of slab geometry using charge coupled device (CCD) [5]. Owing to the large number of spatial samples, the latter leads to a sub-millimeter resolution. However, both the two kinds of systems require matching fluid, which is inconvenient for experiments. In 2007 a free-space FMT with full 360° geometry projections was reported [6], which could capture the 3-D surface shape of the imaging animal and get rid of the matching fluid [7]. Additionally, it overcame the limited projection angles of the slab geometry while keeping the high sampling of photon field. Ultimately, this imaging approach is fairly flexible to integrate with other imaging modalities such as X-ray CT or MRI, supplying anatomic information as well.

With the improvement of physical systems for photon collection, there is a need for researches on new reconstruction methods for FMT [1, 2]. On the one hand, the high sampling measurements and the real animal-shape geometry lead to large datasets, which require long computational time and big storage capacity. So the computational efficiency of the reconstruction method need to be improved [2]. On the other hand, an inversion method, which is appropriate for one system, may not function properly for different systems.

In the past four years adaptive finite element (AFE) and adaptive meshing techniques have been employed for reconstruction in optical tomography. In 2003 Gu, et al. [8] proposed a reconstruction method using adaptive meshing for 2-D diffuse optical tomography. The results using phantom measurements demonstrated both qualitative and quantitative improvement of optical image reconstruction. For fluorescence tomography, Joshi et al. [9] developed in 2003 an AFE method for reconstructing fluorescent targets in a reflectance cube geometry using hexahedral mesh. Recently Lee and Joshi et al. [10] presented a fully adaptive finite element based reconstruction algorithm using tetrahedral elements, which was applied to the reconstruction of human breast geometry using point source illumination and detection. In addition to the above works based on finite element analysis, Wang and Song et al. [11] applied adaptive meshing technique in FMT of slab geometry, using the analytical solution of diffusion equations. Compared with the globally uniform discretization, adaptive methods have been shown to significantly reduce the data size and improve the computational efficiency. Additionally they could stabilize the solution of the inversion problem while providing the necessary resolution for FMT of small animal [9–11].

The adaptive reconstruction methods were realized by locally refining the particular region based on *a priori* information derived from a previous reconstruction procedure. Considering this aspect, the first reconstruction procedure is different from the succeeding ones, in that it does not have *a priori* information from a previous procedure. However, each of the four papers mentioned [8–11] use the same inversion techniques for all their different discretization levels, without considering the differences between the first discretization level and the succeeding ones.

Furthermore, the previous adaptive reconstruction algorithms [8–11] were demonstrated using regular geometry such as cylinder, hemisphere, cube, or a combination of these geometries with spatially homogenous optical parameters. However, for free-space imaging of small animals, the reconstruction should be based on an irregular shape of the animal being imaged, which leads to a more accurate description of the forward model and reconstruction results than the regular shape [12–14]. Additionally, the influence of tissue optical heterogeneity on fluorescence tomography reconstruction has been studied for years [15–19]. It has been proved that *a prior* knowledge of inhomogeneous optical properties improves the reconstructed image quality and quantification [17–19]. The distributions of optical parameters corresponding to different organs could be obtained by diffusion optical tomography, utilizing the anatomic information from co-registered MRI [20], X-ray CT [21] and ultrasound [22].

Dual meshing strategy was employed in [8], [9] and [10], where both the forward and inverse meshes were adaptively refined. In fluorescence tomography the inversion was treated as a nonlinear optimization problem and solved it by a truncated Gauss-Newton (GN) strategy [9] or GN with Steihaug-Conjugate Gradient [10]. However, it required that the forward model with updated parameters had to be computed at each iteration procedure. So the computational cost and programming difficulty was increased especially when using the adaptive algorithm.

In this paper we report the development of a novel hybrid adaptive finite element algorithm for reconstruction of FMT, based on a linear scheme. Two different inversion strategies (Conjugate Gradient and Landweber iterations) are separately applied to the first mesh level and the succeeding levels, considering their different characteristics. The new algorithm was validated using surface data generated by numerical simulations, based on the latest free-space system of a rotating subject. The algorithm was applied to a 3-D digital mouse model, incorporating the known optical heterogeneity between different tissues. The results suggest that we are able to achieve high computation efficiency and spatial resolution for models with irregular shape and inhomogeneous distributions of optical parameters.

This paper is organized as follows. In section 2, the diffusion approximation of photon propagation and its finite element solution are initially introduced. Then the formulation of the linear scheme and the novel hybrid adaptive finite element algorithm are presented. In section 3, numerical simulations with fluorescence targets embedded in a 3-D model from a mouse atlas are shown. Finally, we discuss the results and draw a number of conclusions in section 4.

## 2. Methodology

#### 2.1 The model of diffusion equations

In the near infrared spectral window, diffusion equation is always used to describe the photon migration in biological tissues. For the continuous wave FMT with point excitation sources, the model for photon propagation is usually presented using the following coupled diffusion equations [9, 10, 17, 23]:

where Φ* _{x,m}* denotes the photon density for excitation (subscript

*x*) and fluorescence light (subscript

*m*). In the linear model employed in this paper

*µ*is the absorption coefficient and

_{ax,am}*D*=1/3 (

_{x,m}*µ*+µ

_{ax,am}*′*) is the diffusion coefficient of the tissue [3, 24]. The absorption and scattering of the excitation light, caused by the fluorescence probes, were assumed to have little effect on the distribution of Φ

_{sx,sm}*[3], as the fluorescence probes often occupy a very small volume compared to the total imaging area. The fluorescent yield*

_{x}*ηµ*(r) is the unknown fluorescence parameter to be reconstructed, where

_{af}*η*is the quantum yield and

*µ*is the absorption coefficient of fluorescence probe to excitation light.

_{af}Whereas the FMT with 360° projections is carried out by rotating the subject while the excitation source and the detection device are stationary [6, 7], for the reconstruction we consider the phantom or animal geometry model to be stationary while the source and the corresponding detected boundary part to be moving. On the right hand side of Eq. (1.1) *r _{sl}*(

*l*=1,2,…,

*L*) represents the different excitation point source positions with respect to the subject with the amplitude Θ

*.*

_{s}To solve these equations, the Robin-type boundary conditions are implemented on the boundary ∂Ω of domain Ω[25]:

where *n*⃗ denotes the outward normal vector to the surface and *q* is a constant which is approximated

as *q*≈(1+*R*)/(1-*R*) with *R*=-1.4399*k*
^{-2}+0.7099*k*
^{-1}+0.6681+0.0636*k*. *k* is the ratio of optical reflective index of the inner tissue to that outside the boundary.

#### 2.2 Finite element discretization

The finite element method has been widely used for numerically solving the diffusion equations [18, 23, 24], especially with arbitrary geometries [18]. Based on the finite element theory [26], Φ_{x}(r) and Φ* _{m}*(

*r*) in Eq. (1) can be equivalent obtained by solving the following weak form equations:

by implementing the test function *ψ*(*r*).

Discretizing the domain Ω with *N _{p}* vertex node and

*N*elements and employing

_{e}*ψ*(r) as the shape function, Φ

_{i}*(r) can be represented by*

_{x,m}where Φ* _{xi,mi}* denotes the nodal value of Φ

*(r)on the vertex*

_{x,m}*V*.

_{i}Let $\{{\gamma}_{1},{\gamma}_{2},\cdots ,{\gamma}_{{N}_{p}}\}$ be the basis functions and the unknown fluorescent yield *ηµ _{af}*(r) be denoted as

*x*(r). Then it can be approximately expressed as

where *γ _{i}* could be chosen to be the same as

*ψ*.

_{i}By incorporating Eq. (4) and Eq. (5) into Eq. (3), the following matrix equations are obtained

with

$\{\begin{array}{l}{K}_{\mathrm{xi},\mathrm{xj}}={\int}_{\Omega}({D}_{x}\nabla {\psi}_{\mathrm{xi}}\xb7\nabla {\psi}_{\mathrm{xj}}+{\mu}_{\mathrm{ax}}{\psi}_{\mathrm{xi}}{\psi}_{\mathrm{xj}})dr+\frac{1}{2q}{\int}_{\partial \Omega}{\psi}_{\mathrm{xi}}{\psi}_{\mathrm{xj}}dr\\ {L}_{\mathrm{xi},\mathrm{xj}}={\int}_{\Omega}{\Theta}_{s}\delta \left(r-{r}_{s}\right){\psi}_{\mathrm{xj}}dr\end{array}$

and

$\{\begin{array}{l}{K}_{\mathrm{mi},\mathrm{mj}}={\int}_{\Omega}({D}_{m}\nabla {\psi}_{\mathrm{mi}}\xb7\nabla {\psi}_{\mathrm{mj}}+{\mu}_{\mathrm{am}}{\psi}_{\mathrm{mi}}{\psi}_{\mathrm{mj}})dr+\frac{1}{2q}{\int}_{\partial \Omega}{\psi}_{\mathrm{mi}}{\psi}_{\mathrm{mj}}dr\\ {f}_{i,j}={\int}_{\Omega}{\u0424}_{x}\left(\mathrm{r}\right){\psi}_{\mathrm{mi}}{\psi}_{\mathrm{mj}}dr\end{array}.$

#### 2.3 Generation of the linear scheme

Assuming that the optical properties of both the excitation and the emission are known, a linear relationship between the measured fluorescence data on the boundary and the unknown internal fluorescent yield can be generated by the following steps:

Step1. For each excited source position at *r _{sl}*(

*l*=1,2,…, L) with respect to the subject, the corresponding Φ

*is directly obtained by solving Eq. (6).*

_{x,sl}Step2. As the stiff matrix K* _{m}* is symmetrical and positive definite, Eq. (7) can be changed to

for a specific excitation light distribution Φ* _{x,sl}*.

Step3. {Φ* _{m,sl}*} can be divided into two parts {Φ

^{Meas}*} and {Φ*

_{m, sl}

^{NonM}*} where the former consists of the vertices on the surface for measurement and the latter is the remainder of the vertices. Removing the corresponding rows as {Φ*

_{m,sl}

^{NonM}*} from [*

_{m, sl}*B*], the following matrix equation is formed as follows:

_{sl}Step4. Assembling Eq. (9) for different excitation cases, the final weighted matrix for a particular reconstruction mesh system is generated as:

### 2.4 The hybrid adaptive finite element (AFE) reconstruction algorithm

### 2.4.1 Background

In adaptive finite-element based reconstruction technique [27], the imaging domain Ω is not divided into a fixed, uniformly fine mesh so as to achieve the necessary resolution. Instead, Ω is dynamically discretized in several levels from a coarser mesh to a finer mesh with local refinement of the region of interest. A typical AFE based reconstruction algorithm [8–10] includes the following three steps: (1) The unknown parameter is reconstructed firstly utilizing a uniformly coarse mesh, where a rough estimation of the parameter is obtained; (2) Based on the previous reconstructed result, a particular portion of the imaging domain is selected and refined, which lead to a new mesh with the inherited initial value of the required parameter; (3) A new reconstruction procedure is carried out on the refined mesh. If the quality of the reconstructed image is not satisfied, a new iteration of mesh refinement and reconstruction will be repeated.

It is obvious that the first reconstructed procedure is different from the succeeding ones in two aspects: (1) It is based on a uniformly mesh while others are with a locally refined mesh. (2) It has no *a priori* information for the inversion problem while the others can inherit an initial value of the unknown parameter from a previous reconstruction procedure. Considering the above difference between the first reconstructed procedure and the succeeding ones, we proposed a hybrid adaptive finite element reconstruction algorithm.

### 2.4.2 The hybrid AFE reconstruction algorithm

Supposing there is a sequence of tetrahedral mesh levels {Θ_{1},…Θ* _{k}*,…} in the AFE based reconstruction algorithm, a linear relationship as (10) could be generated as

for each mesh level. In fluorescence tomography these are usually ill-posed problems; That is, the weight matrix A* _{k}* is ill-conditioned and the boundary measurements contains noise. The common used Tikhonov regularization is always employed to find the regularized solution, as the minimizer of the following weighted combination of the residual norm and the side constraint:

Usually iterative regularization methods are used for image reconstruction. In the hybrid adaptive finite element algorithm, two different iterative inversion techniques are separately applied to the first mesh level and the succeeding levels. The procedure of the proposed algorithm is illustrated in Fig. 1.

Firstly the imaging domain is discretized into a uniformly coarse mesh Θ_{1}, where there is no *a priori* knowledge about the distribution of fluorescence targets and all the vertices are included in the inversion. As an effective iterative regularization method, the CG algorithm is employed in the first reconstruction procedure. It is applied to the normal equation *A*
^{T}
_{1}
*A*
_{1}
*X*
_{1}=*A*
^{T}
_{1}Φ* ^{mea}*, where the low-frequency components of the solution tend to converge faster than the high-frequency components. Hence, it has some inherent regularization effect where the number of iterations plays the role of the regularization parameter [28]. The initial value of

*X*

_{1}is set to zero and the optimized iteration number is determined by L-curve criterion [29] or by experience. The CG method converges very fast and is able to find an approximate value close to the real distribution in only a few iteration steps.

After the iteration in the first level stops, mesh refinement is applied. According to *a posteriori* error estimation of maximum selection, the elements with greater reconstructed value are selected to be refined [27]. In the implementation, the elements with the average value of the four vertices that is no less than 60% of the maximum value are selected for refinement each time. The boundary mesh of large value is also selected to be refined. Unlike other previous reports we use the longest refinement method to divide the tetrahedral element to second generation elements as showed in Fig. 2. That is, only the longest element edge is refined.

When switching to a finer mesh, the initial value *X*
^{0}
_{k+1} of the (*k*+1)*th* (*k*>=1) mesh level inherits the final solution *X _{k}* of the

*k th*level by linear interpolation as follows:

Then, the optimization in (*k*+1)*th* mesh level is performed. In the proposed algorithm, the reconstructed result in the previous mesh level not only guides mesh refinement and provides an initial value for the refined mesh, but also determines a permissible region for the location of the fluorescence probes. In our study we choose nodes with the top 60% values in 0 k+1 *X* as the permissible region for the location of fluorescence targets. That is, the nodes values outside the permissible region are each set to zero. Thus, columns in *A*+1(*k*>=1)corresponding to the nodes outside the permissible region are removed. The matrix equation in (*k*+1)*th* (*k*>=1)mesh level becomes *A ^{per}_{k+1}X^{per}_{k+1}*=Φ

^{mea}_{k+1}.

For the (*k*+1)*th*(*k*>=1) mesh level where *a priori* knowledge of initial value and the permissible region are supplied by the previous reconstructed result, the following Landweber iterative regularization method is employed as [30]

where the superscripts of ^{X}
^{n+1}
*k*
^{+1} and *Xnk*
_{+1} denote the iteration numbers of Landweber method and *α*=1/*λ*
_{max} with *λ*
_{max}, the maximum eigenvalue of *A ^{per}k*

_{+1}·(

*A*

^{per}_{k+1})

*. The terminating iteration number of Landweber method is determined by both the discrepancy principle and the known value range of fluorescent yield [31]. Usually Landweber method converges slowly. However, as a permissible region for the location of fluorescence targets is determined in the mesh levels other than the first one and the size of the corresponding weighted matrix is reduced, the convergence process is accelerated in this situation. Additionally, the reconstruction value is also stabilized, due to the small step length between neighboring iterations.*

^{T}## 3. Simulations and results

A series of computerized simulations were designed to test the new algorithm. Most parts of the algorithm were coded in Matlab for flexibility. A commercial finite element software COMSOL Multiphysics^{™} was employed to generate and solve the forward problem. Reconstructions were carried out on a personal computer with 2.8 GHz Pentium 4 processor and 1.5 GB RAM.

#### 3.1 Experimental setup

The synthetic measurements were generated based on an imaging configuration of the free-space FMT with 360° geometry projections [6, 7]. The sketch of the system structure is shown in Fig. 3(a). The mouse being imaged is suspended on a rotating stage, with a laser beam focused on its surface. A highly sensitive CCD camera is placed on the opposite side of the laser, collecting the fluorescent photon fields propagating through the mouse.

A 3-D mouse atlas of CT and cryosection data was employed to provide 3-D anatomical information [32, 33] (http://neuroimage.usc.edu/Digimouse.html). Firstly the mouse atlas was imported to COMSOL Multiphysics™ to form the 3-D geometry model. In our simulations, the abdomen part of the mouse atlas with a height of 15*mm* was chosen as the region to be investigated. The rotational axis of the mouse was defined as the z axis with the bottom plane set as *z*=0. Different views of the geometrical model are shown in Figs. 3(c) and 3(d). One unit in the model equals 0.1 *mm*. To simplify the problem, we considered the optical property outside the kidneys to be homogenous. Optical parameters were assigned as *µ _{a}*=0.12

*mm*

^{-1}and µ

*′*=1.2

_{s}*mm*

^{-1}inside the kidneys (the blue part in Figs. 3(c) and 3(d)) and

*µ*=0.23

_{a}*mm*

^{-1}and µ

*′*=1.0

_{s}*mm*

^{-1}outside the kidneys (the green part in Figs. 3(c) and 3(d)) [34].

In this simulation fluorescence tomography of 360° full view was performed using 15 projections with 150° field of view detected for each projection [35]. The collimated laser beam was modeled as an isotropic point source, located one mean free path of photon transport (1/µ*′ _{s}*) beneath the surface [25]. For each projection there were 2 point source positions of 18° interval at the

*z*=7.5

*mm*plane, as shown in Fig. 3(b). Therefore, two images were collected for each projection.

To generate the synthetic measurements on the surface, the mouse geometry model was discretized into a fine tetrahedral-element mesh with the maximum element size of 0.57*mm*
^{3}.Then finite element method was applied to solve Eq. (1) and calculate the photon density distribution. Figure 4 shows an example of the geometry and the mesh when a fluorescence probe was placed in the left kidney. The geometry of the mouse abdomen part was discretized into 12832 nodes and 66047 tetrahedral elements. Figure 4(c) is the boundary view of the forward mesh and Fig. 4(d) is a different view for the internal part with *z*<=8 *mm*. To simulate the real case, one percent of random Gaussian noise was then added to the calculated photon density distributions before applying the inverse algorithm.

#### 3.2 Reconstruction for a single target

At first reconstruction for a single fluorescent target was attempted. To evaluate the performance of the new algorithm in a heterogeneous medium, the fluorescent target was designed to be located at two positions with different optical properties.

In the first case, a cylinder target with a 0.7 *mm* diameter and a 0.6 *mm* height was embedded in the left kidney. The fluorescent yield was set to 0.5. The reconstruction was carried out using the proposed algorithm. The maximum mesh level was set to *K*
_{max}=3. The geometrical model in Fig. 3(c) and Fig. 3(d) was initially discretized into 2263 nodes and 10715 elements, as shown in Figs. 5(a) and 5(b). There was a total of 8606 nodes of measurement in the first mesh level for all the 30 different excitation source positions. Six conjugate gradient (CG) iterations were performed in the first mesh level. The top 30% of the calculated node values were kept as the reconstructed result.

The reconstructed result in the first mesh level is shown in Figs. 6(a) and 6(b) using different views. Figure 6 also shows the evolution of the reconstruction in the second and third mesh levels. The computation time for the reconstruction was approximately 135s in the 1^{st} mesh level, 173s in 2^{nd} mesh level and 197s in 3^{rd} mesh level.

Figure 5(c) illustrates the Residual Error (RE) of ‖*A _{k}X_{k}*-Φ

*‖ for the iterations in each mesh level. The residual error during the computation of the first mesh level descends fastest, when using the conjugate gradient algorithm. In the succeeding mesh levels the RE decreases slowly and smoothly, indicating that the step length between the two neighbor Landweber iterations becomes much smaller. In Fig. 5(c) the first RE value in the second mesh level and the third mesh level is much larger than the rest values in the same level because of the error brought by the linear interpolation from the reconstructed result of the previous level.*

^{meak}Reconstruction for the same fluorescent target embedded in tissues outside the kidney was also performed. All the experimental configurations were the same as those for single target in the left kidney. The results are shown as Fig. 7. Table 1 summarizes the reconstruction results for a single fluorescent target in the above two cases. Additionally, we have also added 5% Gaussian noise to the calculated photon density. The reconstruction results obtained using the proposed algorithm were very similar to those with 1% Gaussian noise except the 4% decrease in the reconstructed value of fluorescent yield.

#### 3.3 Reconstruction of double targets

The numerical experiment with double fluorescent targets was also performed to further demonstrate the validity of the algorithm. Each target was a cylinder of 0.7 *mm* in diameter and 0.6 *mm* in height with the fluorescent yield of 0.5. Both of the targets were embedded outside the kidney region with an edge-to-edge distance of 2 *mm*. Figure 8 shows the reconstruction results in each mesh level. In the third level the two targets could be resolved clearly although their positions are nearer to the surface than that of the real targets. The reconstruction result is also summarized in Table 2.

#### 3.4 Comparison with a fixed mesh

In fluorescence tomography reconstruction based on domain discretization, the position and shape of the reconstructed fluorescence targets as well as the achievable resolution are related to the scale of discretization. For example, for the reconstructed result based on the initial coarse mesh, the inversed fluorescence targets are much larger than the real one, as shown in Fig. 6(a), Fig. 6(b), Fig. 7(a) and Fig. 7(b). However, the reconstructed results utilizing the proposed method show significantly improvement on the reconstructed position and shape, as in Fig. 6(e), Fig. 6(f), Fig. 7(e) and Fig. 7(f). Also the quantification of the reconstruction is improved. When the spatial resolution is concerned, reconstruction under a coarse mesh could not distinguish between the two targets (Fig. 8(a) and Fig. 8(b)), while the final results using the proposed algorithm could (Fig. 8(e) and Fig. 8(f)).

In order to evaluate the proposed hybrid adaptive finite element reconstruction algorithm, reconstruction using conventional finite element-based method was also performed. Figure 9 shows the reconstructed results utilizing the CG method in a fixed fine mesh, with the same fluorescence targets located in the same positions as in Fig. 6 and Fig. 8. For the single fluorescence target, it can be seen that reconstruction using the hybrid AFE method could achieve a more accurate shape and position than when using a fixed mesh. For double targets, the reconstruction using the proposed algorithm can clearly resolve the two targets, while the reconstructed results in a fixed finer mesh can not tell the two targets apart totally and moreover, the reconstructed positions are more inaccurate than that shown in Fig. 8(f). The corresponding time costs and data sizes are shown in Table 3, where the proposed algorithm shows significant decrease in both the time cost and the number of reconstructed nodes. Additionally, we can observe in the last column of Table 3 that the fluorescence yield value reconstructed by the regular finite element method is only about 7% of that by the proposed algorithm.

## 5. Discussion and conclusion

We presented a novel hybrid adaptive finite element algorithm for reconstruction in fluorescence tomography and validated it through numerical simulations. The proposed algorithm was evaluated with a 3-D mouse atlas model of inhomogeneous optical properties using synthetic surface measurements from a non-contact FMT system utilizing 360° geometrical projections. Firstly, reconstruction for a single fluorescent target located in the kidneys or outside the kidneys was studied. For both cases, the difference between the reconstructed target center and the real one was less than 0.5 *mm*. Two fluorescent targets with an edge-to-edge distance of 2 *mm* were also reconstructed. The results show that the two targets were resolved clearly. Compared with the regular finite element-based method using a fixed mesh, the proposed algorithm can reconstruct the shapes, the positions and the yield values more accurately. For reconstruction of double targets, our algorithm shows an improved spatial resolution compared with that of the regular method using a fixed mesh. Moreover, the reconstructed image quality and the spatial resolution could be further improved by optimizing a number of the system setup parameters such as the projection numbers, distribution of excitation sources and the detection range [35].

The improved computational efficiency of the suggested algorithm was also demonstrated, by comparing the dataset size and time cost with a conventional finite element method. It can be seen in Table 3 that the time cost of the proposed algorithm is only 61.5% of that used by the reconstruction in the fixed finer mesh. The number of reconstruction nodes is also reduced significantly, which should lead to a corresponding decrease in memory cost. The employed adaptive mesh strategy makes the reconstruction more intelligent, providing finer mesh around the targets’ location and coarser mesh in other regions. Therefore, compared with a uniformly fine mesh, the total number of reconstructed nodes is significantly reduced, accompanied by improved performance of reconstructed quality and spatial resolution. Apart from the superiority in dataset size, the reconstruction speed of the new algorithm is improved further by the linear scheme employed. For each mesh level our reconstruction algorithm establishes a steady linear relationship between surface measurements and the unknown fluorescent distribution. By comparison, in the nonlinear optimization method, the forward model needs to be re-computed with updated parameters for each iteration procedure. For example, a time cost of 241.2s, 993.6s and 3261.6s, is reported in Ref. [10]. This result was based on the forward/inverse meshes of 1955/554, 6279/794, 12660/945 with 3456 measurements, using a faster PC (3.6 GHz Pentium 4 processor with 3.25 GB RAM) than the one used in our project. This provides additional evidence to the significant improvement that we achieved in computational efficiency and thus, the superiority of our algorithm.

In the hybrid adaptive finite element algorithm, two different inversion techniques are separately applied to the first mesh level and the succeeding levels. In the initial mesh level, all the vertices in the imaging domain are included in the inversion because there is no *a priori* knowledge. In this case, the conjugate gradient (CG) method converges very fast and is able to find an approximate value close to the real distribution of fluorescence yield in only a few iteration steps. Differently, the reconstructions in the succeeding mesh levels inherit an initial value from the previous reconstructed result. Based on *a posteriori* error estimation of maximum selection, a small permissible region where the fluorescent targets are located can be determined. In this study, Landweber iteration algorithm was employed for the inversion of the mesh levels other than the initial level. It has a small step length between neighboring iterations, which could lead to improved reconstruction [36, 37]. Simulations in this manuscript also indicate that it does stabilize the inversion and achieves an estimated image very near to the real distribution. In addition, the permissible region reduces the dataset size. Thus the computing time for Landweber iterations remains short. The residual error for iterations at all three levels is shown in Fig. 5(c). This demonstrates that the RE value drops very quickly in the first level while it decreases much more slowly in the other two levels. The simulation results suggest that the combination of the two methods is efficient and results in a good reconstructed image quality.

With the development of *in-vivo* free-space fluorescence molecular imaging for small animals, two major issues exist: there is a need for reconstruction using the real animal-shape model and there is a need to take into consideration the heterogeneous optical properties of tissue. Several reports considered the importance of the heterogeneous nature of optical properties for the reconstruction for fluorescence tomography [17–19]. A recent evaluation of the influence of optical property maps on inversion of the fluorescence data from mice [19] suggests that heterogeneous optical properties provide significant improvements over homogeneous assumptions for highly absorbing regions (5x background) of sufficient size. Although a normalized scheme has been shown to reduce the influence of heterogeneous background [16], the excitation light needs to be detected in the physical imaging system and the quantitative accuracy of the recovered fluorescence yield is yet to be characterized. The algorithm described here addresses these two issues. It is based on a geometrical model generated from a mouse atlas, which is more accurate in describing the forward model and reconstruction result than previous models using regular shapes [12–14, 18]. In addition, the algorithm employs the anatomic information of internal organs, applying different optical properties to different organs and tissues. Therefore we have been able to simulate the photon propagation in small animals much more closely to the real situation.

The reconstruction errors in fluorescent yield value are mainly caused by the ill-posed nature of FMT, the mesh refinement scale, the effectiveness of the optimization method, the existence of noise, and the absence of a priori knowledge [38]. Nevertheless, while the new algorithm underestimates the values of the fluorescence yield for both single and double targets (Table 3), the reconstructed value using our new algorithm improves the computed yield by more than an order of magnitude when compared to the computations using a fixed mesh. With the known shape and location of fluorescent targets obtained by MRI, Srinivasan et al. [18] reconstructed the value of fluorescent yield accurately using numerical experiments. However, in practical applications, accurate information about location and shape is difficult to obtain. Therefore, further research is required to improve the reconstruction.

In summary, using numerical simulations, we have demonstrated the feasibility and efficiency of a novel hybrid adaptive finite element algorithm, suggesting improved reconstruction of target location and shape while simultaneously reducing the size of the data set and the accompanying computation time. We are planning to extend our proposed method to the whole mouse and test the performance of the algorithm using phantom and mouse experiments.

## Acknowledgments

This work is supported in part 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. **V. Ntziachristos, “Fluorescence molecular imaging,” Annu. Rev. Biomed. Eng. **8**, 1–33 (2006). [CrossRef] [PubMed]

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

**4. **J. Chang, H. L. Graber, and R. L. Barbour, “Luminescence optical tomography of dense scattering media,” JOSA A **14**, 288–299 (1997). [CrossRef] [PubMed]

**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. **N. Deliolanis, T. Lasser, D. Hyde, A. Soubret, J. Ripoll, and V. Ntziachristos, “Free-space fluorescence molecular tomography utilizing 360° geometry projections,” Opt. Lett. **32**, 382–384 (2007). [CrossRef] [PubMed]

**7. **H. Meyer, A. Garofalakis, G. Zacharakis, S. Psycharakis, C. Mamalaki, D. Kioussis, E. N. Economou, V. Ntziachristos, and J. Ripoll, “Noncontact optical imaging in mice with full angular coverage and automatic surface extraction,” Appl. Opt. **46**, 3617–3627 (2007). [CrossRef] [PubMed]

**8. **X. Gu, Y. Xu, and H. Jiang, “Mesh-based enhancement schemes in diffuse optical tomography,” Med. Phys. **30**, 861–869 (2003). [CrossRef] [PubMed]

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

**10. **J. H. Lee, A. Joshi, and E. M. Sevick-Muraca, “Fully adaptive finite element based tomography using tetrahedral dual-meshing for fluorescence enhanced optical imaging in tissue,” Opt. Express **15**, 6955–6975 (2007). [CrossRef] [PubMed]

**11. **D. Wang, X. Song, and J. Bai, “A novel adaptive mesh based algorithm for fluorescence molecular tomography using analytical solution,” Opt. Express **15**, 9722–9730 (2007). [CrossRef] [PubMed]

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

**13. **R. B. Schulz, J. Ripoll, and V. Ntziachristos, “Noncontact optical tomography of turbid media,” Opt. Lett. **28**, 1701–1703 (2003). [CrossRef] [PubMed]

**14. **D. Hyde, A. Soubret, J. Dunham, T. Lasser, E. Miller, D. Brooks, and V. Ntziachristos, “Analysis of reconstructions in full view fluorescence molecular tomography,” Proc. SPIE **6498**, 649803 (2007). [CrossRef]

**15. **X. Li, B. Chance, and A. G. Yodh, “Fluorescent heterogeneities in turbid media: limits for detection, characterization, and comparison with absorption,” Appl. Opt. **37**, 6833–6844 (1998). [CrossRef]

**16. **A. Soubret, J. Ripoll, and V. Ntziachristos, “Accuracy of fluorescent tomography in the presence of heterogeneities:study of the normalized born ratio,” IEEE Trans. Med. Imaging **24**, 1377–1386 (2005). [CrossRef] [PubMed]

**17. **A. B. Milstein, S. Oh, K. J. Webb, C. A. Bouman, Q. Zhang, D. A. Boas, and R. P. Millane, “Fluorescence optical diffusion tomography,” Appl. Opt. **42**, 3081–3094 (2003). [CrossRef] [PubMed]

**18. **S. Srinivasan, B. W. Pogue, S. Davis, and F. Leblond, “Improved quantification of fluorescence in 3-D in a realistic mouse phantom,” Proc. SPIE **6434**, 64340S (2007). [CrossRef]

**19. **S. Bjoern, S. V. Patwardhan, and J. P. Culver, “The influence of Heterogeneous optical properties upon fluorescence diffusion Tomography of small animals,” Springer Proc. in Physics **114**, 361–365 (2007). [CrossRef]

**20. **B. Brooksby, B.W. Pogue, S. Jiang, H. Dehghani, S. Srinivasan, C. Kogel, J. Weaver, S.P. Poplack, and K. D. Paulsen, “Imaging breast adipose and fibroglandular tissue molecular signatures using hybrid MRI-guided near-infrared spectral tomography,” Proceedings of the Natl. Acad. Sci. **103**, 8828–8833 (2006). [CrossRef]

**21. **Q. Zhang, T. J. Brukilacchio, A. Li, J. J. Stott, T. Chaves, E. Hillman, T. Wu, M. Chorlton, E. Rafferty, R. H. Moore, D. B. Kopans, and D. A. Boas, “Coregistered tomographic x-ray and optical breast imaging: initial results,” J. Biomed. Opt. **10**, 024033–0240339 (2005). [CrossRef] [PubMed]

**22. **Q. Zhu, E. B. Cronin, A. A. Currier, H. S. Vine, M. Huang, N. Chen, and C. Xu, “Benign versus malignant breast masses: optical differentiation with US-guided optical imaging reconstruction,” Radiology **237**, 57–66 (2005). [CrossRef] [PubMed]

**23. **H. Jiang, “Frequency-domain fluorescent diffusion tomography: a finite element based algorithm and simulations,” Appl. Opt. **37**, 5337–5343 (1998). [CrossRef]

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

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

**26. **S. S. Rao, “The finite element method in engineering,” (Butterworth-Heinemann, Boston, 1999).

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

**28. **M. Hanke and P. C. Hansen, “Regularization methods for large-scale problems,” Surv. Math. Ind. **3**, 253–315 (1993).

**29. **P. C. Hansen, “Analysis of Discrete ill-posed problems by means of the L-curve,” SIAM Review **34**, 561–580 (1992). [CrossRef]

**30. **L. Landweber, “An iteration formula for Fredholm integaral equations of the first kind,” Am. J. Math. **73**, 615–624 (1951). [CrossRef]

**31. **G. A. Latham, “Best L^{2} Tikhonov Analogue for Landweber Iteration,” Inverse Probl. **14**, 1527–1537 (1998) [CrossRef]

**32. **B. Dogdas, D. Stout, A. Chatziioannou, and R. M. Leahy, “Digimouse: A 3D Whole Body Mouse Atlas from CT and Cryosection Data,” Phys. Med. Biol. **52**, 577–587 (2007). [CrossRef] [PubMed]

**33. **D. Stout, P. Chow, R. Silverman, R. M. Leahy, X. Lewis, S. Gambhir, and A. Chatziioannou, “Creating a whole body digital mouse atlas with PET, CT and cryosection images,” Molecular Imaging and Biology. **4**, S27 (2002).

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

**35. **T. Lasser and V. Ntziachristos, “Optimization of 360° projection fluorescence molecular tomography,” Med. Image Anal. **11**, 389–399 (2007). [CrossRef] [PubMed]

**36. **W. Q. Yang, D. M. Spink, T. A. York, and H. McCann, “An image-reconstruction algorithm based on Landweber’s iteration method for electrical-capacitance tomography,” Meas. Sci. Technol. **10**, 1065–1069 (1999). [CrossRef]

**37. **L. H. Peng, G. Lu, and W. Q. Yang, “Image reconstruction algorithms for electrical capacitance tomography: state of the art,” J. Tsinghua Univ. (Sci & Tech) **44**, 478–484 (2004).

**38. **S. C. Davis, B. W. Pogue, H. Dehghani, and K. D. Paulsen, “Contrast-detail analysis characterizing diffuse optical fluorescence tomography image reconstruction,” J. Biomed. Opt. **10**, 050501-1:3(2005). [CrossRef]