1. Introduction

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

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

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

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

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

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

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

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

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

2. Methodology

2.1 The model of diffusion equations

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

where Φ_x,m denotes the photon density for excitation (subscript x) and fluorescence light (subscript m). In the linear model employed in this paper µ_ax,am is the absorption coefficient and D_x,m=1/3 (µ_ax,am+µ′_sx,sm) is the diffusion coefficient of the tissue [3, 24]. The absorption and scattering of the excitation light, caused by the fluorescence probes, were assumed to have little effect on the distribution of Φ_x [3], as the fluorescence probes often occupy a very small volume compared to the total imaging area. The fluorescent yield ηµ_af(r) is the unknown fluorescence parameter to be reconstructed, where η is the quantum yield and µ_af is the absorption coefficient of fluorescence probe to excitation light.

Whereas the FMT with 360° projections is carried out by rotating the subject while the excitation source and the detection device are stationary [6, 7], for the reconstruction we consider the phantom or animal geometry model to be stationary while the source and the corresponding detected boundary part to be moving. On the right hand side of Eq. (1.1) r_sl(l=1,2,…, L) represents the different excitation point source positions with respect to the subject with the amplitude Θ_s.

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

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

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

2.2 Finite element discretization

The finite element method has been widely used for numerically solving the diffusion equations [18, 23, 24], especially with arbitrary geometries [18]. Based on the finite element theory [26], Φ_x(r) and Φ_m(r) in Eq. (1) can be equivalent obtained by solving the following weak form equations:

by implementing the test function ψ(r).

Discretizing the domain Ω with N_p vertex node and N_e elements and employing ψ_i(r) as the shape function, Φ_x,m(r) can be represented by

where Φ_xi,mi denotes the nodal value of Φ_x,m(r)on the vertex V_i.

Let $\{{\gamma }_1,{\gamma }_2,\cdots ,{\gamma }_{N_p}\}$ be the basis functions and the unknown fluorescent yield ηµ_af(r) be denoted as x(r). Then it can be approximately expressed as

where γ_i could be chosen to be the same as ψ_i.

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

with

