PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of bmengonBioMed CentralBiomed Central Web Sitesearchsubmit a manuscriptregisterthis articleBioMedical Engineering OnLine
 
Biomed Eng Online. 2010; 9: 20.
Published online 2010 May 20. doi:  10.1186/1475-925X-9-20
PMCID: PMC2893184

Image reconstruction of fluorescent molecular tomography based on the tree structured Schur complement decomposition

Abstract

Background

The inverse problem of fluorescent molecular tomography (FMT) often involves complex large-scale matrix operations, which may lead to unacceptable computational errors and complexity. In this research, a tree structured Schur complement decomposition strategy is proposed to accelerate the reconstruction process and reduce the computational complexity. Additionally, an adaptive regularization scheme is developed to improve the ill-posedness of the inverse problem.

Methods

The global system is decomposed level by level with the Schur complement system along two paths in the tree structure. The resultant subsystems are solved in combination with the biconjugate gradient method. The mesh for the inverse problem is generated incorporating the prior information. During the reconstruction, the regularization parameters are adaptive not only to the spatial variations but also to the variations of the objective function to tackle the ill-posed nature of the inverse problem.

Results

Simulation results demonstrate that the strategy of the tree structured Schur complement decomposition obviously outperforms the previous methods, such as the conventional Conjugate-Gradient (CG) and the Schur CG methods, in both reconstruction accuracy and speed. As compared with the Tikhonov regularization method, the adaptive regularization scheme can significantly improve ill-posedness of the inverse problem.

Conclusions

The methods proposed in this paper can significantly improve the reconstructed image quality of FMT and accelerate the reconstruction process.

Background

Near-infrared (NIR) light can travel several centimeters through biological tissue, and hence has the potential to qualify the molecular information by fluorochromes in tissue [1]. Recently, there has been increasing interest in the molecularly-based medical imaging method, such as fluorescent molecular tomography (FMT) [2-4], in which the injected fluorophore may accumulate in diseased tissue. During the imaging process, the tissue surface is illuminated with excitation light. Then, the fluorophores are excited to emit the light, which is detected as fluorescence [5]. The process of fluorescent light generation and transportation through tissues can be described by a forward model, so that the surface measurements can be predicted on the basis of a guess of the system parameters and the given source positions. To reconstruct an image, it is necessary to calculate the internal optical and fluorescent properties with the given measured data and sources [6].

One of the major challenges in the reconstruction of FMT is its high computational complexity resulted from extremely large-scale matrix manipulations. Generally, the iterative solution approaches, such as CG method [7] and Gauss-Newton (GN) method [8], are more efficient than the direct solution approaches. Additionally, the iterative methods based on the reduced system can be more efficient than those based on the global system. One of such systems is the Schur complement system, which was firstly used by Haynsworth [9]. The condition number of the Schur complement of a matrix is never greater than that of the given matrix, and hence the convergence properties of iterative solving of linear systems can be significantly improved [7,10]. In this paper, we propose to adapt this idea for the FMT reconstruction. The most important innovation of our method lies in its tree structured level-by-level decomposition strategy, where decompositions in each level are performed in two ways. This strategy is quite different from that in [10] where only one component of the global solution is derived in the Schur complement system. The advantages of our method are obvious because a further improvement in the reconstruction accuracy and speed can be achieved with level-by-level Schur complement decomposition. Another contribution of this paper is that we propose a modified spatially variant regularization method incorporating the objective function to tackle the ill-posed nature of the inverse problem.

Methods

Forward Model and Finite Element Formulation

FMT acquisitions are obtained through a two-step image formation model [11]. In the first step, sources at several locations are used to illuminate the tissue. This step, in frequency domain, is driven by the diffusion equation [12]

equation image
(1)

where the subscript x denotes the excitation wavelength; [nabla] is the gradient operator; Sx(W/cm3) is the excitation light source; Φx(W/cm2), Dx(cm), and kx (cm-1) represent the photon fluence, the diffusion coefficient, and the decay coefficient, respectively; Ω denotes the bounded domain of reconstruction.

In the second step, the fluorophores are excited to emit the fluorescence. The second step can be modelled by a second diffusion equation

equation image
(2)

