- Open Access
Image reconstruction of fluorescent molecular tomography based on the tree structured Schur complement decomposition
© Zou et al; licensee BioMed Central Ltd. 2010
- Received: 31 December 2009
- Accepted: 20 May 2010
- Published: 20 May 2010
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.
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.
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.
The methods proposed in this paper can significantly improve the reconstructed image quality of FMT and accelerate the reconstruction process.
- Inverse Problem
- Mean Square Error
- Singular Value Decomposition
- Regularization Parameter
- Tikhonov Regularization
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 . 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 . 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 .
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  and Gauss-Newton (GN) method , 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 . 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  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.
Forward Model and Finite Element Formulation
where the subscript x denotes the excitation wavelength; ∇ is the gradient operator; S x (W/cm3) is the excitation light source; Φ x (W/cm2), D x (cm), and k x (cm-1) represent the photon fluence, the diffusion coefficient, and the decay coefficient, respectively; Ω denotes the bounded domain of reconstruction.
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; 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.
where n is a vector normal to the boundary ∂Ω, b x,m is the Robin boundary coefficient.
with Ω h and Γ h being the bounded domain and its boundary, respectively.
Inverse Process of FMT
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.
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 .
where λ is a regularization parameter, which can be determined by the Morozov discrepancy principle , I is an identity matrix.
Adaptive Regularization Scheme
The problem of image reconstruction for FMT is ill-posed . The Tikhonov regularization technique, as mentioned above, is one of the major methods to reduce the ill-posedness of the problem . 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 . 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 . 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.
with Δ x i (i = 1, 2, ..., N) being the component of the vector Δ x. It can be easily seen that each node p i (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 = λ.
where r i is the position of node p i , r s and r m respectively denote the positions of the source and detector closest to the node p i , c 1 and c 2 are two positive parameters determined empirically in our paper.
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.
where k = J T J + λ and b = J T Δ y.
To solve the inverse problem of FMT in the Schur complement system, the solution space R n 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 R n . The basis of the m-dimensional coarse subspace U is formed by the columns of Γ ∈ R n × m and the columns of Ψ ∈ R n×(n-m) form the basis of the (n - m) dimensional subspace V.
where u and v are the projections of Δ x on the subspaces U and V, respectively.
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.
where , , , and , while the two components on the right-hand side (RHS) of equation (26) are , and . 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.
It can be found that the condition number of matrix S (i+1,2j) is smaller than that of matrix S (i,j) . Hence, solving the inverse problem in the Schur complement system at the i+1th level will be more efficient than solving it at the i th level. We herein solve equation (29) using the biconjugate gradient method . Its advantage is that it does not square the condition number of the original equations . 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 . The algorithm for solving equation (29) can be summarized as follows
Input an initial guess Δ x (i+1,2j)0;
Initialize d 0 = f 0 = r 0 = p 0 ← b (i+1,2j) - S (i+1,2j) Δ x (i+1,2j)0;
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.
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 . 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 i th 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 i th level. The whole reconstruction algorithm is summarized as follows
Set x 0 to an initial guess;
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 i th level for i =1,..., L and j =1,..., 2 i , the subspaces at the i th level are formed by the columns of Γ (i,j) and Ψ (i,j), respectively;
Set i = L;
For j = 1,..., 2 i 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;
While i ≥ 0
For j = 1,..., 2 i 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 i th level;
i = i - 1;
x 0 ← x 0 + Δ x (0,1);
If ||Δ x (0,1)|| > ε
go to 2;
x ← x 0, output x.
As mentioned before, the Schur complement system has a smaller condition number than that of the system from which it is constructed . 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.
Optical and fluorescent properties of one-object phantom
Optical and fluorescent properties of two-object phantom
Comparison of performance of methods
Adaptive regularization scheme
2.973 × 10-4
2.860 × 10-4
Tikhonov regularization method
5.352 × 10-4
4.892 × 10-4
Computation time of FMT image reconstruction for 30 Iterations
Performance comparison of reconstruction methods for 3D case
Computation time (s)
3.629 × 10-3
1.241 × 10-3
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.
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.
- 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. 10.1109/79.962278View ArticleGoogle Scholar
- Ntziachristos V: Fluorescence Molecular Imaging. Annual Review of Biomedical Engineering 2006, 8: 1–33. 10.1146/annurev.bioeng.8.061505.095831View ArticleGoogle Scholar
- 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. 10.1038/nbt1074View ArticleGoogle Scholar
- 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, 1107–1110.Google Scholar
- 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. 10.1111/j.1751-1097.1999.tb01953.xView ArticleGoogle Scholar
- 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. 10.1016/j.jcp.2005.04.003View ArticleGoogle Scholar
- Axelsson O: Iterative Solution Methods. Cambridge University Press; 1994.View ArticleGoogle Scholar
- 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. 10.1016/j.jqsrt.2008.08.006View ArticleGoogle Scholar
- Haynsworth E: On the Schur complement. Basel Math Notes 1968., 20: Google Scholar
- 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. 10.1088/0957-0233/18/11/017View ArticleGoogle Scholar
- 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, 1335–1340. full_textGoogle Scholar
- Gibson AP, Hebden JC, Arridge SR: Recent advances in diffuse optical imaging. Phys Med Biol 2005, 50: R1-R43. 10.1088/0031-9155/50/4/R01View ArticleGoogle Scholar
- Lin Q: Basic text book of numerical solution method for differential equation. Volume 2003. 2nd edition. Science Press;Google Scholar
- 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. 10.1364/AO.36.002260View ArticleGoogle Scholar
- Arridge SR, Schweiger M: A general framework for iterative reconstruction algorithms in optical tomography using a finite element method. In Computational Radiology and Imaging: Therapy and Diagnosis. Volume 110. Edited by: Borgers C, Natterer F. IMA Volumes in Mathematics and Its Applications, Springer-Verlag; 1998:45–70.View ArticleGoogle Scholar
- Magnus JR, Neudecker H: Matrix Differential Calculus with Applications in Statistics and Econometrics. In Wiley Series in Probability and Statistics. New York: Wiley; 1988.Google Scholar
- Kaipio J, Somersalo E: Statistical and Computational Inverse Problems. Springer; 2005.Google Scholar
- 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. 10.1109/TMI.2003.815072View ArticleGoogle Scholar
- Tikhonov AN, Arsenin VY: Solution of Ill-Posed Problems. Washington, DC: Winston; 1977.Google Scholar
- Pogue BW, Patterson MS, Farrell TJ: Forward and inverse calculations for 3-D frequency-domain diffuse optical tomography. In In Optical Tomography, Photon Migration, and Spectroscopy of Tissue and Model Media: Theory, Human Studies, and Instrumentation. Proc SPIE 2389 Edited by: Chance B, Alfano RR. 1995, 328–339.View ArticleGoogle Scholar
- Pogue BW, McBride TO, Prewitt J, Osterberg UL, Paulsen KD: Spatially variant regularization improves diffuse optical tomography. Appl Opt 1999, 38: 2950–2961. 10.1364/AO.38.002950View ArticleGoogle Scholar
- Jennings A: Matrix Computation. Wiley; 1992.Google Scholar
- Sarkar TK: On the Application of the Generalized BiConjugate Gradient Method. J Electromagnetic Waves and Applications 1987, 1: 223–242. 10.1163/156939387X00036View ArticleGoogle Scholar
- 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. 10.1016/S0045-7825(97)00217-XView ArticleGoogle Scholar
- Langtangen HP, Tveito A: A numerical comparison of conjugate gradient-like methods. Comm Appl Numer Methods 1988, 4: 793–798. 10.1002/cnm.1630040613MathSciNetView ArticleGoogle Scholar
- Jacobsen M: Two-grid iterative methods for ill-posed problems. In Master Thesis. Technical University of Denmark, Denmark; 2000.Google Scholar
- Colton D, Kress R: Inverse Acoustic and Electromagnetic Scattering Theory. Springer; 1998.View ArticleGoogle Scholar
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.