$\{\begin{array}{l}{K}_{\mathrm{xi},\mathrm{xj}}={\int }_{\Omega }(D_x\nabla {\psi }_{\mathrm{xi}}·\nabla {\psi }_{\mathrm{xj}}+{\mu }_{\mathrm{ax}}{\psi }_{\mathrm{xi}}{\psi }_{\mathrm{xj}})dr+\frac{1}{2q}{\int }_{\partial \Omega }{\psi }_{\mathrm{xi}}{\psi }_{\mathrm{xj}}dr\\ L_{\mathrm{xi},\mathrm{xj}}={\int }_{\Omega }{\Theta }_s\delta \left(r-r_s\right){\psi }_{\mathrm{xj}}dr\end{array}$

and

$\{\begin{array}{l}{K}_{\mathrm{mi},\mathrm{mj}}={\int }_{\Omega }(D_m\nabla {\psi }_{\mathrm{mi}}·\nabla {\psi }_{\mathrm{mj}}+{\mu }_{\mathrm{am}}{\psi }_{\mathrm{mi}}{\psi }_{\mathrm{mj}})dr+\frac{1}{2q}{\int }_{\partial \Omega }{\psi }_{\mathrm{mi}}{\psi }_{\mathrm{mj}}dr\\ f_{i,j}={\int }_{\Omega }{Ф}_x\left(\mathrm{r}\right){\psi }_{\mathrm{mi}}{\psi }_{\mathrm{mj}}dr\end{array}.$

2.3 Generation of the linear scheme

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

Step1. For each excited source position at r_sl(l=1,2,…, L) with respect to the subject, the corresponding Φ_x,sl is directly obtained by solving Eq. (6).

Step2. As the stiff matrix K_m is symmetrical and positive definite, Eq. (7) can be changed to

for a specific excitation light distribution Φ_x,sl.

Step3. {Φ_m,sl} can be divided into two parts {Φ^Meas _{m, sl}} and {Φ^NonM _m,sl} where the former consists of the vertices on the surface for measurement and the latter is the remainder of the vertices. Removing the corresponding rows as {Φ^NonM _{m, sl}} from [B_sl], the following matrix equation is formed as follows:

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

2.4 The hybrid adaptive finite element (AFE) reconstruction algorithm

2.4.1 Background

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

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

2.4.2 The hybrid AFE reconstruction algorithm

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

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

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

Firstly the imaging domain is discretized into a uniformly coarse mesh Θ₁, where there is no a priori knowledge about the distribution of fluorescence targets and all the vertices are included in the inversion. As an effective iterative regularization method, the CG algorithm is employed in the first reconstruction procedure. It is applied to the normal equation A ^T ₁ A ₁ X ₁=A ^T ₁Φ^mea, where the low-frequency components of the solution tend to converge faster than the high-frequency components. Hence, it has some inherent regularization effect where the number of iterations plays the role of the regularization parameter [28]. The initial value of X ₁ is set to zero and the optimized iteration number is determined by L-curve criterion [29] or by experience. The CG method converges very fast and is able to find an approximate value close to the real distribution in only a few iteration steps.

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

When switching to a finer mesh, the initial value X ⁰ _k+1 of the (k+1)th (k>=1) mesh level inherits the final solution X_k of the k th level by linear interpolation as follows:

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

 

Fig. 1. Flow diagram for the hybrid adaptive finite element reconstruction algorithm

Download Full Size | PDF

 

Fig. 2. Mesh refinement method. (a) is the original tetrahedral element. (b) is the second generation elements after refinement.

Download Full Size | PDF

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

where the superscripts of ^{X n+1} k ⁺¹ and Xnk ₊₁ denote the iteration numbers of Landweber method and α=1/λ _max with λ _max, the maximum eigenvalue of A^perk ₊₁·(A^per _k+1)^T. The terminating iteration number of Landweber method is determined by both the discrepancy principle and the known value range of fluorescent yield [31]. Usually Landweber method converges slowly. However, as a permissible region for the location of fluorescence targets is determined in the mesh levels other than the first one and the size of the corresponding weighted matrix is reduced, the convergence process is accelerated in this situation. Additionally, the reconstruction value is also stabilized, due to the small step length between neighboring iterations.

3. Simulations and results

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

 

Fig. 3. Experimental setup and the geometrical model for a mouse abdomen part. The sketch of the free-space FMT system is illustrated in (a), which is similar to that in Ref. [6]. The plane of excitation sources is shown in (b) where the black points represent positions of the isotropic point sources. The points 01–y015 are the first excitation sources for the 15 projections and the points 016–030 are the second ones. The field of view (FOV) of 150° for detection with respect to the excitation source point 01 is also shown. (c) and (d) are different views of the mouse geometry model used in the reconstruction. In order to reduce the boundary artifacts, which interfere with the finite element computation, the model was generated by sampling the original atlas data (intersections of the vertical and horizontal lines) and then approximating the curves using spline function to form the kidney surface and the body surface.

Download Full Size | PDF

3.1 Experimental setup

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

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

In this simulation fluorescence tomography of 360° full view was performed using 15 projections with 150° field of view detected for each projection [35]. The collimated laser beam was modeled as an isotropic point source, located one mean free path of photon transport (1/µ′_s) beneath the surface [25]. For each projection there were 2 point source positions of 18° interval at the z=7.5 mm plane, as shown in Fig. 3(b). Therefore, two images were collected for each projection.

 

Fig. 4. The forward geometry model for generating surface measurements is shown in (a) and (b) with a fluorescent probe in the left kidney, where the red part represents the fluorescent probe. Different views for the mesh of tetrahedral elements are shown in (c) and (d).

Download Full Size | PDF

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

3.2 Reconstruction for a single target

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

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

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

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

 

Fig. 5. (a). and (b). show different surface views of the reconstruction mesh in the initial mesh level.(c) plots of the Residual Error in the three mesh levels. See text for details.

Download Full Size | PDF

 

Fig. 6. Reconstruction results for a single fluorescence target embedded in the left kidney using the proposed algorithm. (a) the 3D reconstructed results obtained from the initial uniformly coarse mesh, with a threshold of 70% of the maximum value. (The red cylinders denote the real target). (b) the transverse view of the reconstruction at z=77 in the same mesh, where the white circle indicates the real fluorescence target. (c) and (d) the corresponding reconstruction using the refined, second mesh level. (e) and (f) the final reconstruction utilizing our algorithm, which is based on the third mesh level.

Download Full Size | PDF

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

 

Fig. 7. Reconstruction results for a single fluorescence target outside the kidneys using the proposed algorithm. (a) the 3D reconstructed results obtained from the initial uniformly coarse mesh, with a threshold of 70% of the maximum value. (The red cylinders denote the real target). (b) the transverse view of the reconstruction at z=77 in the same mesh, where the white circle indicates the real fluorescence target. (c) and (d) the corresponding reconstruction using the refined, second mesh level. (e) and (f) the final reconstruction utilizing our algorithm, which is based on the third mesh level.

Download Full Size | PDF

Table 1. Summary of the reconstruction results for a single target

View Table | View all tables in this article

3.3 Reconstruction of double targets

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

Table 2. Summary of the reconstruction results for double targets

View Table | View all tables in this article

3.4 Comparison with a fixed mesh

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

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

 

Fig. 8. Reconstruction results of double fluorescent targets in different discretized levels using the proposed algorithm. (a) the 3D reconstructed results obtained from the initial uniformly coarse mesh, with a threshold of 70% of the maximum value. (The red cylinders denote the real target). (b) the transverse view of the reconstruction at z=77 in the same mesh, where the white circle indicates the real fluorescence target. (c) and (d) the corresponding reconstruction using the refined, second mesh level. (e) and (f) the final reconstruction utilizing our algorithm, which is based on the third mesh level.

Download Full Size | PDF

 

Fig. 9. The transverse view of the reconstructed results at z=77, based on a fixed fine mesh with 7339 nodes and 34147 elements. For a single fluorescence target with the same size and location as in Fig. 6, the reconstruction results using CG is shown in (a), with a threshold of 70% of the maximum value. The corresponding results for two fluorescence targets with the same sizes and locations as in Fig. 8, with the same threshold as (a) is shown in (b),

Download Full Size | PDF

Table 3. Comparison of reconstruction in a fixed mesh and using the proposed algorithm

View Table | View all tables in this article

5. Discussion and conclusion

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

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

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

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

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

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