where the subscript m indicates the emission wavelength, ω(rad/s) denotes the modulation frequency of the source. Sm is the emission light source. The diffusion coefficient Dx,m(cm), and the decay coefficient kx,m (cm-1) are defined, respectively, as[6]

equation image
(3)
equation image
(4)

where μax,mi (cm-1) denotes the absorption coefficient due to endogenous chromophores; μax,mf (cm-1) represents the absorption coefficient due to exogenous fluorophores; An external file that holds a picture, illustration, etc.
Object name is 1475-925X-9-20-i5.gif is the reduced scattering coefficient; q is the quantum efficiency of the fluorophore; τ(s) is the lifetime of fluorescence; and finally, c(cm/s) is the speed of light in the medium.

Here, the Robin-type boundary conditions are implemented on the boundary [partial differential]Ω of domain Ω to solve the above diffusion equations

equation image
(5)
equation image
(6)

where n is a vector normal to the boundary [partial differential]Ω, bx,m is the Robin boundary coefficient.

To solve the forward problem within the finite element method (FEM) framework, the domain Ω is divided into P elements and joined at N vertex nodes. The solution Φx,m is approximated by the piecewise linear function An external file that holds a picture, illustration, etc.
Object name is 1475-925X-9-20-i8.gif, with [var phi]i (i = 1...N) being basis functions [13]. Hence, equations (1) and (2) can be rewritten as

equation image
(7)
equation image
(8)

where

equation image
(9)
equation image
(10)

The elements of finite element matrix Ax,m can be obtained from the formula

equation image
(11)

with Ωh and Γh being the bounded domain and its boundary, respectively.

Inverse Process of FMT

The inverse process of FMT is to estimate the spatial distribution of the optical or fluorescent properties of the tissues from measurements [14]. In the discrete case, the reconstruction problem can be defined as the optimization of the objective function

equation image
(12)

where G is the forward operator, || || is L2-norm, x and y are the calculated optical or fluorescent properties of the tissues and the detector readings, respectively.

Suppose that the objective function E attains its extremum at x + Δx, expanding the gradient of the objective function E' about x in a Taylor series and keeping up to the first-order term leads to

equation image
(13)

Equation (13) can be further written as [15]

equation image
(14)

where T denotes the transpose, Δy = y - G(x) is the residual data between the measurements and the predicted data. The Jacobian matrix J is a measure of the rate of change in measurement with respect to the optical parameters. It describes the influence of a voxel on a detector reading. H is the Hessian matrix, whose entries are the second-order partial derivatives of the function with respect to all unknown parameters describing the local curvature of the function with respect to many variables [16].

Introducing the Tikhonov regularization term to tackle the ill-posedness of the inverse problem and ignoring the Hessian matrix, the solution to the linearized reconstruction problem can be described as follows

equation image
(15)

where λ is a regularization parameter, which can be determined by the Morozov discrepancy principle [17], I is an identity matrix.

Adaptive Regularization Scheme

The problem of image reconstruction for FMT is ill-posed [18]. The Tikhonov regularization technique, as mentioned above, is one of the major methods to reduce the ill-posedness of the problem [19]. However, there exists one protrudent difficulty for this technique in the determination of the regularization parameter. A general unexpected characteristic of the NIR imaging is that the resolution and contrast of the reconstructed images degrade with the increased distance from the sources and the detectors [20]. Considering the fact that the value of the regularization parameter has important effect on the contrast and resolution of the resultant images, one strategy to solve this problem is to use a spatially variant regularization parameter. Meanwhile, it can be inferred that the objective function is related to the regularization parameters [15]. During the process of minimizing the objective function, decreasing λ will speed up the convergence if the value of objective function is decreasing, otherwise increasing λ can enlarge the searching area (trust-region). Upon the basis of these considerations, we propose a modified regularization method both adaptive to the spatial variations and the objective function.

Suppose that the number of measurements and the number of the vertex nodes are M and N, respectively. Thus, we have for the matrices in equation (15): Δx[set membership]RN × 1, Δy[set membership]RM × 1, J[set membership]RM × N, and I[set membership]RN × N. To construct a spatially variant regularization framework, the inverse term of (JTJ + λI)-1 in equation (15) is replaced with (JTJ + λ)-1, which results in the following equation

equation image
(16)

where An external file that holds a picture, illustration, etc.
Object name is 1475-925X-9-20-i19.gif is a diagonal matrix. Equation (16) can be rewritten as

