## Abstract

As a new mode of molecular imaging, bioluminescence tomography (BLT) has become a hot topic over the past two years. In this paper, a multilevel adaptive finite element algorithm is developed for BLT reconstruction. In this algorithm, the mesh is adaptively refined according to *a posteriori* error estimation, which helps not only to improve localization and quantification of sources but also to enhance the robustness and efficiency of reconstruction. In the numerical simulation, bioluminescent signals on the body surface of a heterogeneous phantom are synthesized in a molecular optical simulation environment (MOSE) that we developed to model the photon transportation via Monte Carlo simulation. The performance of the algorithm is evaluated in numerical tests involving single and multiple sources in various arrangements. The results demonstrate the merits and potential of the multilevel adaptive approach for BLT.

©2006 Optical Society of America

## 1. Introduction

Molecular imaging, especially small-animal molecular imaging, is a rapidly developing bio-medical imaging field [1]. The goal of molecular imaging is to depict non-invasive *in vivo* cellular and molecular processes sensitively and specifically, such as monitoring multiple molecular events, cell trafficking and targeting. It is well recognized that molecular imaging may be instrumental for tumorigenesis studies, cancer diagnosis, metastasis detection, drug discovery and development, and gene therapies [2, 3, 4]. Due to the high performance and low cost of the photonics-based imaging modalities, optical molecular imaging has attracted much attention, which includes bioluminescence tomography (BLT) [5, 6] and fluorescence molecular tomography (FMT) [7].

Bioluminescence has been extensively applied in biology for decades [8], but bioluminescence tomography (BLT) was only recently developed to reconstruct an underlying source distribution in a small animal such as a mouse from external optical measurement. Traditionally, planar bioluminescent imaging detects only surface light signals and cannot generate a depth location [9]. Although the high detection ability of bioluminescence imaging can’t be provided with fluorescence imaging due to external exciting light sources, BLT is more ill-posed than FMT in theory [6, 10]. Generally, a theoretical study shows that the BLT solution is not unique without sufficient *a priori* information [6]. Earlier BLT algorithms were developed for a single spectral band [11, 12]. Then, based on the wavelength dependence of the optical properties of the tissues and the spectral characteristics of an underlying bioluminescent source, hyper- and multi-spectral BLT methods were also developed and obtained promising results in phantom experiments and mouse studies with single or multiple sources [13, 14, 15]. Furthermore, in the optical-PET (OPET) system development, Alexandrakis *et al*. [15] demonstrated the feasibility and limitation of BLT and underlined that the assumption of a homogeneous optical background was inadequate for BLT reconstruction. At the same time, the importance of a permissible source region in BLT reconstruction was well recognized [16].

For performance evaluation of BLT algorithms, the simulation of photon transportation in the biologic tissues plays an important role [17]. The analytic approach is suitable for the simple homogeneous cases. The numerical approach is preferable in dealing with the complex geometries. Among various numerical methods, Monte Carlo method is the most precise but the least efficient. To synthesize realistic optical data, Monte Carlo simulation is usually the method of choice. Also, the *inverse crime* is likely to occur when closely correlated computational ingredients are used in the forward solver and the inversion scheme [18]. The strategy to avoid the *inverse crime* is that the forward solver can’t be relevant to the reconstruction method. In view of the finite element analysis, the quality of BLT reconstruction is not only related to the signal-noise-ratio (SNR) of measured data but also to the discretization of the domain. To a large extent, the finer the discretized mesh becomes, the better the reconstruction is. However, the overly detailed mesh may aggravate the ill-posedness of the BLT problem and increase the numerical instability and computational cost. Hence, the optimization of finite element meshing is critical in the BLT reconstruction. In fluorescence tomography, an adaptive finite element method was proposed to solve the aforementioned dilemma using a dual meshing scheme [19]. Although the dual meshing strategy improves the flexibility of the tomographic algorithm, it also increases the computational cost and programming difficulty especially in the adaptive algorithm.

In this paper, a BLT algorithm is developed based on the adaptive finite element methods (FEMs), which utilizes the single initial coarse volumetric mesh to recover the source distribution quantitatively. This diffusion equation model based method also establishes the linear relationship between the measured data and the unknown source distribution through the permissible source region concept. Then, the adaptive mesh refinement is performed via decomposing the selected mesh elements by virtue of two different *a posteriori* error estimation techniques in the forbidden and permissible source regions. In the adaptive procedure, the optimization is done using a modified Newton method coupled with multilevel and active-set strategies. In the numerical simulation, bioluminescent signals on the body surface of a heterogeneous phantom are synthesized in a molecular optical simulation environment (MOSE) that we developed to model photon transportation via Monte Carlo simulation. In the next section, we present the BLT framework based on the multilevel adaptive finite element analysis. In the third section, we evaluate the performance of the proposed algorithm using the MOSE based synthetic data. In the last section, we discuss the relevant issues and conclude the paper.

