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 600nm 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 ²]; S(x) is the source energy density [Watts/mm ³]; μ_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.4399n ^-2 + 0.7099n ^-1 + 0.6681 +0.0636n [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 ¹(Ω) 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 ¹(Ω) and S_λ ∈L ²(Ω); Φ^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}$₁,… $\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}$₁ ⊂… ⊂ $\mathcal{V}$_k ⊂… ⊂ H ¹(Ω). Now, we only consider the kth 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 kth discretized level:

where ${\varphi }_i^k$ (x) is the ith node value. Let {${\gamma }_1^k$, …,γ^k_NSk} be the interpolation basis functions, S(x) is approximated by S_k (x):

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

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:

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 ⁱ _Θk ∥ and the distance between the last two steps ∥ dⁱ_Sk ∥ as the indexes, where g ⁱ _Θk =∇Θ_k:($S_k^i$ ) and dⁱ_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 ⁱ _Θk ∥ is less than ${\epsilon }_g^{\left(k\right)}$ or ∥ dⁱ_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 kth 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$\mathcal{V}$k and 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:

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.

 

Fig. 1. Irregular refinements for the green closure in 3D.

Download Full Size | PDF

3. Numerical studies

A series of computational experiments was designed to evaluate our multilevel adaptive BLT algorithm. A heterogeneous cylindrical phantom of 30mm 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.5mm.

 

Fig. 2. Heterogeneous phantom. (a) A heterogeneous phantom with a single light source, composed of muscle (green), lungs (blue), heart (carmine), bone (white), liver (pink) and source (red); (b) The discretized mesh of the phantom used in MOSE; (c) The initial mesh used in the adaptive reconstruction; (d) The middle cross-section, including muscle, lungs and bone.

Download Full Size | PDF

Table 1. Optical parameters of the heterogeneous phantom.

View Table | View all tables in this article

3.2. Light source reconstruction

 

Fig. 3. Four views of the phantom surface with an angular increment of 90 degrees. Red lines denote the isoline of the surface light power.

Download Full Size | PDF

First, a single bioluminescent source was embedded in the right lung. In the MOSE simulation, a spherical source with 1.0mm 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 × 106 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)=XTX 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 Lmax 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

as the permissible source region.

 

Fig. 4. Comparison between the actual and reconstructed sources. (a) BLT reconstruction in the case of a single light source, and (b) BLT reconstruction in the case of four light sources.

Download Full Size | PDF

 

Fig. 5. Mesh evolution for the permissible source region in the single light source case. The green mesh denotes the permissible source region; The red sphere is the actual source. (a) The first-level reconstruction, (b) and (c) The second- and third-level results, respectively.

Download Full Size | PDF

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 ³, 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 ³. 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 4mm, 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.

Table 2. Comparison between the fixed mesh of three different discretized scales and the adaptive mesh. Mesh size denotes the average element diameter of reconstructed results.

View Table | View all tables in this article

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.0mm 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.

 

Fig. 6. BLT reconstruction in the fixed mesh of different descretized scales. The green mesh denotes the permissible source region; The red sphere is the actual source. (a) and (b) The average element diameter of the mesh is about 2mm and 1mm respectively.

Download Full Size | PDF

Table 3. Comparison between the actual and reconstructed sources in the four-source case.

View Table | View all tables in this article

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.0mm radius and 0.238nano-Watts/mm ³ 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.

 

Fig. 7. Numerical study on the spatial resolution of the multilevel adaptive algorithm. (a) The BLT reconstruction in the case of 1.5mm separation; (b) The counterpart in the case of 1.0mm separation; and (c) That in the case of 0.5mm separation.

Download Full Size | PDF

Table 4. Comparison between the actual and reconstructed source centers and energy densities.

View Table | View all tables in this article

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.