equation image
(17)

with Δxi (i = 1, 2, ..., N) being the component of the vector Δx. It can be easily seen that each node pi (i = 1, 2,...,N) in the reconstructed domain is regularized by a corresponding regularization parameter λi (i = 1, 2,...,N) respectively. Obviously, the above mentioned Tikhonov regularization can be regarded as a special case of equation (17) when λ1 = λ2 = (...)λN = λ.

It was pointed out in [21] that, the resolution and contrast of the images decrease with the increment of the regularization parameters and vice versa. Therefore, the adaptive regularization parameter λi can be defined as follows to compensate the decrease of the resolution and contrast with the increased distance from the sources and detectors:

equation image
(18)

where ri is the position of node pi, rs and rm respectively denote the positions of the source and detector closest to the node pi, c1 and c2 are two positive parameters determined empirically in our paper.

To make the regularization parameter adaptive to the objective function as defined in equation (12), we propose to incorporate it in the regularization as follows

equation image
(19)

In equation (19), the arctan function is used to guarantee a relatively small fluctuation range of the regularization parameters and avoid too large values of them. Obviously, regularization parameters determined from equation (19) relate to the objective function in a similar manner to that as pointed out before. In such a way, the regularization parameters are adaptive not only to the spatial variations but also to the variations of the objective function to accelerate the convergence.

Reconstruction Based on the Schur Complement System

As has been pointed out previously, the iterative methods based on the Schur complement system can be more efficient to solve large-scale problems. Hence, we propose to reconstruct the tomographic image of FMT with level-by-level decomposition in the Schur complement system.

For convenience of discussions, equation (16) can be further rewritten in a more compact form as

equation image
(20)

where k = JTJ + λ and b = JT Δy.

To solve the inverse problem of FMT in the Schur complement system, the solution space Rn is firstly decomposed into two subspaces of U and V with dimensions m and n-m, respectively. Let [Γ Ψ] be an orthonormal basis of the solution space Rn. The basis of the m-dimensional coarse subspace U is formed by the columns of Γ [set membership] Rn × m and the columns of Ψ [set membership] Rn×(n-m) form the basis of the (n - m) dimensional subspace V.

Therefore, the solution to equation (20) can be expressed with the bases of the two subspaces as follows

equation image
(21)

where u and v are the projections of Δx on the subspaces U and V, respectively.

Because both the condition number and the scale of the system can be reduced after Schur complement decomposition, we propose to further decompose both the projections u and v level by level with the Schur complement decomposition along two paths in a tree structure, and then solve the subsystems in the Schur complement systems. Our approach is different from that proposed in [10], where only the projection v is solved in the Schur complement system. The level-by-level Schur complement decomposition can be schematically illustrated as in Figure Figure1.1. We derive the iterative system in the following discussions.

Figure 1
Schematic illustration of Schur complement decomposition with a tree structure. The global solution Δx is decomposed with the Schur complement system level by level along the two paths in the tree structure.

Suppose that the subsystem at the ith level is as follows

equation image
(22)

where S(i, j) is the Schur complement matrix with the subscript (i, j) being the jth (j = 0, 1,..., 2i) term at the ith (i = 0, 1,..., L) level in the tree structure as illustrated in Figure Figure1.1. Particularly, S(0,0) is the global matrix k as defined in equation (20). To solve this system in the Schur complement system, equation (22) will be further decomposed at the i+1th level. Thus, the solution Δx(i,j) is firstly expressed with the bases of the two subspaces as

equation image
(23)

where Δx(i+1,2j-1) and Δx(i+1,2j) are the projections of Δx(i,j) on the subspaces formed by the columns of Γ(i,j) and Ψ(i,j), respectively.

Substituting equation (23) into equation (22) yields

equation image
(24)

Multiplying both sides of equation (24) from the left by [Γ(i,j)Ψ(i,j)]T, we can obtain

equation image
(25)

Thus, equation (25) can be further rewritten into a two-by-two block system

equation image
(26)