## 2. Methods

#### 2.1. Diffusion approximation and boundary condition

In a bioluminescence imaging experiment, luciferase can be introduced into various types of cells, organisms, and genes in a living mouse. Then, luciferin is combined with the luciferase in the presence of oxygen and ATP to generate bioluminescent signals of about 600*nm* in wavelength. The bioluminescent photon in the biological tissue is subject to both scattering and absorption. In this wavelength range, photon scattering predominates over photon absorption and the diffusion equation has been widely used in bio-photonics [16, 17, 20]. When the bioluminescence imaging experiment is carried out in a dark environment, the propagation of bioluminescent photons in the biological tissue can be well modeled by the steady-state diffusion equation and Robin boundary condition [16, 21]:

where Ω and ∂Ω are the domain and its boundary respectively; Φ(**x**) denotes the photon flux density [*Watts*/*mm*
^{2}]; *S*(**x**) is the source energy density [*Watts*/*mm*
^{3}]; *μ*_{a}
(**x**) is the absorption coefficient [*mm*
^{-1}]; *D*(**x**)=1/(3(*μ*_{a}
(**x**)+ (1 - *g*)/is(**x**))) is the optical diffusion coefficient [*mm*], *μ*_{s}
(**x**) the scattering coefficient [*mm*
^{-1}], and *g* the anisotropy parameter; **v** is the unit outer normal on ∂Ω. Given the mismatch between the refractive indices *n* for Ω and *n*′ for the external medium, *A*(**x**;*n*,*n*′) can be approximately represented:

where *n*′ is close to 1.0 when the mouse is in air; *R*(**x**) can be approximated by *R*(**x**) ≈ -1.4399*n*
^{-2} + 0.7099*n*
^{-1} + 0.6681 +0.0636*n* [21]. The measured quantity is the outgoing flux density *Q*(**x**) on ∂Ω, that is:

#### 2.2. Reconstruction based on multilevel adaptive FEMs

Based on the finite element theory [22], the weak solution of the flux density Φ(**x**) is given through Eqs. (1-1) and (1-2):

where *H*
^{1}(Ω) is the Sobolev space. Considering the approximation Φ^{h} of the exact solution Φ and the uniqueness solution of the BLT problem, the following error bound was derived [23]:

