1. Introduction

The field of molecular imaging has gained considerable attention over the past few years due to its ability to reveal the molecular and cellular information in vivo. Hence, it becomes an increasingly important instrument for biomedical researchers to diagnose diseases, evaluate therapies, and facilitate drug development with small animals such as mouse models [1, 2, 3, 4]. As a mode of molecular imaging, optical molecular imaging, especially bioluminescence tomography (BLT), is outstanding because of its high performance, low cost and non-invasion [3, 5]. BLT is an emerging imaging technique which has been only recently developed to recover bioluminescent source distribution inside a living small animal in 3D [4, 6].

The bioluminescent signal is emitted when luciferin is combined with luciferase in the presence of oxygen and ATP, and luciferase enzymes are generally from firefly (FLuc), click beetle (CBGr68, CBRed), and Renilla reniformins (hRLuc) [7]. So signal has different emission spectra, roughly from 400nm to 750nm, which can be detected using a sensitive low light imaging systems, typically based on charge-coupled device (CCD) camera. Bioluminescent signal observed on the surface of the small animal forms the basis for tomographic reconstruction of internal bioluminescent source. However, three-dimensional (3D) bioluminescent reconstruction from boundary data is not unique and a highly ill-posed inverse problem in the general case [5].

The uniqueness research of BLT shows that a priori information has a quite effect on source reconstruction [5]. The commonly used priori information is the permissible source region strategy. So far, a priori permissible source region strategy has been developed for BLT reconstruction [8], which is inferred by the surface light power distribution and the heterogeneous structure of the phantom. Although promising results have been reported [4, 6, 8], it is not always reliable to infer such a permissible region especially when a single or multiple bioluminescent sources locate at half-radius or deeper depth from the small animal surface [9]. Recently, a posteriori permissible source region method has also been developed [9]. However, it takes several hours to select the permissible source region in the case of complex heterogenous phantom. Therefore, there is a critical need to establish a fast mathematical method for selecting the permissible source region.

Aside from the permissible source region strategy, the spectral characteristic of bioluminescent source can be considered as a priori information to improve BLT reconstruction [9, 10, 11, 12, 13, 14, 15], so there is an increasing interest in multispectral bioluminescence tomography. However, bioluminescent signal acquisition generally takes much longer time than that for fluorescent imaging, one exposure normally needs 5 to 10 minutes if the bioluminescent source is deep inside a small animal [16]. Although the bioluminescent signal is generally supposed to be stable in previous studies; in fact, it decays when the exposure time exceeds over one hour [16]. Traditionally, to collect the whole surface images around the imaged object, a rotating mechanism is employed to take four views (front, back, and two sides) [4, 6, 8, 16, 17]. So if we take all images in the multispectral case, one hour is not sufficient and it is possible that the observed boundary data is inaccurate. Therefore, in the multispectral case, the exposure time should be reduced with some methods. In addition, despite multispectral data acquisition could enhance the reconstruction result, dimensional disaster also arises from multispectral measurement and the large-scale data set seriously affects the reconstruction speed.

The development of a fast and robust BLT reconstruction algorithm based on practical application to biological research is necessary. In practice, in view to physical limitations as X-ray computed tomography (CT), the measurement for BLT is angle-limited and the acquired data are often incomplete [18, 19]. In this situation, partial data acquisition methods such as the single-view image technique could work. Recently, the feasibility of BLT with partial measured data has been proved theoretically and two reconstruction methods are extended to this case, but they are time-consuming and need about 5–6 hours in a single wavelength case [19]. The reconstruction time should be longer in the multispectral imaging because the computation data had increased rapidly.