where An external file that holds a picture, illustration, etc.
Object name is 1475-925X-9-20-i30.gif, An external file that holds a picture, illustration, etc.
Object name is 1475-925X-9-20-i31.gif, An external file that holds a picture, illustration, etc.
Object name is 1475-925X-9-20-i32.gif, and An external file that holds a picture, illustration, etc.
Object name is 1475-925X-9-20-i33.gif, while the two components on the right-hand side (RHS) of equation (26) are An external file that holds a picture, illustration, etc.
Object name is 1475-925X-9-20-i34.gif, and An external file that holds a picture, illustration, etc.
Object name is 1475-925X-9-20-i35.gif. From equation (26), it can be seen that S(i,j)11 and S(i,j)22 correspond to the equations for the unknowns of Δx(i+1,2j-1) and Δx(i+1,2j), respectively, while S(i,j)12 and S(i,j)21 define the coupling between these two sets, which will be eliminated in the following discussions.

Applying block Gaussian elimination to equation (26) leads to [22]

equation image
(27)

where An external file that holds a picture, illustration, etc.
Object name is 1475-925X-9-20-i37.gif, which is called the Schur complement with respect to S(i,j)11[7], An external file that holds a picture, illustration, etc.
Object name is 1475-925X-9-20-i38.gif. From equation (27), we have

equation image
(28)
equation image
(29)

It can be found that the condition number of matrix S(i+1,2j) is smaller than that of matrix S(i,j)[9]. Hence, solving the inverse problem in the Schur complement system at the i+1th level will be more efficient than solving it at the ith level. We herein solve equation (29) using the biconjugate gradient method [23]. Its advantage is that it does not square the condition number of the original equations [24]. Basically, the biconjugate gradient method can be used to solve the large-scale systems with the fastest speed among all the generalized conjugate gradient methods in many cases [25]. The algorithm for solving equation (29) can be summarized as follows

Algorithm 1

1. Input an initial guess Δx(i+1,2j)0;

2. Initialize d0 = f0 = r0 = p0b(i+1,2j) - S(i+1,2j)Δx(i+1,2j)0;

3. For k = 0, 1, 2... until convergence do

equation image

   End for

After the derivation of Δx(i+1,2j) from equation (29) with algorithm 1, the next task is to obtain the other component of Δx(i+1,2j-1) for the synthesis of the solution Δx(i,j). Here, Δx(i+1,2j-1) is also solved in the Schur complement system due to its low condition number.

Eliminating the block S(i,j)12 in equation (26) using block Gaussian elimination with S(i,j)22 as pivot block, we have

equation image
(30)

where An external file that holds a picture, illustration, etc.
Object name is 1475-925X-9-20-i43.gif and An external file that holds a picture, illustration, etc.
Object name is 1475-925X-9-20-i44.gif. Thus S(i+1,2j-1) is the Schur complement with respect to S(i,j)22.

From equation (30), we can obtain

equation image
(31)

Thus, the solution Δx(i+1,2j-1) can be obtained in a same manner as in Algorithm 1, and the only difference is that Δx(i+1,2j), S(i+1,2j), and b(i+1,2j) should be replaced with Δx(i+1,2j-1), S(i+1,2j-1), and b(i+1,2j-1), respectively. Solving equation (31) is computationally efficient because of the reduced condition number in the Schur complement system [7]. Moreover, such a strategy of deriving both Δx(i+1,2j-1) and Δx(i+1,2j) in the Schur complement system can be implemented in a parallel manner, since equations (29) and (31) are decoupled. Therefore the subsystem at the ith level as in equation (22) can be decomposed into the two linear subsystems at the i+1th level, i.e., Schur complement systems as in equations (29) and (31). After obtaining Δx(i+1,2j-1) and Δx(i+1,2j), they are then substituted into equation (23) to yield the solution Δx(i,j) at the ith level. The whole reconstruction algorithm is summarized as follows

Algorithm 2

1. Set x0 to an initial guess;

2. x x0, calculate b and k at x in equation (20) with the adaptive regularization scheme as in equation (19);

3. The global system of equation (20) is decomposed with the Schur complement system level by level in a same manner as the decomposition of equation (22) into equations (29) and (31) to obtain the subsystem S(i,j)Δx(i,j) = b(i,j) at the ith level for i =1,..., L and j =1,..., 2i, the subspaces at the ith level are formed by the columns of Γ(i,j) and Ψ(i,j), respectively;