where Φ(*S*_{λ}
)∈/*H*
^{1}(Ω) and *S*_{λ}
∈*L*
^{2}(Ω); Φ^{h}(${S}_{\lambda}^{h}$
) and ${S}_{\lambda}^{h}$
corresponding to Φ(*Sλ*) and *Sλ* are obtained by introducing a regular triangulation $\mathcal{T}$_{h}
of Ω and the finite element space $\mathcal{V}$_{h}
;*λ* is the regularization parameter; *h* denotes the maximum element size of the triangulation and *c* is a constant independent of *λ* and *h*. Based on Eq. (5), it is beneficial to the improvement of the BLT solution to reduce the maximum element size *h*.

The adaptive finite element analysis has been studied over about two decades [24], and involves iterative acceleration, *a posteriori* error estimation and mesh refinement. In the framework of adaptive finite element analysis, let {$\mathcal{T}$_{1},… $\mathcal{T}$_{k},…} be a sequence of nested triangulation of the given domain Ω based on the adaptive mesh refinement, where the sequence gradually changes from coarse to fine along with the increase in *k*. The spaces of linear finite elements $\mathcal{V}$_{k} are introduced on the discretized levels $\mathcal{T}$_{k}, satisfying $\mathcal{V}$_{1} ⊂… ⊂ $\mathcal{V}$_{k} ⊂… ⊂ *H*
^{1}(Ω). Now, we only consider the *k*th discretized level which includes *N*
_{$\mathcal{T}$k} elements and *N*
_{$\mathcal{P}$k} vertex nodes. {${\psi}_{1}^{k}$,…, *ψ*
^{k}
_{N$\mathcal{P}$k}} is the nodal basis of the space $\mathcal{V}$_{k}. Φ_{k}(**x**) is an approximation of Φ(**x**) on the *k*th discretized level:

where ${\varphi}_{i}^{k}$
(**x**) is the *i*th node value. Let {${\gamma}_{1}^{k}$, …,*γ ^{k}_{NSk}*} be the interpolation basis functions,

*S*(

**x**) is approximated by

*S*

_{k}(x

**):**

**$${\mathit{S}}_{k}\left(\mathbf{x}\right)=\sum _{i=1}^{{{N}_{S}}_{k}}{s}_{i}^{k}\left(\mathbf{x}\right){\gamma}_{i}^{k}\left(\mathbf{x}\right)$$**

**where N_{Sk} and ${s}_{i}^{k}$
are the number of the interpolation basis functions and the interpolation nodal values respectively. The selection of interpolation basis functions ${\gamma}_{i}^{k}$ may be the same with or different from that of nodal basis functions ${\psi}_{i}^{k}$
, which depends on the choice of source variables ${s}_{i}^{k}$
. Here, we select element as the basis source variable and piecewise constant functions as the interpolation basis functions, which are different from the piecewise linear nodal basis functions ${\psi}_{i}^{k}$
. Incorporating Eqs. (6) and (7) with Eq. (4), we have**

**$$\left({\displaystyle \sum _{i=1}^{{N}_{{P}_{k}}}\left({\displaystyle {\int}_{{\tau}_{i}}\left(D(\nabla {\psi}_{m}^{k})\cdot (\nabla {\psi}_{n}^{k})+{\mu}_{a}{\psi}_{m}^{k}\cdot {\psi}_{n}^{k}\right)\text{d}\mathbf{x}+{\displaystyle {\int}_{\partial {\tau}_{i}}\frac{{\psi}_{m}^{k}\cdot {\psi}_{n}^{k}}{2{A}_{n}}}\text{d}\mathbf{x}}\right)}\right)\cdot {\Phi}_{k}=\left(\sum _{i=1}^{{{N}_{P}}_{k}}{\int}_{{\tau}_{i}}{\psi}_{m}^{k}\bullet {\psi}_{n}^{k}d\mathbf{x}\right)\bullet {S}_{k}\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}$$**

**where τ_{i}
is the ith element; ∂τ_{i}
, the ith boundary element, on which the integration is equal to 0 if τ_{i}
is not on the boundary. Then, the matrix form of Eq. (8) can be obtained as follows:**

**Because M_{k}
is a sparse positive definite matrix, we have Φ_{k} = ${M}_{k}^{-1}$
F_{k}
S_{k}
. The aim of BLT is to reconstruct source distribution S_{k}
from the boundary measured flux data ${\mathrm{\Phi}}_{k}^{m}$, which can be obtained through removing the interior discretized points values from Φ_{k}. In order to improve the BLT reconstruction and reduce the ill-posed nature, a priori knowledge can be obtained through the surface measurement and anatomical and optical information. Therefore, S_{k}
can be divided into the permissible source region ${S}_{k}^{p}$
and the forbidden region ${S}_{k}^{*}$. Deleting those columns of M
^{-1}
F_{k}
corresponding to ${S}_{k}^{*}$ and retaining those rows corresponding to ${\mathrm{\Phi}}_{k}^{m}$ in ${M}_{k}^{-1}$
F_{k}
, a linear relationship between the measurable boundary flux ${\mathrm{\Phi}}_{k}^{m}$ and the unknown source density ${S}_{k}^{p}$
is established:**

**Then, we define the following kth level optimization problem to compute source distribution based on the Tikhonov regularization method:**

**$$\underset{{S}_{\mathit{inf}}^{k}\le {S}_{k}^{p}\le {S}_{\mathit{sup}}^{k}}{min}={\Theta}_{k}\left({S}_{k}^{p}\right)=\left\{{\parallel {A}_{k}{S}_{k}^{p}-{\Phi}_{k}^{m}\parallel}_{\Lambda}+{\lambda}_{k}{\eta}_{k}\left({S}_{k}^{p}\right)\right\}$$**

**where ${S}_{\mathit{\text{inf}}}^{k}$
and ${S}_{\mathit{\text{sup}}}^{k}$
are the kth level lower and upper bounds of the source density; Λ is the weight matrix, ∥V∥_{Λ} = V^{T}
ΛV; λ_{k}
the regularization parameter; and η_{k}
(∙) the penalty function.**

**A modified Newton method with the active set strategy was used to deal with the minimization problem effectively [25]. However, the ill-posedness of the inverse problem results in the sensitivity of the solution to the initial guess ${S}_{k}^{0}$ of ${S}_{k}^{p}$
. Fortunately, our proposed multilevel adaptive BLT algorithm can improve the reconstruction stability and accelerate the convergence speed. First, the optimization procedure is activated using a small initial value ${S}_{1}^{0}$ on the coarsest mesh level. Let the prolongation operator be denoted by ${I}_{k}^{k+1}$ from the level k to k+1 and ${S}_{k}^{r}$
be the optimized results on the kth level, we have**

**where ${I}_{k}^{k+1}$ can be realized through the inheritance of son elements from father one. Due to the adaptive mesh refinement, several types of prolongations are utilized depending on the refinement form of the single element, which will be thoroughly described in the last paragraph of this subsection. Through correctness of the results on the coarse level, the optimization on the fine mesh will become more robust and effective than BLT without such a multilevel meshing.**

**The switch condition from the kth level to the (k+1)th level and the stopping criterion of the program are important. In the multilevel adaptive algorithm, we select the norm of the gradient ∥g
^{i}
_{Θk} ∥ and the distance between the last two steps ∥ d^{i}_{Sk} ∥ as the indexes, where g
^{i}
_{Θk} =∇Θ_{k}:(${S}_{k}^{i}$
) and d^{i}_{Sk} =${S}_{k}^{i}$
-${S}_{k}^{i-1}$. Let ${\epsilon}_{g}^{\left(k\right)}$ and ${\epsilon}_{d}^{\left(k\right)}$ be the corresponding thresholds. When ∥g
^{i}
_{Θk} ∥ is less
than ${\epsilon}_{g}^{\left(k\right)}$ or ∥ d^{i}_{Sk} ∥ is less than ${\epsilon}_{d}^{\left(k\right)}$, the local mesh refinement will be triggered. Note that BLT with a coarsely discretized mesh means less unknown variables, higher computational efficiency and better numerical stability than that with a finely discretized counterpart. Hence the optimization of the objective function Θ_{k}(${S}_{k}^{p}$
) is indispensable on the coarse mesh. As far as the stopping criterion is concerned, we utilize the discrepancy between the measured and computational boundary nodal flux data or the maximum number of mesh refinement L_{max}
to evaluate if the whole reconstruction procedure should be terminated, that is, ∥ ${\mathrm{\Phi}}_{k}^{m}$- Φ^{c}
_{k}∥< ε_{Φ} or k≥ L_{max}
, where ∥${\mathrm{\Phi}}_{k}^{m}$ -${\mathrm{\Phi}}_{k}^{c}$∥ is the error energy norm and ε
_{Φ} denotes the stopping threshold.**

*A posteriori* error estimation plays a significant role in the multilevel adaptive algorithm [26]. To guarantee a quasi-uniform mesh of the whole triangulation *$\mathcal{T}$*_{k}
and increase the measured points on the domain boundary where the flux density Φ_{k} has a sharp gradient, two different error estimates are used in the forbidden and permissible source regions, which are based on the hierarchical defect correction technique [26] and the direct maximum selection method. The permissible source region can be decided through incorporating various types of obtainable *a priori* knowledge, such as optical, biomedical, physiological and anatomical information and so on, which may reduce the non-uniqueness of the BLT problem. Therefore, on a coarse mesh which brings few discretized elements in the permissible source region, the elements with higher values of the optimization results most likely represent actual source locations despite they may look larger and darker than the actual. By selecting these elements for refinement, the BLT solution will be improved in general. In the forbidden region, the quadratic finite element spaces {$\mathcal{W}$_{1},…$\mathcal{W}$_{k},…} are defined corresponding to the linear finite element spaces {$\mathcal{V}$_{1}, …$\mathcal{V}$_{k} …} to perform the hierarchical defect correction based error estimation. On the *k*th mesh level, let *e _{$\mathcal{V}$k}* be an error between an approximated solution $\overline{\Phi}$

_{$\mathcal{V}$k}on the linear finite element space $\mathcal{V}$

_{k}and the real one Φ. After the approximated solution $\overline{\Phi}$ $\mathcal{W}$

_{k}is obtained on the quadratic finite element space $\mathcal{W}$

_{k}, the total error

*e$\mathcal{W}$*

_{k}between Φ and $\overline{\Phi}$

_{$\mathcal{W}$k}is approximately equal to the sum of

*e*and

_{$\mathcal{V}$k}*e*

^{$\mathcal{W}$k}

_{$\mathcal{V}$k}which denotes the error between $\overline{\Phi}$

_{$\mathcal{W}$k}and $\overline{\Phi}$

_{$\mathcal{V}$k}. In the space ℰ

_{k}obtained by the decomposition of $\mathcal{W}$

_{k}into $\mathcal{V}$

_{k}and ℰ

_{k}, i.e. $\mathcal{W}$

_{k}= $\mathcal{V}$

_{k}⊕ ℰ

_{k}, the approximation

*ẽ*

^{$\mathcal{W}$k}

_{$\mathcal{V}$k}to

*e*

^{$\mathcal{W}$k}

_{$\mathcal{V}$k}can be made as:

**$${\tilde{e}}_{{\mathcal{V}}_{k}}^{{\mathcal{W}}_{k}}={D}_{{\mathcal{E}}_{k}}^{-1}r{\mathcal{E}}_{k}$$**

**where D_{ℰk} is the approximation to M_{ℰk}; and r_{ℰk} the iterative residual in the space ℰ_{k}
. With the error indicator ẽ^{$\mathcal{W}$k}
_{$\mathcal{V}$k}, we may select the portion of the elements in the forbidden region for refinement.**

**A reasonable triangulation to the given domain facilitates an implementation of the local mesh refinement. In the biomedical research field, triangular and tetrahedral elements are popular for representing the anatomical surface topology and volumetric features [27, 28, 29]. After the triangulation is specified, the local mesh refinement can be done via the so-called red refinement, which divides a selected tetrahedron into eight son tetrahedral elements [30]. Following the red refinement, the irregular refinement green closure is used to close the triangulation for consistence. In three dimensions, there are 64 possible edge refinement patterns to deal with. Using symmetry consideration, we restrict ourselves to four types of irregular refinements shown in Fig. 1(a)–1(d). The red refinement coupled with the green closure reasonably accomplishes the local mesh refinement.**

**3. Numerical studies**

**A series of computational experiments was designed to evaluate our multilevel adaptive BLT algorithm. A heterogeneous cylindrical phantom of 30 mm height and 10mm radius was set up. It consisted of four ellipsoids and one cylinder to represent muscle (M), lungs (Lu), heart (H), bone (B) and liver (Li), as shown in Fig. 2(a). Optical parameters from the literature [15] and the references therein were assigned to each of the five components, as listed in Table 1.**

**3.1. Photon transportation simulation using MOSE**

**Monte Carlo (MC) method remains a gold standard for photon transportation simulation because of its accuracy and flexibility [31, 32]. Using MC method, we have been developing a molecular optical simulation environment (MOSE), which takes 2D/3D analytical models and microCT/MRI slices to define the object geometry. Using an optimized facet searching method, the simulation speed of our MOSE is faster by an order of magnitude than some popular programs, such as the commercial software TracePro [33]. In addition, MOSE is written in Visual C++ and OpenGL with a user-friendly interface, through which bioluminescent sources, lens and detectors can be flexibly arranged [34]. In the MOSE simulation, the aforementioned heterogeneous phantom was discretized by triangular elements to produce the synthetic data. The mesh of the phantom in Fig. 2(b) consisted of 34072 triangles and 11499 surface measurement points with an average element diameter of about 0.5 mm.**

**3.2. Light source reconstruction**

**First, a single bioluminescent source was embedded in the right lung. In the MOSE simulation, a spherical source with 1.0 mm radius and 0.238nano-Watts/mm
^{3} power density was centered at (-3.0,5.0,15.0), as shown in Fig. 4(a). The source was sampled by 1.0 × 10^{6} photons and was assumed to obey the uniform distribution. In the reconstruction procedure, the initial discretization of the phantom was the coarse volumetric mesh in Fig. 2(c), which was rather different from the mesh used in the forward simulation. The flux density values on the surface points were obtained by the interpolation between the available values obtained using MOSE. We set the lower bound ${S}_{\mathit{\text{inf}}}^{k}$
= 0.0, the upper bound ${S}_{\mathit{\text{sup}}}^{k}$
a sufficiently large positive number, the weighted matrix ∧ the unit matrix, the penalty function η_{k}
(X)=X^{T}X and the regularization parameter λ_{k}
= 1.16 × 10^{-10} at each level. In addition, the gradient tolerance ${\epsilon}_{g}^{\left(k\right)}$ and the distance tolerance ${\epsilon}_{d}^{\left(k\right)}$ were invariable at each level, being equal to 1.0 × 10^{-17} and 1.0 × 10^{-16} respectively, and the stopping threshold ε
_{Φ} and the maximum number of mesh refinement L_{max}
were set to 2.1 × 10^{-9} and 4 respectively. In each adaptive mesh refinement, we refined 50% of the reconstructed elements based on the direct maximum selection method in the permissible source region and decomposed 1.5% elements in the forbidden region according to the error indicator ẽ^{$\mathcal{W}$k}
_{$\mathcal{V}$k}.**

**When the bioluminescent source is placed in the phantom or small animal, the peak intensity attenuation and full width at half maximum (FWHM) of the emitting light varies with the depth change of source [9]. Despite that the quantitative change is affected by optical property parameters and light wavelength, the decision of permissible source region is beneficial from a priori information above. Figure 3 shows four views of the surface light power distribution with an angular increment of 90 degrees, and red lines as the isoline of light power distribution depict the diffusive results of source on the phantom surface. Through the difference of surface light power distribution in four views, we can roughly infer the source position and outline a permissible source region. Furthermore, we may utilize the obtainable a priori knowledge including the anatomical information of the phantom, the corresponding optical parameters and the dependent relationship between FWHM and source depth [9] to specify the domain**

**$${P}_{s}=\{\left(x,y,z\right)\mid 13.0<z<17.0,\left(x,y,z\right)\in \mathit{Right}\phantom{\rule{.2em}{0ex}}\mathit{lung}\}$$**

**as the permissible source region.**

**To improve the computation speed, our multilevel adaptive algorithm was coded in C++. When the initial guess ${S}_{1}^{0}$ was set to 1.0 × 10 ^{-6}
nano
-Watts/mm
^{3}, the reconstruction procedure was terminated after two adaptive iterations at ∥${\mathrm{\Phi}}_{k}^{m}$-${\mathrm{\Phi}}_{k}^{c}$∥ = 2.05 × 10^{-9}, which took about 150 seconds on an desktop computer with Pentium 4 2.8GHz and 1GB RAM. The final reconstruction is shown in Fig. 4(a), which has the maximum source density 0.219nano-Watts/mm
^{3}. The relative error (RE) between the actual source and reconstructed results was 8.0%, which was computed according to RE = ∣S_{recons}
- S_{real}
}∣/S_{real}
}. Although the reconstructed source was slightly closer to the phantom surface than the actual source, its center was located at (-3.20,5.77,15.48) and very close to the actual one. The mesh evolution for the permissible source region and the reconstructed source is shown in Fig. 5.**

**In order to compare the reconstructed results from the fixed mesh and the adaptive evolution mesh, we selected three fixed meshes with an average element diameter of about 4 mm, 2mm and 1mm, respectively. The corresponding reconstructed results on different discretized meshes are shown in Fig. 5(a) and 6. The quantitative comparison between the fixed mesh of different discretized scales and the adaptive mesh are shown in Table 2. Incorporating the permissible source region into the tomographic algorithm, the BLT reconstruction on the coarse and fine meshes could all localize the bioluminescent source, but the reconstructed results were severely affected by the size and distribution of the fixed discretized elements, especially on the coarse mesh. Despite that the position and shape of reconstructed source was improved as the mesh became finer, the time cost also increased correspondingly. In addition, based on the same selection of permissible source region, it is noted that the quantitative information of source density was not remarkably improved as the mesh became finer due to the absence of multilevel strategy.**

**Furthermore, the multilevel adaptive algorithm was employed in the case of multiple light sources. In the experiment, four light sources were placed in the lungs, that is, two sources in the left lung and the other two in the right lung. The power density of each source was the same as that in the single source case. The initial permissible source region was set to contain the lungs between 11.0 mm and 19.0mm along the Z-axis direction. The position of the actual light sources and the final results are shown in Fig. 4(b) and the quantitative information also summarized in Table 3. In the four-source case, the average relative error between the actual and reconstructed source densities were 32.2%, but the refined tetrahedral elements representing the reconstructed sources were around the centers of the actual light sources.**

**3.3. Spatial resolution evaluation**

**To evaluate the spatial resolution of our multilevel adaptive algorithm, we designed a set of experimental models with two spatially close light sources. The two spherical sources of 1.0 mm radius and 0.238nano-Watts/mm
^{3} power density were placed at various edge-to-edge distances from 1.5mm to 0.5mm. The same initial permissible source region was used in all the experiments, which was set to include the right lung between 11.5mm and 18.5mm along the Z-axis. Figure 7 compares the actual and reconstructed light sources separated by 1.5mm, 1.0mm and 0.5mm, respectively. Although the tetrahedral elements of the reconstructed sources were adjacent to each other, a distinct separation was visible for the 1.5mm and 1.0mm separation. The comparative data for the 1.5mm and 1.0mm instances is also summarized in Table 4. The reconstructed source centers and power densities all indicated that the performance of the multilevel adaptive algorithm was very satisfactory.**

**4. Discussions and conclusion**

**We have developed the novel multilevel adaptive finite element algorithm to reconstruct the bioluminescent source distribution, and evaluated its performance in numerical simulation. Our results indicate satisfactory localization and quantification in terms of source center and power, as well as a spatial resolution in the millimeter domain is feasible with aid of a priori information on the permissible source region. Our work demonstrates the merits and potential of the multilevel adaptive approach for optical molecular tomography. Taking into account the whole computational experiment framework, three novel features have been applied to establish and evaluate the multilevel adaptive finite element algorithm.**

**First, the multilevel adaptive strategy makes the BLT reconstruction intelligent. The elements for the reconstructed source become more and more refined through a posteriori error estimation, which generally leads to better location and quantification to the actual source. The consumption time, robustness and stability of the optimization procedure are beneficial from multilevel tactics with the accompaniment of the adaptive mesh evolution. Although the source reconstruction method of BLT based on the uniform refinement mesh has the multilevel merits to probably acquire the preferable results than the proposed algorithm, the reconstruction process will likely suffer from its expensive computation cost and instability which arises from the tremendous dimension of the unknown variables and the severely ill-posed characteristic due to the overly detailed discretization to the given domain, to say nothing of the direct reconstruction method on the overly detailed mesh without possession of the adaptive and multilevel features.**

**Secondly, a priori knowledge, especially the utility of permissible source region, improves the BLT reconstruction significantly. The more a priori information we have, the more precise and stable the BLT reconstruction becomes. In this work, four factors including the light power distribution of the phantom surface, the anatomical structure of the phantom, the peak intensity attenuation of the photon transportation and the FWHM of the bioluminescent source in the different biological tissues help choose the permissible source region. Recently, a visual hulls based method has been proposed to estimate a tumor shape from a series of 2D bioluminescence views, which can be also used to specify the permissible source region [35].**

**Lastly, the theoretical and computational ingredients in the forward process are made different from that in the source inversion to avoid the inverse crime. In view of the precise simulation of Monte Carlo method to photon transportation in the biological tissues, it is more reasonable than other methods to synthesize the measured data. Different meshes are used to further decouple the forward and inverse processes.**

**Although the proposed multilevel adaptive algorithm has been evaluated through a series of computational experiments, many practical factors such as the detection method are not fully considered. Recently, a non-contract detection mode with a highly sensitive charge-coupled device (CCD) camera has been developed [1], which improves the SNR and the spatial sampling of the detected signal and simplifies the experimental equipment and operation. We believe that the utility of the MC method based synthetic measured data insures the appearance of the same performance with the computational experiments if the equipment and calibration techniques are precisely designed in the real BLT experiment. Our future work will focus on the source reconstruction using the multilevel adaptive algorithm and the physical phantoms. Relevant results will be reported later.**

**Acknowledgments**

**This work is supported by NBRPC(2006CB705700), NSFDYS (60225008) and NSFC (30370418, 30500131 and 60532050), China. WX. Cong and G. Wang are supported in part by an NIH/NIBIB grant EB001685.**

**References and links**

**
1
. **
V.
Ntziachristos
,
J.
Ripoll
,
L. V.
Wang
, and
R.
Weisslder
, “
Looking and listening to light: the evolution of whole body photonic imaging
,”
Nature Biotechnology
**
23
**
(3),
313
–
320
(
2005
).
[CrossRef] [PubMed]

**
2
. **
C.
Contag
and
M. H.
Bachmann
, “
Advances in bioluminescence imaging of gene expression
,”
Annu. Rev. Bio-med. Eng.
**
4
**
,
235
–
260
(
2002
).
[CrossRef]

**
3
. **
S.
Bhaumik
and
S. S.
Gambhir
, “
Optical imaging of Renilla luciferase reporter gene expression in living mice
,”
Proc. Natl. Acad. Sci. USA
**
99
**
,
377
–
382
(
2002
).
[CrossRef]

**
4
. **
T. F.
Massoud
and
S. S.
Gambhir
, “
Molecular imaging in living subjects: seeing fundamental biological processes in a new light
,”
Genes Dev.
**
17
**
,
545
–
580
(
2003
).
[CrossRef] [PubMed]

**
5
. **
G.
Wang
,
E. A.
Hoffman
,
G.
McLennan
,
L. V.
Wang
,
M.
Suter
, and
J. F.
Meinel
, “
Development of the first bioluminescence ct scanner
,”
Radiology
**
229(P)
**
,
566
(
2003
).

**
6
. **
G.
Wang
,
Y.
Li
, and
M.
Jiang
, “
Uniqueness theorems in bioluminescence tomography
,”
Med. Phys.
**
31
**
,
2289
–
2299
(
2004
).
[CrossRef] [PubMed]

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

**
8
. **
Thérèse
and
J. W.
Hastings
, “
Bioluminescence
,”
Annu. Rev. Cell Dev. Bi.
**
14
**
,
197
–
230
(
1998
).
[CrossRef]

**
9
. **
W.
Rice
,
M. D.
Cable
, and
M. B.
Nelson
, “
In vivo imaging of light-emitting probes
,”
J. Biomed. Opt.
**
6
**
,
432
–
440
(
2001
).
[CrossRef] [PubMed]

**
10
. **
M.
Jiang
and
G.
Wang
, “
Image reconstruction for bioluminescence tomography
,”
Proc. SPIE
**
5535
**
,
335
–
351
(
2004
).
[CrossRef]

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

**
12
. **
X.
Gu
,
Q.
Zhang
,
L.
Larcom
, and
H.
Jiang
, “
Three-dimensional bioluminescence tomography with model-based reconstruction
,”
Opt. Express
**
12
**
,
3996
–
4000
(
2004
),
http://www.opticsexpress.org/abstract.cfm?URI=OPEX-12-17-3996
.
[CrossRef] [PubMed]

**
13
. **
C.
Kuo
,
O.
Coquoz
,
D. G.
Stearns
, and
B. W.
Rice
, “
Diffuse luminescence tomography of in vivo bioluminescence markers using multi-spetral data
,” in
*
Proceedings of the 3rd International Meeting of the Society
*
(
MIT Press
,
2004
), p.
227
.

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

**
16
. **
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
),
http://www.opticsexpress.org/abstract.cfm?URI=OPEX-13-18-6756
.
[CrossRef] [PubMed]

**
17
. **
A. P.
Gibson
,
J. C.
Hebden
, and
S. R.
Arridge
, “
Recent advances in diffuse optical imaging
,”
Phys. Med. Biol.
**
50
**
,
R1
–
R43
(
2005
).
[CrossRef] [PubMed]

**
18
. **
S.
Holder
,
*
Electrical Impedance Tomography
*
, (
Institute of Physics Publishing, Bristol and Philadelphia
,
2005
).

**
19
. **
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
),
www.opticsexpress.org/abstract.cfm?URI=OPEX-12-22-5402
.
[CrossRef] [PubMed]

**
20
. **
S. R.
Arridge
,
M.
Schweiger
,
M.
Hiraoka
, and
D. T.
Delpy
, “
A finite element approach for modeling photon transport in tissue
,”
Med. Phys.
**
20
**
,
299
–
309
(
1993
).
[CrossRef] [PubMed]

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

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

**
23
. **
G.
Wang
,
M.
Jiang
,
J.
Tian
,
W.
Cong
,
Y.
Li
,
W.
Han
,
D.
Kumar
,
X.
Qian
,
H.
Shen
,
T.
Zhou
,
J.
Cheng
,
Y.
Lv
,
H.
Li
, and
J.
Luo
, “
Recent development in bioluminescence tomography
,” presented in the third IEEE International Symposium on Biomedical Imaging (ISBI 2006),
Virginia, USA
, 6–9 Apr.
2006
.

**
24
. **
W.
Bangerth
and
R.
Rannacher
,
*
Adaptive finite element methods for differential equations
*
, (
Birkhäuser Verlag
,
2003
).

**
25
. **
P. E.
Gill
,
W.
Murray
, and
M.
Wright
,
*
Practical optimization
*
, (
Academic Press, New York
,
1981
).

**
26
. **
M.
Ainsworth
and
J. T.
Oden
,
*
A posteriori error estimation in finite element analysis
*
, (
Wiley
,
2000
).
[CrossRef]

**
27
. **
W. E.
Lorensen
and
H. E.
Cline
, “
Marching cubes: a high resolution 3D surface construction algorithm
,” in
*
Proceedings of the 14th Annual Conference on Computer Graphics and Interactive Techniques
*
(
ACM Press
,
1987
), pp.
163
–
169
.

**
28
. **
Z.
Wu
and
J. M.
Sullivan
Jr.
, “
Multiple material marching cubes algorithm
,”
Int. J. Numer. Methods Eng.
**
58
**
,
189
–
207
(
2003
).
[CrossRef]

**
29
. **
M.
Garland
and
P. S.
Heckbert
, “
Surface simplification using quadric error metrics
,” in
*
Proceedings of the 24th Annual Conference on Computer Graphics and Interactive Techniques
*
(
ACM Press
,
1997
), pp.
209
–
216
.

**
30
. **
J.
Bey
, “
Tetrahedral grid refinement
,”
Computing
**
55
**
,
355
–
378
(
1995
).
[CrossRef]

**
31
. **
L. H.
Wang
,
S. L.
Jacques
, and
L. Q.
Zheng
, “
MCML-Monte Carlo modeling of photon transport in multi-layered tissues
,”
Comput. Meth. Prog. Biomed.
**
47
**
,
131
–
146
(
1995
).
[CrossRef]

**
32
. **
D.
Boas
,
J.
Culver
,
J.
Stott
, and
A.
Dunn
, “
Three dimensional Monte Carlo code for photon migration through complex heterogeneous media including the adult human head
,”
Opt. Express
**
10
**
,
159
–
169
(
2002
),
http://www.opticsinfobase.org/abstract.cfm?URI=OPEX-10-3-159
.
[PubMed]

**
33
. **
User’s Manual (release 3.0) of TracePro (Software for Opto-Mechanical Modeling), Lambda Research Corporation, Littleton, MA.

**
34
. **
H.
Li
,
J.
Tian
,
F.
Zhu
,
W.
Cong
,
L. V.
Wang
,
E. A.
Hoffman
, and
G.
Wang
, “
A mouse optical simulation enviroment (MOSE) to investigate bioluminescent phenomena in the living mouse with the Monte Carlo Method
,”
Acad. Radiol.
**
11
**
,
1029
–
1038
(
2004
).
[CrossRef] [PubMed]

**
35
. **
J.
Huang
,
X.
Huang
,
D.
Metaxes
, and
D.
Banerjee
, “
3D tumor shape reconstruction from 2D bioluminescence images
,” presented in the third IEEE International Symposium on Biomedical Imaging (ISBI 2006),
Virginia, USA
, 6–9 Apr.
2006
.