Based on the consideration of above problems, we propose a reconstruction algorithm for BLT from multispectral partial boundary measurement. The most distinct difference between current and earlier tomography algorithms is the development of an optimal permissible source region strategy for reducing the ill-posedness of BLT, which is selected using iterative approach in combination of multispectral measurement without human intervention. In the algorithm, partial boundary measured data is used to reduce the exposure time which arise from multispectral measurement. Adaptive finite element algorithm is employed to reconstruct the underlying source distribution based on diffusion equation approximation. Then the linear relationship between the measured data and unknown source is established through the optimal permissible source region concept. In the optimization procedure, a spectral projected gradient-based optimization is used to reduce the reconstruction time and improve the robustness. In the numerical experiments, a heterogeneous phantom is designed to evaluate the performance of the algorithm. In the end, we discuss the relevant issues and conclude the paper.

2. Methods

2.1. Formulation of BLT

In bioluminescence imaging, it is important to depict the propagation of the photon transport accurately in biological tissue. The bioluminescent photon is subject to both scattering and absorption. However, in the 400nm–800nmwavelength range, photon scattering predominates over photon absorption. When the bioluminescence imaging experiment is performed in a dark environment, the propagation of photons can be described by the steady-state diffusion equation and Robin boundary condition [8, 20, 21]. Taking the influence of light wavelength λ on tissue optical property into account, the following model is given [9, 12]:

where Ω is a bounded smooth domain in the three-dimensional Euclidean space R ³ that contains an object to be imaged; ∂Ω is the corresponding boundary; Φ(x, λ) denotes the photon flux density [Watts/mm ²]; S(x, λ) is the bioluminescent source density [Watts/mm ³]; µ _a(x, λ) is the absorption coefficient [mm ^-1]; D(x, λ)=1/(3(µ _a(x, λ)+(1-g)µ _s(x, λ))) is the optical diffusion coefficient, µ _s(x, λ) the scattering coefficient [mm ^-1], and g the anisotropy parameter; ν(x) 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 [22]. Supposed that the measurement is conducted on some disjoint surface patches γ_j ⊂ ∂Ω, j=1,2,…,M, each is smooth and connected [19]. Let γ=⋃^M _j=1 γ_j, the measured quantity is the outgoing flux density Q(x, λ) on partial boundary γ and it can be expressed:

2.2. The optimal permissible source region strategy

In the practical experiment, the light which can reach the body surface of the small animals is separated into m bands τ ₁,…,τ _m using appropriate filters, with τ _l=[λ _l-1, λ _l), l=1,2,…, m-1, τ _m=[λ _m-1, λ _m]. Here λ ₀ <λ ₁ <…<λ _m is a partition of the spectrum range. Based on the finite element theory [23], the weak solution of the flux density Φ(x, τ _l) ∈ H ¹(Ω) on each band τ _l is given through the Eqs. (1) and (2):

where H ¹(Ω) is the Sobelev space and Ψ(x, τ _l) is an arbitrary piece-wise test function. In the framework of adaptive finite element analysis, let {T ₁,…,T _k,…} be a sequence of nested triangulation of the given domain Ω based on adaptive mesh refinement, where the sequence gradually changes from coarse to fine along with the increase in k [6]. The space of linear finite elements V _k are introduced on the discretized level T _k, satisfying V ₁ ⊂ … V _k ⊂ … … H ¹(Ω). Now, we only consider the kth discretized level which includes $N_{T_k}$ elements and $N_{P_k}$ vertex nodes. ${\psi }_1^k$,…, ${\psi }_{N_{P_K}}^k$ is the nodal basis of the space V _k. Then Eq. (5) can be simplified as following matrix form on the single-band τ _l [6, 9]:

where the components of the above matrices K _k(τ_l), C _k (τ _l), B _k(τ _l), F _k(τ _l) are given by

Let M _k(τ _l)=K _k(τ _l)+C _k(τ _l)+B _k(τ _l), in view to M _k(τ _l) is a sparse positive definite matrix, we have:

Let A _k(τ _l)=M ^-1 _k (τ _l)F _k(τ _l), by substituting A _k(τ _l) into the Eq. (8), we obtain the Eq. (9):