4. Set i = L;

   For j = 1,..., 2i do

      Combining equations (26), (27), and (30), solve S(i+1,2j)Δx(i+1,2j) = b(i+1,2j) and S(i+1,2j-1)Δx(i+1,2j-1) = b(i+1,2j-1) with Algorithm 1, where Δx(i+1,2j-1), S(i+1,2j-1), and b(i+1,2j-1) are used instead of Δx(i+1,2j), S(i+1,2j), and b(i+1,2j) when Algorithm 1 is employed for the latter case;

   End for

5. While i ≥ 0

   {

   For j = 1,..., 2i do

      Substitute the solutions Δx(i+1,2j) and Δx(i+1,2j-1) into equation (23) to obtain the solution Δx(i,j) at the ith level;

   End for

   i = i - 1;

   }

6. x0x0 + Δx(0,1);

7. If ||Δx(0,1)|| >ε

         go to 2;

   Else

         x x0, output x.

As mentioned before, the Schur complement system has a smaller condition number than that of the system from which it is constructed [7]. As a result, iterative methods based on the Schur complement systems can be more efficient than the methods based on the global matrix as in equation (20) due to its reduced scale and the smaller condition number. Therefore, the proposed algorithm can be expected to be more efficient than the conventional ones, as the results demonstrated in the next section.

Results and Discussion

In this work, assuming that the scattering coefficients are known, we focus on the reconstruction of the absorption coefficient μaxf. Two phantoms as illustrated in Figure Figure22 are used to evaluate the proposed algorithm. Figure 2(a) contains one object, and Figure 2(b) contains two objects of different shapes. Table Table11 and Table Table22 outline the optical and fluorescent parameters in different regions of the simulated phantoms corresponding to Figures 2(a) and 2(b), respectively. Four sources and thirty detectors are equally distributed around the circumference of the simulated phantom. The simulated forward data are obtained from equations (1) and (2), in which Gaussian noise with a signal-to-noise ratio of 10dB is added to evaluate the noise robustness of the algorithms. The parameters c1 and c2 in equation (19) are, respectively, set to 0.2 and 2. The initial guesses for solutions Δx(i+1,2j) and Δx(i+1,2j-1) of equations (29) and (31) are set to 0. The initial value of x0 is set to 5 mm-1. The subspace created from the right singular vectors of the singular value decomposition (SVD) is optimal. Since SVD is computationally expensive, it is expected that a subspace close to SVD subspace will do almost as good. Thus, the choice of an oscillatory basis can be a basis created by sine or cosine functions with increasing frequency [26]. Here discrete cosine basis is employed in the simulations. To reliably evaluate the performance of different methods for the inverse problem, the best way is to use an independent forward model, which is different from the one employed in the inverse problem, to generate the synthetic data [27]. Therefore, in our case, a finer mesh as shown in Figure Figure33 with 169 nodes and 294 triangular elements is used to generate the forward simulated data.

Figure 2
Simulated phantoms for FMT. (a) One object with absorption coefficient μaxf of 0.2 mm-1 on the background medium with μaxf of 0.06 mm-1, and (b) two objects of different shapes, one of which with high absorption coefficient of 0.2 mm-1 ...
Table 1
Optical and fluorescent properties of one-object phantom
Table 2
Optical and fluorescent properties of two-object phantom
Figure 3
Mesh used for forward solver of FMT. The mesh contains 169 nodes and 294 triangular elements.

It is well known that the most significant superiority of the anatomical imaging modality lies in the high spatial resolution. Hence, it will be helpful to improve the image quality and accelerate the reconstruction process if we use the anatomical image as prior information for mesh generation. The reconstructed domain is firstly uniformly discretized according to the Delaunay triangulation scheme, after which the uniform mesh is refined only for the areas with large variations of the pixel values. To simulate this idea, we employ the images shown in Figures 4(a) and 4(b) with a resolution of 100 × 100 pixels as the prior images corresponding to Figures 2(a) and 2(b), respectively. The meshes are generated as shown in Figure Figure55 for the inverse problem of FMT. The mesh with 122 nodes and 212 triangular elements (Figure 5(a)), and the mesh with 148 nodes and 264 triangular elements (Figure 5(b)) are generated incorporating the prior information as shown in Figures 4(a) and 4(b), respectively.

