1. Introduction

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

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

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

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

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

2. Methods

2.1 Forward problem

When fluorescent objects are embedded in turbid media, the fluorescence measurements U_m(r_s,r_d) detected at position r_d by a steady light source at position r_s with the emission wavelength λ_m can be formulated as [4,18]

where $D^{{\lambda }_m}$ is the diffusion coefficient with the emission wavelength λ_m. $k^{{\lambda }_x}$ is the wave number with the excitation wavelength λ_x, and $k^{{\lambda }_m}$ is the wave number with the emission wavelength λ_m. $U_x(r_s,r,k^{{\lambda }_x})$ is the theoretically calculated average photon intensity at position r, and $G(r_d-r,k^{{\lambda }_m})$ is the Green’s function that describes photon propagation from a point r to the detector at location r_d. S ₀ is a calibration factor that accounts for the system error such as laser power and the unknown gain and attenuation factors of the system. n(r) is the fluorescence distribution representing the fluorochrome concentration at location r multiplied by the fluorescent yield. For a single source-detector pair (s ₁,d ₁), Eq. (1) can be discretized in vector form as follows:

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

where W_ij is an element of the weight matrix W representing the weight of voxel j to measurement i.

2.2 Reconstruction method

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

Let {Θ₁,Θ₂,…,Θ_k} represent a sequence of hexahedron mesh levels of the volume of investigation, in which the mesh changes from coarse to fine with the increase in the mesh level k. Therefore, Eq. (3) on each level can be written as

where the index k represents the kth level. Practically, weight matrix W_k at the kth level is usually ill-posed, because of the inherent property of diffusion equation and the existence of noise. Therefore, it is usually impractical to solve for fluorescence distribution n_k from the linear system directly. In the algorithm proposed, an optimization approach is employed by minimizing objective function Ψ_k(n_k) based on the Tikhonov regularization method as follows:

where U_m(r_s,r_d) is the fluorescence measurements and λ_k is the regularization parameter. n^L_k and n^U_k are the lower and upper bounds of the fluorescence distribution at the kth level. Ψ_k(n_k) is optimized using Interior/CG [17]. Then, Eq. (5) can be written as

where s is a vector of slack variables and ν>0. The optimization consists of finding solutions of Eq. (6) for a sequence of positive barrier parameters {ν_l} that converges to zero. When switching to a finer mesh, the final solution n^r_k of the kth level is used to generate the initial value n ⁰ _k+1 of the (k+1)th level as follows:

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

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

2.3 Software implementation

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

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

3. Computational experiments

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

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

 

Fig. 1. Parallel plate FMT system setup. The schematic [5] of major components is shown in (a). Optical fibers’ arrangements are shown in (b) with black points, and the detection points’ arrangements are shown in (c).

Download Full Size | PDF

3.1 Single fluorescent target

In the experiment, fluorescent target of Cy5.5 fluorochrome was investigated. A laser source with light wavelength 675 nm was chosen, and the corresponding emission wavelength was 694 nm. The matching fluid had a reduced scattering coefficient µ′_s of 6 cm^-1 and an absorption coefficient µ_a of 0.3 cm^-1 for both the excitation and emission wavelengths. The optical properties used were consistent with estimates of the bulk optical properties of mouse tissues [20] for different organs or body regions, and they had been demonstrated successfully in phantom experiments [5]. A single small cubic fluorescent target was embedded at the position (x=0 cm, y=0 cm, z=0.75 cm), which was at a depth of 0.75 cm from the illumination surface. The size of the cubic fluorescent target was 0.01×0.01×0.01 cm. The fluorescence distribution n was set to unit 1 in the target and unit 0 in the background.

3.2 Dual fluorescent targets

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

4. Results

4.1 Single fluorescent target

A centrally located volume of size 1.8×1.2×1.0 cm was selected as the reconstruction region, as shown by the gray region of the imaging chamber in Fig. 1(a). In the reconstruction the lower bound was set to n^L_k=0 and the upper bound was set to n^U_k=100. Tomographic reconstructions were obtained as illustrated in Fig. 2(d) with zero initial guess n ⁰ ₁=n^L ₁ and regulation parameter λ_k=10^-9.

 

Fig. 2. Reconstruction of single fluorescent target. Mesh evolution in the reconstruction is demonstrated in (a), (b) and (c). Top 70% of the contour levels of reconstructed fluorescence distribution are shown in (d), in which the red cube represents the real target.

Download Full Size | PDF

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

In Fig. 2(d), the top 70% of the contour levels of reconstructed fluorescence distribution are shown, and the red cube represents the real target. The target center was reconstructed at the position (x=-0.02 cm, y=0.01 cm, z=0.75 cm), which is close to the real one. The reconstructed target fluorescence distribution n_recons was 1.1, with relative error (RE=|n_real-n_recons|/|n_real|) 10% to the real target fluorescence distribution n_real. The relative error is reasonable, because fluorescent targets with different fluorescence distribution and different size may have similar measurements, especially in the presence of noise. It could be illustrated by the 2.0% measurement error RE_meas=‖U _m,meas-U _m,cal‖/‖U _m,meas‖ between the measurements U _m,meas and the theoretically calculated data U _m.cal, which is equal to the Gaussian noise percentage added to the measurement data.

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

4.2 Dual fluorescent targets

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

 

Fig. 3. Reconstruction of two fluorescent targets for different edge-to-edge spacing: (a) 0.3 cm; (b) 0.15 cm; (c) 0.1 cm. The top 70% of the contour levels of reconstructed fluorescence distribution are shown, and red cubes represent the real targets.

Download Full Size | PDF

Table 1. Summary of results for dual fluorescent targets reconstructions. d_real is the edge-to-edge targets’ separation, and n_recons is the reconstructed targets’ fluorescence distribution.

View Table

5. Discussion and conclusion

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

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

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