Generally, BLT problem is an ill-posed problem. One may consider to utilize the permissible source region as a priori knowledge to improve the BLT reconstruction [6, 8]. In order to reduce the ill-posedness of BLT, a novel optimal permissible source region strategy is developed by using iterative approach based algorithm before BLT reconstruction. In the permissible source region, the actual source may present and in the forbidden region, the actual source is zero by definition. Commonly, the energy contribution of each spectral band τ _l can be determined by performing a beforehand spectral analysis, that is S(τ _l)=ω(τ _l)S, where ω(τ _l) ≥ 0 and ∑^m _{l =1} ω(τ _l) ≈1, S denotes the total photon density. At the coarsest level (k=1), the whole reconstruction object Ω is assumed to be the permissible source solution S ₁, taking into account the above spectral distribution, based on Eq. (9), we have:

Where

$A=\left[\begin{array}{c}\omega \left({\tau }_1\right)A_1\left({\tau }_1\right)\\ \omega \left({\tau }_2\right)A_1\left({\tau }_2\right)\\ ⋮\\ \omega \left({\tau }_m\right)A_1\left({\tau }_m\right)\end{array}\right],\Phi =\left[\begin{array}{c}{\Phi }_1\left({\tau }_1\right)\\ {\Phi }_1\left({\tau }_2\right)\\ ⋮\\ {\Phi }_1\left({\tau }_m\right)\end{array}\right]$

A ^T is the transpose of A and A ^T A is an $m{N}_{T_1}\times m{N}_{T_1}$ symmetric matrix, then Eq. (10) has become a standard linear equation. Starting from an initial guess $S_1^0$ , the solution S ₁ can be updated iteratively by:

Where n is iteration number and β _n is the gradient of each iteration,

α _n is the step size, in order to get α _n, we construct a nonlinear least-squares optimization problem:

And

α _n can be gotten by directly solving the Eq. (14), that is:

In the iteration process, we directly use the measured data Φ^meas ₁ (τ _l) to substitute Φ₁(τ _l), because the effect of noise is low originally. With the increase of the number of iterations, the permissible source region will become smaller and smaller. When ‖β‖ _n≤ δ or iteration number n> N _max, the iteration is terminated and the rough region of the optimal solution S* is obtained, that is Ω*. We name it the optimal permissible source region and the sketch is shown in Fig. 1. Otherwise, S ₁ is updated by the Eq. (11). Here, N _max is maximum iteration number and the experiential typical value for δ is between 10^-2 and 10^-5. Generally, for N _max and δ, we choose 500 and 5×10^-3 for practical use, respectively.

 

Fig. 1. The sketch of searching optimal permissible source region.

Download Full Size | PDF

In BLT, taking into account the real physical meaning, the source density constrain can be taken as a priori information. Therefore, the nonnegative penalty is adopted during the iteration process.

2.3. Source reconstruction

After establishing the optimal permissible source region, reformulate the Eq. (8): Φk(τ _l)=M ^-1 _ei (τ _l)F _k(τ _l)S _k(τ _l), taking the linear relation between the unknown source value S _k and boundary measured points ${\Phi }_k^{meas}$ into consideration, we have:

Where G _k(τ _l) can be established by deleting the rows of [M _k(τ _l)^-1 F _k(τ _l)] corresponding to unmeasured points. Incorporating the optimal permissible source region, we have:

Where

${\Phi }_k^{\mathrm{meas}}=\left[\begin{array}{c}{\Phi }_k^{\mathrm{meas}}\left({\tau }_1\right)\\ {\Phi }_k^{\mathrm{meas}}\left({\tau }_2\right)\\ ⋮\\ {\Phi }_k^{\mathrm{meas}}\left({\tau }_m\right)\end{array}\right],G_k=\left[\begin{array}{c}\omega \left({\tau }_1\right)G_k\left({\tau }_1\right)\\ \omega \left({\tau }_2\right)G_k\left({\tau }_2\right)\\ ⋮\\ \omega \left({\tau }_m\right)G_k\left({\tau }_m\right)\end{array}\right]$