Figure 4
Model of prior image. (a) The prior image corresponds to Figure 2(a), and (b) the prior image corresponds to Figure 2(b).
Figure 5
Meshes used for inverse problem of FMT. (a) Mesh with 122 nodes and 212 triangular elements, and (b) mesh with 148 nodes and 264 triangular elements..

Figures 6(a) and 6(b) show, respectively, the reconstructed images of μaxf for one object phantom using the adaptive regularization scheme and Tikhonov regularization method. Figures 7(a) and 7(b) depict the results for two objects phantom from the above two different algorithms. Here, both of the results from Figures Figures66 and and77 are based on the CG method. As seen from Figures Figures66 and and7,7, better reconstructed results can be achieved from the adaptive regularization scheme. To quantitatively assess the accuracy of the different algorithms, the mean square error (MSE) is introduced

equation image
(32)
Figure 6
Reconstructed images of absorption coefficient due to exogenous fluorophores μaxf for one object. (a) Reconstructed result with the adaptive regularization scheme, and (b) reconstructed result with the Tikhonov regularization method.
Figure 7
Reconstructed images of absorption coefficient due to exogenous fluorophores μaxf for two objects. (a) Reconstructed result with the adaptive regularization scheme, and (b) reconstructed result with the Tikhonov regularization method.

where N is the total number of nodes in the domain. The superscript calc denotes the values obtained using reconstruction algorithms; and actual denotes the actual distribution of μaxf which is used to generate the synthetic image data set. Table Table33 lists the performance of the reconstruction algorithms in terms of MSE. It can be seen that the adaptive regularization scheme can significantly improve the quality of the reconstructed images and achieve a smaller MSE in either case.

Table 3
Comparison of performance of methods

Figure Figure88 shows the reconstructed images of μaxf for one object phantom using the different algorithms after 1, 15, and 30 iterations, respectively. After 30 iterations, the reconstructed image from the proposed algorithm has a relatively higher contrast than those obtained from the other two algorithms. Figure Figure99 depicts the reconstructed images of μaxf for two objects phantom using the different algorithms, from which it can be seen that the proposed method can reconstruct the images more accurately than the other two methods even after the first iteration. According to the third column of Figure Figure9,9, the reconstructed image quality based on our algorithm is significantly improved as compared with that based on the other two methods.

Figure 8
Reconstructed images of absorption coefficient due to exogenous fluorophores μaxf for one object. (a)-(c) with the CG method, (d)-(f) with the Schur CG method, and (g)-(i) with our method. The reconstruction results after one iteration are shown ...
Figure 9
Reconstructed images of absorption coefficient due to exogenous fluorophores μaxf for two objects.(a)-(c) with the CG method, (d)-(f) with the Schur CG method, and (g)-(i) with our method. The reconstruction results after one iteration are shown ...

We investigated how the MSE changed against the number of iterations for different algorithms. Figure Figure1010 shows a fast convergence of our algorithm with a less MSE than the other two algorithms. In addition, the CG method converges slower than the Schur CG method and our method, which means that solving the inverse problem based on the Schur complement system is superior to that based on the global system. The computation time of different algorithms is further investigated in our work to evaluate the convergence rate. Table Table44 lists the computation time after 30 iterations for different algorithms. From this table, it can be seen that the time needed for our algorithm after 30 iterations is less than that of the Schur CG method. Although the former is a little bit longer than the time needed for the CG method, our algorithm needs only less than 5 iterations to achieve the precision of the CG method after 30 iterations. As a result, the CG method needs much more iterations to achieve a given precision of reconstruction than our method. Therefore, compared with the other two methods, the proposed algorithm is more efficient and stable.

Figure 10
Number of iterations versus MSE between the original and the reconstructed images. (a) Reconstruction of one object, and (b) reconstruction of two objects.
Table 4
Computation time of FMT image reconstruction for 30 Iterations

To further validate the proposed algorithm for 3D reconstruction, a phantom as illustrated in Figure Figure1111 is used for simulations. Within this phantom, a small cylindrical object is suspended. In Figure Figure11,11, the dashed curves represent the planes of measurements. Four sources and sixteen measurements are used for each plane in the simulations. The mesh for reconstructing the 3D image is shown in Figure Figure12,12, which contains 858 nodes and 3208 tetrahedral elements. Figures Figures1313 and and1414 depict the reconstructed 2D cross sections of the 3D phantom shown in Figure Figure1111 using the Schur CG method and the proposed algorithm, respectively. Table Table55 lists the performance of the above two methods for a quantitative comparison. From this table, we can conclude that our proposed algorithm can also speed up the reconstruction process and achieve high accuracy for the 3D case.

Figure 11
Simulated phantom for 3D reconstruction. The phantom of radius 10 mm and height 40 mm with a uniform background of μaxf = 0.005mm-1 , which is positioned at x = 10mm, y = 0mm and z = 20mm. The small cylindrical anomaly has a radius of 2 mm and ...
Figure 12
3D mesh for image reconstruction. 3D mesh for image reconstruction with 858 nodes and 3208 tetrahedral elements.
Figure 13
Reconstructed images using the Schur CG method. Reconstructed images using the Schur CG method, which are 2D cross sections through the reconstructed 3D volume. The right-hand side corresponds to the top of the cylinder (z = 40 mm), and the left corresponds ...
Figure 14
Reconstructed images using the proposed algorithm. Reconstructed images using the proposed algorithm, which are 2D cross sections through the reconstructed 3D volume. The right-hand side corresponds to the top of the cylinder (z = 40 mm), and the left ...
Table 5
Performance comparison of reconstruction methods for 3D case

Conclusion

In this paper, we developed a novel image reconstruction method of FMT, based on the tree structured Schur complement decomposition in combination with the adaptive regularization scheme. The proposed approach decomposes the global inverse problem level by level with the Schur complement decomposition, and the resultant subsystems are solved with the biconjugate gradient method. The spatially variant regularization parameter is determined adaptively according to the objective function. Simulation results demonstrate that the proposed method outperforms the previous methods, such as the CG and the Schur CG methods, in both reconstruction accuracy and speed.

Competing interests

The authors declare that they have no competing interests.

Authors' contributions

WZ conceived the study, implemented the algorithm, and drafted the initial manuscript. JJW participated in the design of the study, analyzed the simulation results, and helped to draft the manuscript. DDF critically reviewed the manuscript, helped to analyze the results, revised the manuscript, and provided the valuable advice.

Acknowledgements

This research is supported by the National Natural Science Foundation of China, No.60871086, the Natural Science Foundation of Jiangsu Province China No.BK2008159, the CSC Scholarship, PolyU, and ARC grants.