And W _k is a diagonal matrix for selecting the permissible region, that is:

$W_k=\mathrm{Diag}\left(w_{k\left(11\right)},w_{k\left(22\right)},\dots ,w_{k\left(\mathrm{ii}\right)},\dots ,w_{k\left(N_{P_l}{N}_{P_l}\right)}\right)$ $w_{k\left(\mathrm{ii}\right)}=\{\begin{array}{cc}1& \{s_{k\left(i\right)}\ge {\gamma }^k{s}_k^{\max }\}\\ 0& \{s_{k\left(i\right)}<{\gamma }^k{s}_k^{\max }\}\end{array}$

When k=1, s _k(i) is the optimal solution S ^* and s ^max _k is its maximum. k ≥ 2, s _k(i) and s ^max _k are the reconstructed results prolonged from k-1th level and the corresponding maximum, that is

where I^k _k-1 is the ratio operator from k-1 to k, γ ^k is the ratio factor, S _k-1 is the reconstruction result on (k-1)th level. By retaining the columns of G _k W _k corresponding to the permissible region S ^p _k, the final form of the linear system between the measurable boundary flux Φ_meas^k and S^p_k is obtained:

Since the matrix A _k is a severely ill-conditioned matrix because of the ill-posedness of BLT, the surface measured data is corrupted by noise, it is not practical to directly solve for S ^p _k from linear system (19). Then the following kth objective function is established:

where s ^k _inf and s ^k _suf are the kth level lower and upper bounds of the source density; Λ is the weight matrix, VΛ=V ^TΛV; ρ _k the regularization parameter. S ^init _k is initial value at the kth level. For BLT, Θk(S ^p _k) is a large-scale optimization problem with box-bound constrains, there-fore, a spectral projected gradient-based large-scale optimization algorithm is modified to solve the least square problem [9, 24, 25, 26]. In the algorithm, as for posteriori error estimation and local mesh refinement, similar disposal methods are adopted according to [6, 9]. The current gradient norm $\parallel g_{\underset{.}{\Theta }k}^i\parallel$ and the iteration number N ⁱ _k on each level are selected as switch condition indexes for triggering local mesh refinement. The discrepancy between the measured and computational boundary nodal flux data and the number k of mesh refinement are used as stop criterion. Finally, the flow chart of the algorithm is shown in Fig. 2.

 

Fig. 2. The flow chart of the proposed algorithm.

Download Full Size | PDF

3. Numerical simulation

As far as phantom is concerned, it plays an important role in evaluating reconstruction algorithm and various phantoms have been designed for bioluminescence tomography. In previous studies, the assumption of a homogeneous optical background had demonstrated the disadvantage and limitation for BLT reconstruction [13]. It’s not proper and accurate to evaluate the reconstruction algorithm.

Therefore, a heterogeneous phantom of 30mm height and 10mm radius was designed to evaluate the proposed algorithm [6]. It was made up of four ellipsoids and one cylinder to represent muscle (M), lungs (L), heart (H), bone (B) and liver (Li), as shown in Fig. 3(a). Based on the emission spectral distribution, the spectrum [500nm–750nm] can be divided into five discrete bins with steps of 50nm. In the experiments, the optical properties of each component are assumed as priori information and optical property parameters of each bin are compiled in Table 1 [9].

 

Fig. 3. Heterogeneous phantom. (a) A heterogeneous phantom with Muscle, Bone, Heart, Lungs, Liver and a single source in right lung. (b) The discretized mesh used in MOSE. (c) The initial mesh used in reconstruction algorithm. (d) The cross-section, including Muscle, Bone, Lungs and a source. The four arrows show the directions of the CCD camera.

Download Full Size | PDF

Table 1. Optical property parameters in different wavelength range (Unit: mm^-1) [9].

View Table | View all tables in this article

Because BLT is an inverse source problem, it is necessary to acquire the surface measured data objectively and reliably. In view to inverse crime, the synthetic data was produced by molecular optical simulation environment (MOSE) which was developed using Monte Carlo method [27]. In the single source simulation, the solid spherical source with 1mm radius and source intensity of 0.238nano-Watts/mm ³ was centered at (-3,5,15) inside the right lung as shown in Fig. 3(a). In MOSE, the heterogeneous phantomwas discritized by triangular elements with an average element diameter of about 1mm, shown in Fig. 3(b). In the simulation, the bioluminescent source of each band was sampled by 1.0×10 ⁶ photons and the default initial guess S ⁰ ₁ for searching the optimal permissible source region was set 0.01. When k=1, the ratio factor γ ^k=0; k > 1, γ ^k was initially set 10^-1 for selecting permissible source region and was changed by multiplying a factor of 10.0 with the increase of level k. In the optimization process, s ^k _inf, s ^k _suf and S ^init _k are equal to 0.0, 10000 and 1.0×10^-5, respectively. The stopping threshold ε_Φ and the maximum number of mesh refinement K _max were set to 1.0×10^-8 and 3, respectively. Noted that in the reconstruction procedure, the initial discretization of phantom was coarse volumetric mesh shown in Fig. 3(c). In terms of the present condition, a single-view image technique is feasible for acquiring partial measured data. So in our simulations, the data only from the front view shown in Fig. 3(d) is used in case of partial measurement.

3.1. Experiment Results

 

Fig. 4. Reconstruction results. (a) Reconstructed result using monochromatic synthetic data in [500nm,550nm] with partial measurement; (c) Spectrally resolved BLT reconstruction with proposed algorithm; (e) Multispectral BLT reconstruction with the optimal permissible source region strategy and complete measurement; (g) Reconstruction result with the method proposed in literature [9]; (b), (d), (f) and (h) are the magnified right lung images of (a), (c), (e) and (g), respectively. The sphere is the actual source.

Download Full Size | PDF

In the partial measurement experiment, the default initial guess S ⁰ ₁ was used. As for switch condition indexes, the terminated gradient norm ε _g and the maximum iteration number N ^max _k were set to 1.0×10^-7 and 5000, respectively. ρ _k was 1.16×10^-9 at each level. First, using partial measured data in [500nm,550nm], the optimal permissible source region was established and the reconstruction was performed, corresponding result was shown in Fig. 4(a). Despite the use of optimal permissible source region, the monochromatic measurement-based tomographic reconstruction result was still undesirable no matter how we changed the initial guess S ⁰ ₁, δ, ρ _k or maximum iteration number N _max. However, using multispectral partial measured data, BLT reconstruction showed preferable result with the proposed algorithm after one level refinement and the result was shown in Fig. 4(c). By contrast, when complete measured data from four views was used for selecting the optimal permissible source region, the corresponding result was demonstrated in Fig. 4(e).

In order to qualify the results, both the absolute source position error (ε) and the relative source density error (ξ) are defined: $\epsilon =\sqrt{{\left(x-x_0\right)}^2+{\left(y-y_0\right)}^2+{\left(z-z_0\right)}^2}$, ξ=|S _recon-S _real|/S _real, where (x, y, z) is the reconstructed center of source and (x ₀, y ₀, z ₀) is the actual center of source. S _recon and S _real are reconstructed source density and actual source density, respectively, the unit is nano-Watts/mm ³. Quantitative results about both the location and density of source were shown in Table 2.

In order to compare with other algorithm, BLT reconstruction was also conducted with the

Table 2. Quantitative comparison between the actual and reconstructed source centers and energy densities with different methods. The unit of ε is mm.

View Table | View all tables in this article

method proposed in literature [9], the result was demonstrated in Fig. 4(g) and quantitative result in Table 2. Quantitative comparison shows that the proposed algorithm has a better performance. Otherwise, the reconstruction speed is also a standard for evaluating algorithm. In order to get the result in Fig. 4(g), the total time cost used the algorithm mentioned in literature [9] was about 11 hours, but our method took much less time, about 3 hours. This was because the system matrix at first level was reduced from 3350×6878 to 1690×5397 with the employment of partial measurement and the optimal permissible source region. These work was all implemented on the same desktop computer with Pentium 4 3.4GHz and 1GB RAM.

From the above comparisons, it is obvious that the proposed algorithm can recover the bioluminescent source distribution relative accurately. This means that the measured data can be acquired only on part body of a small animal, and then the time for acquiring data can be dropped. Meanwhile, the system matrix dimension can also be reduced. Furthermore, our algorithm needs relative shorter reconstruction time and it is more appropriate for practical application. Nevertheless, spectrally resolved BLT is indispensable and provides more sufficient a priori information.

3.1.2. BLT reconstruction with noisy data

In order to evaluate the stability and robust of the proposed algorithm, noise effect was considered. In the experiment, white Gaussian noise with different level was added to the synthetic data and corresponding reconstructions were carried out. The corresponding results were shown in Fig. 5 and Table 3. The results reveal the stability of proposed algorithm under noise. When the noise level is lower (5% and 10%), the method has a better performance not only in reconstruction position but also source density. However, the noise level is higher (20% and 50%), although the reconstructed source density is perturbed, the position is not affected. After we have carried out intensive numerical simulation, we find two reasons to explain this. One reason is that the noise effect is low originally for selecting the optimal permissible source region. The other is that source reconstruction benefits from the adaptive from-coarse-to-fine mesh sequence.

Table 3. Reconstruction results with the proposed algorithm under different noise data.

View Table | View all tables in this article

 

Fig. 5. Reconstruction results with the proposed algorithm under different noise data. (a) Noise level was 5% noise, and for (c), (e), (g) was 10%, 20%, 50%, respectively; (b), (d), (f) and (h) are the magnified right lung images of (a), (c), (e) and (g), respectively. The sphere is the actual source.

Download Full Size | PDF

3.1.3. Consideration of different initial guesses

It is well known that a big disadvantage of most iterative approach based algorithm is sensitive to the initial guess (S ⁰ ₁). It is regarded as a criterion to test the stability of the proposed algorithm. Due to the limited space, representative results using different initial guesses of iterative approach are given. When the initial guess was 0.01, the result had been reported, shown in Fig. 4(c). When the initial guesses were set 0 and 0.1, the reconstructed results all indicated the algorithm was stable and robust to initial guess. The results were shown in Fig. 6 and Table 4. When the initial guess was equal to zero, the maximum reconstruction density was 0.142nano-Watts/mm ³. Although the density had a great difference with real density, the recovered center position of bioluminescence source was satisfactory, that was (-1.54,5.35,14.57).

Table 4. Quantitative reconstruction results using different initial guesses

View Table | View all tables in this article

 

Fig. 6. Reconstruction results with the proposed algorithm in view of different initial guesses. (a) initial guess S ⁰ ₁ is zero and for (c) is 0.1; (b) and (d) are magnified right lung images of (a) and (c), respectively. The sphere is the actual source.

Download Full Size | PDF

3.1.4. Optical property errors consideration

Currently, with further investigation of BLT, multimodality imaging fusion has been become a hot spot [4, 9, 13, 28]. Besides those, Diffusion Optical Tomography (DOT) is also appropriate for BLT which can be used to acquire optical properties of mouse tissues. In previous algorithms, optical property parameters of biomedical tissues were generally taken as a priori information to reduce the ill-posedness of BLT and deal with the nonuniqueness of BLT [5, 9, 12]. With DOT in combination with BLT, the problem can be solved. In view to the error range of DOT reconstruction quantification accuracy [9, 29], ±50% errors for all tissues were used to evaluate the proposed algorithm and reconstruction results were demonstrated in Fig. 7. The reconstructed position of source for -50% and +50% were (-3.34,3.87,14.08) and (-1.64,5.55,15.18), respectively. The maximum absolute source position error is 1.50mm, which accounts for the capability for tolerating optical property errors.

 

Fig. 7. Optical property errors consideration in the proposed algorithm. (a) The reconstruction result with +50% errors; (c) The reconstruction result with -50% errors; (b) and (d) are the magnified reconstructed results with +50% and -50% errors, respectively. The sphere is the actual source.

Download Full Size | PDF

Furthermore, the possible effect of two sources closer togetherwas studied. The two spherical sources of 1.0mm radius and 0.238nano-Watts/mm ³ power density were placed in the right lung with edge-to-edge distance of 2mm. Reconstruction results were shown in Talbe 5. In the case of +50% error for all tissues, the center position of two sources both had an offset of about 2mm, but the sources could be distinguished. However, when the error for all tissues is -50%, the reconstructed center position of a source had a large offset and the two sources could not be resolved. The effect of two sources with different edge-to-edge distances is necessary to further explore.

Table 5. Comparison between the actual and reconstructed source centers with optical property errors. The unit of ε is mm.

View Table | View all tables in this article

4. Discussions and conclusion

Given the difficulty that there is no unique solution to the BLT problem in general case and the need for fast BLT reconstruction algorithm, we have developed the novel optimal permissible source region strategy and present a BLT reconstruction algorithm to identify a 3D bioluminescent source distribution, finally evaluate its performance in numerical simulation with quantitative results. Simulation studies have indicated that the algorithm can exhibit very favorable performance. In the reconstruction process, the optimal permissible source region has played an important role to reduce the ill-posedness of BLT as well as enhance the numerical stability and computational efficiency. In addition, in the case of hyper- and multi-spectral measurement, the exposure time increases and dimension disaster arises. With the proposed algorithm, the problems can be dealt successfully.

A challenge problem for BLT is the reconstruction speed. In the algorithm, an optimal permissible source region is firstly established before reconstruction, so reconstruction time and memory cost are reduced. Secondly, adaptive finite element algorithm improves the reconstruction speed. Furthermore, the finite-element based reconstruction method makes it possible to handle a more complex geometrical model such as the real mouse phantom. The heterogeneous phantom experiments have illustrated that the method is fairly robust with respect to both strong Gaussian noise and initial guess. More importantly, the data in the algorithm were produced using Monte Carlo-based method, and hence they are totally free of the well-known inverse crime [8, 30]. Therefore, the presented algorithm can be properly evaluated. In the BLT prototype, multimodality fusion such as DOT, Micro-CT and BLT has become a research hot spot. Taking into account the high tolerance of optical property errors, the presented algorithm is suitable for future multimodality BLT system. However, because the algorithm depends on multiple factors, it is possible to get favorite results with try and error method and the optimal setting is still under investigation. As measure of performance, we considered the difference of locations of the center as well as the difference in intensity between the reconstructed and real source. Our algorithm demonstrates a satisfactory performance in this case, however, considering the difference in volume as well would be as important. In other words, there are remaining issues needed to be considered.

In summery, we have presented a novel strategy for selecting the permissible source region, and proposed a BLT reconstruction algorithm for practical application. It can provide a better performance in view of reconstruction quality, robustness and speed. Now, we are making an effort to construct a multimodality BLT system based on the proposed reconstruction algorithm, the new system will provide the potential to enhance our understanding of disease and drug activity during preclinical and clinical drug development, which could aid decisions to select candidates that seem most likely to be successful or to halt the development of drugs that seem likely to ultimately fail [31]. Corresponding results will be reported later.