References

  • Boas DA, Brooks DH, Miller EL, DiMarzio CA, Kilmer M, Gaudette RJ, Zhang Q. Imaging the body with diffuse optical tomography. IEEE Signal Processing Mag. 2001;18:57–75. doi: 10.1109/79.962278. [Cross Ref]
  • Ntziachristos V. Fluorescence Molecular Imaging. Annual Review of Biomedical Engineering. 2006;8:1–33. doi: 10.1146/annurev.bioeng.8.061505.095831. [PubMed] [Cross Ref]
  • Ntziachristos V, Ripoll J, Wang LHV, Weissleder R. Looking and listening to light: the evolution of whole-body photonic imaging. Nat Biotechnol. 2005;23:313–320. doi: 10.1038/nbt1074. [PubMed] [Cross Ref]
  • Milstein AB, Oh S, Webb KJ, Bouman CA. Direct reconstruction of kinetic parameter images in fluorescence optical diffusion tomography. IEEE International Symposium on Biomedical Imaging. 2004. pp. 1107–1110.
  • Reynolds JS, Troy TL, Mayer RH, Thompson AB, Waters DJ, Cornell KK, Snyder PW, Sevick-Muraca EM. Imaging of spontaneous canine mammary tumors using fluorescent contrast agents. Photochem Photobiol. 1999;70:87–94. doi: 10.1111/j.1751-1097.1999.tb01953.x. [PubMed] [Cross Ref]
  • Fedele F, Eppstein MJ, Laible JP, Godavarty A, Sevick-Muraca EM. Fluorescence photon migration by the boundary element method. J Comput Phys. 2005;210:109–132. doi: 10.1016/j.jcp.2005.04.003. [Cross Ref]
  • Axelsson O. Iterative Solution Methods. Cambridge University Press; 1994.
  • Tarvainen T, Vauhkonen M, Arridge SR. Gauss Newton reconstruction method for optical tomography using the finite element solution of the radiative transfer equation. J Quant Spectrosc Radiat Transfer. 2008;109:2767–2778. doi: 10.1016/j.jqsrt.2008.08.006. [Cross Ref]
  • Haynsworth E. On the Schur complement. Basel Math Notes. 1968;20
  • Zhao B, Wang H, Chen X, Shi X, Yang W. Linearized solution to electrical impedance tomography based on the Schur conjugate gradient method. Meas Sci Technol. 2007;18:3373–3383. doi: 10.1088/0957-0233/18/11/017. [Cross Ref]
  • Silva AD, Dinten JM, Peltié P, Coll JL, Rizo P. In vivo fluorescence molecular optical imaging: from small animal towards clinical applications. 16th Mediterranean Conference on Control and Automation, Congress Centre, Ajaccio, France. 2008. pp. 1335–1340. full_text.
  • Gibson AP, Hebden JC, Arridge SR. Recent advances in diffuse optical imaging. Phys Med Biol. 2005;50:R1–R43. doi: 10.1088/0031-9155/50/4/R01. [PubMed] [Cross Ref]
  • Lin Q. Basic text book of numerical solution method for differential equation. 2. Vol. 2003. Science Press;
  • Paithankar DY, Chen AU, Pogue BW, Patterson MS, Sevick-Muraca EM. Imaging of fluorescent yield and lifetime from multiply scattered light reemitted from random media. Appl Opt. 1997;36:2260–2272. doi: 10.1364/AO.36.002260. [PubMed] [Cross Ref]
  • Arridge SR, Schweiger M. In: Computational Radiology and Imaging: Therapy and Diagnosis. Borgers C, Natterer F, editor. Vol. 110. IMA Volumes in Mathematics and Its Applications, Springer-Verlag; 1998. A general framework for iterative reconstruction algorithms in optical tomography using a finite element method; pp. 45–70.
  • Magnus JR, Neudecker H. Wiley Series in Probability and Statistics. New York: Wiley; 1988. Matrix Differential Calculus with Applications in Statistics and Econometrics.
  • Kaipio J, Somersalo E. Statistical and Computational Inverse Problems. Springer; 2005.
  • Roy R, Godavarty A, Sevick-Muraca EM. Fluorescence-Enhanced Optical Tomography Using Referenced Measurements of Heterogeneous Media. IEEE Trans Med Imaging. 2003;22:824–836. doi: 10.1109/TMI.2003.815072. [PubMed] [Cross Ref]
  • Tikhonov AN, Arsenin VY. Solution of Ill-Posed Problems. Washington, DC: Winston; 1977.
  • Pogue BW, Patterson MS, Farrell TJ. In: In Optical Tomography, Photon Migration, and Spectroscopy of Tissue and Model Media: Theory, Human Studies, and Instrumentation. Proc SPIE 2389. Chance B, Alfano RR, editor. 1995. Forward and inverse calculations for 3-D frequency-domain diffuse optical tomography; pp. 328–339.
  • Pogue BW, McBride TO, Prewitt J, Osterberg UL, Paulsen KD. Spatially variant regularization improves diffuse optical tomography. Appl Opt. 1999;38:2950–2961. doi: 10.1364/AO.38.002950. [PubMed] [Cross Ref]
  • Jennings A. Matrix Computation. Wiley; 1992.
  • Sarkar TK. On the Application of the Generalized BiConjugate Gradient Method. J Electromagnetic Waves and Applications. 1987;1:223–242. doi: 10.1163/156939387X00036. [Cross Ref]
  • Liu SW, Sung JC, Lin MS. Scattering of antiplane impact waves by cracks in a layered elastic solid. Computer Methods in Applied Mechanics and Engineering. 1998;156:335–349. doi: 10.1016/S0045-7825(97)00217-X. [Cross Ref]
  • Langtangen HP, Tveito A. A numerical comparison of conjugate gradient-like methods. Comm Appl Numer Methods. 1988;4:793–798. doi: 10.1002/cnm.1630040613. [Cross Ref]
  • Jacobsen M. Master Thesis. Technical University of Denmark, Denmark; 2000. Two-grid iterative methods for ill-posed problems.
  • Colton D, Kress R. Inverse Acoustic and Electromagnetic Scattering Theory. Springer; 1998.

Articles from BioMedical Engineering OnLine are provided here courtesy of BioMed Central