Research | Open | Published:
Compressed-sensing-based fluorescence molecular tomographic image reconstruction with grouped sources
BioMedical Engineering OnLinevolume 13, Article number: 119 (2014)
Although the quality of reconstructed results can be improved with the increment of the number of measurements, the scale of the matrices involved in the reconstruction of fluorescence molecular tomography (FMT) will become larger, which leads to the poor efficiency of the process of tomographic image reconstruction. In this paper, we proposed a new method for image reconstruction of FMT based on compressed sensing, in which a scheme of grouped sources is incorporated.
The forward equations are implemented using the finite element method (FEM). The reconstruction model is formulated under the framework of compressed sensing theory. The regularization term and the total variation penalty are incorporated in the objective function. During the reconstruction of FMT, the sources are divided into two groups for iteration in turn. One group of sources is employed in the first iteration of inverse problem, and the other group is employed in the next iteration.
Simulation results demonstrate that the computation time and mean square error (MSE) of the reconstruction with our algorithm are less than those with the traditional method. The proposed algorithm can reconstruct the target with enhanced contrast and more accurate shape.
The proposed algorithm can significantly improve the speed and accuracy of the reconstruction of FMT. Furthermore, our compressed-sensing-based method can reduce the number of measurements.
Optical molecular imaging has important applications in many fields such as molecular biology and clinical diagnosis, where disease-specific tracers are combined with the optical methods to detect and localize the abnormalities at their molecular stage . Optical molecular imaging techniques further impart the potential ability to prevent and treat the lethal diseases. As one of the emerging optical molecular imaging techniques, fluorescence molecular tomography (FMT) plays a significant role and attracts considerable new interest due to its high sensitivity, low cost, and high-throughput capability . FMT has been developed as a tomographic method to yield a robust modality for fluorescent reporters [3, 4]. FMT can be utilized in the studies of drug discovery, tumor diagnosis and therapy assessment. In this imaging modality, an external source of excitation light is needed to irradiate the tissues and then the injected fluorescent markers absorbs the incident light [5, 6]. Upon releasing the energy, the fluorophore emits the light at a longer wavelength and the fluorescent molecule decays to its ground state. The emission light is measured by measurement devices at the surface of the tissue. Monitoring the fluorescent markers is able to collect a large amount of functional information. The images of the absorption coefficient, fluorescent lifetime and the fluorescent yield can be reconstructed from the measured data and a mathematical model of light propagation .
Two processes are required to reconstruct the image of FMT. First, a forward model is utilized to map the parameters to the measurable data. Second, an inverse problem is used to calculate the spatial distribution of the optical and fluorescent properties when the measured data and sources are given. Ntziachristos et al. presented a normalized Born approximation to reconstruct the distribution of fluorochromes with different concentrations embedded in the media . In image reconstruction of FMT, the quality of reconstructed results can be improved with the increasing number of measurement data. However, the scale of the matrices involved in the reconstruction will become larger, which may slow down the process of solving the tomographic inverse problem. Therefore, a new method based on compressed sensing (CS) is proposed in this paper to tackle such a problem and accelerate the reconstruction process. CS is an innovative information theory, which is able to recover sparse signals with the under-sampled measurements [9, 10]. Therefore, the CS-based method of FMT is capable of reconstructing the tomographic images with high speed and low cost from less number of measurements. This will provide potential benefit for biomedical applications. Besides the measurements, source is another important factor to determine the efficiency of the reconstruction. Since the fluorescence molecule is excited to emit the fluorescence light by the source, the position and number of sources can greatly affect the reconstructed results. The quality of reconstruction can be improved with the increased number of sources, which may lead to larger scale of matrix and higher computational requirements involved in the reconstruction process . In , a model-order reduction approach is proposed for reduction of the system complexity. However, the transformation matrix is needed to be constructed with the Wilson–Yuan–Dickens basis vectors or the Lanczos basis vectors in the Krylov subspace. Therefore, we propose to implement the CS-based reconstruction with the grouped sources, which can improve the efficiency of the reconstruction process. Simulation results demonstrate that our proposed method can significantly speed up the reconstruction process with high reconstruction accuracy.
In general, light propagation in the near infrared spectral window is well modelled by the radiative transfer equation (RTE). In order to reduce the computational complexity, the diffusion equation is employed instead. The diffusion equation is the simplest nontrivial approximation that results from the P 1 approximation to RTE . The diffusion equation is widely utilized in the modelling of the light propagation in tissue. Mathematically, the following diffusion equation depicts the process of excitation light propagation in the frequency domain .
In equation (1), the subscript x denotes the properties at the excitation wavelength; D x (r) and k x (r,ω) represent the optical diffusion coefficient and decay coefficient, respectively; Φ x (r, ω) and Q x (r,ω) denote the photon density and excitation light source, respectively; ∇ represents the gradient operator.
The other diffusion equation, which depicts the generation and propagation of fluorescent light, is of the form:
In equation (2), the subscript m denotes the parameters at the emission wavelength; ; Φ m (r,ω) represents the photon density; Other quantities that are included in equation (2) are the fluorescence lifetime τ(r), the fluorescence quantum efficiency ϕ, the diffusion coefficients D m (r), and the decay coefficients k m (r,ω). Further, k x,m (r,ω) and D x,m (r) can be obtained below:
where μ ax,mf (r) are the absorption coefficients due to fluorophore; μ ax,mi (r) are the absorption coefficients due to nonfluorescing chromophore; are the isotropic scattering coefficients, and c n = c/n represents the velocity of light in tissue.
where n(r) is outward facing normal vector for the boundary ∂Ω, b x,m (r) are the Robin boundary coefficients at excitation and emission wavelengths, which depend on the optical refractive index mismatch at the boundary .
The forward equations can be computationally implemented using the analytical methods or the finite element method (FEM). Basically, analytical methods may lead to inaccurate results due to simplified assumptions of properties or geometry . The most significant superiority of FEM is versatility. Thus FEM is applicable to inhomogeneous property distributions and complex geometries. Herein we solve the forward equations with FEM. To facilitate a finite element solution of the forward equations, the reconstruction domain Ω is discretized into P elements with N vertex nodes. By means of the basis functions ψ i (i = 1,2,…,N), the solution Φ x,m can be represented by . Hence, making use of finite element discretization of equations (1) and (2), we can yield the matrix equation as shown below:
The inverse problem is to recover the distribution of tissue parameters x with a series of boundary measurements y as well as several sources. Let us denote the forward operator by F. Thus, the inverse problem can be represented by:
According to the Taylor series expansion , the linear inverse problem can be expressed by the following matrix form:
where Δx denotes the changes in the optical properties, Δy denotes the residual data between the predicted and measured data, and J denotes the Jacobian matrix.
One important condition of CS is sparsity, which means that either the signal itself or its representation in an appropriate basis is sparse or compressible . A discrete signal x ∈ ℝN is said to be r-sparse, if x contains r nonzero entries (r < <N). A discrete signal x ∈ ℝN is said to be r-compressible, if x can be represented by r large coefficients and other coefficients are small in magnitude. Basically, most medical images have sparse representations in some orthonormal basis (e.g., wavelet, Fourier). That is:
where x is the original image vector, θ is the transformed vector, and Ψ denotes the orthonormal basis.
A central idea in the CS theory is to acquire the signal x with the projection measurement vector y ∈ ℝM as follows:
where Φ is a measurement matrix with a size of M × N.
According to the CS theory, the sparse solution θ in equation (12) can be obtained by solving the constrained optimization problem as follows:
min ||θ||1 subject to:
where ||||1 denotes l 1-norm.
The unknown signal x can be ultimately recovered by equation (10).
Compressed-sensing-based image reconstruction of FMT with grouped sources
For the linearized reconstruction, the perturbation from a homogeneous background or reference medium is relatively small in volume . Therefore, sparsity can be pre-assumed for the linearized tomographic imaging problem. Typically in biological media, sparsity may be true if a very good reference medium is available. Thus, a transformation basis such as the Fourier basis can be utilized to reinforce the sparsity of the solution to equation (9) (i.e., Δx) as follows
Based on the discussion in the above two sections, the problem of image reconstruction of FMT under the framework of CS theory can be expressed by:
To implement CS-based reconstruction of FMT, the reconstruction model is formulated by:
where E is the objective function, ||||2 denotes l 2-norm, λ denotes the regularization parameter , and α represents the parameter determining the sparsity. Basically, it will be useful to incorporate a total variation penalty into equation (17), which can be solved via the nonlinear conjugate-gradient method . After obtaining the sparse solution , we can recover the solution Δx in equation (9) by the transform as equation (14).
As mentioned before, the measurements are generated with the excitation source for FMT. Therefore, the source placed in different position can provide different information to reconstruction. Furthermore, the number of sources can also affect the reconstruction results. More number of sources will improve quality of reconstruction, whereas the scale of the Jacobian matrix will increase. Consequently, the computational burden of the whole reconstruction may increase. Thus, in order to tackle such a problem, a scheme of grouped sources is proposed to improve a basic CS-based reconstruction of FMT. In this scheme, the sources are divided into two groups. One group of sources is employed in the first iteration of inverse problem, and the other group is employed in the next iteration. Two groups of sources can provide more information for image reconstruction than only one group of sources, which may improve the quality of reconstruction. More importantly, in this way, the number of rows in the Jacobian matrix can be reduced to only one half of that in the Jacobian matrix with all sources. Suppose the number of vertex nodes, sources, and measurements are N, S, and M, respectively. The computational complexity of solving the forward equations of FMT will be O (N3) + O (N3). Basically, the direct calculation of the Jacobian matrix with all sources needs S × (N + 1) FEM forward calculations. With the scheme of grouped sources, the calculation of the Jacobian matrix is reduced to only FEM forward calculations. Thus, the computational requirements for the Jacobian matrix can obviously be reduced using the proposed algorithm. This is helpful for reducing the computational burden of the tomographic inverse problem. In addition, the iterative reconstructed results from one group of sources will offer a good initial guess for the next iteration of reconstruction from the other group of sources. Figure 1 shows the illustration of the grouped sources strategy. In Figure 1(a), the sources without being grouped are illustrated with black rectangles. Figure 1(b) and (c) display the two groups of sources with blue rectangles and red rectangles, respectively.
In fact, the solution to the inverse problem can be achieved by minimizing the following objective function:
where J (x) denotes the objective function measuring the discrepancy between the measured data y and predicted data F(x) with regards to the forward model. Therefore, the overall reconstruction algorithm consists of following steps
Set x = x 0 with x 0 being an initial guess;
Set m = 0;
if (m%2==0) then
Compute Δy and J at x with the first group of sources;
Compute Δy and J at x with the second group of sources;
Set m = m + 1;
Solve equation (17) with the nonlinear conjugate-gradient method;
Recover Δx using equation (14);
Update x with x = x + Δx;
Calculate the corresponding objective function J (x) by equation (18);
if (J(x) < δ) then
Stop and output x;
Return to step (2);
Results and discussion
In this section, two simulation experiments were performed to test the efficacy of the proposed algorithm: image reconstruction of a single fluorescence target, and two closely spaced fluorescence targets. The synthetic measurement data are generated from the diffusion equations as equations (1) and (2). In addition, random Gaussian noise with a signal-to-noise ratio of 10 dB is applied to these synthetic measurement data to simulate measurement error. The initial guess x 0 in this study is set to 5 mm-1.
Single fluorescent target
Figure 2 displays the simulated phantom with a single fluorescent target. The phantom has eight sources and thirty detectors for the measurement.
Implementation of a finite-element simulation needs an appropriate mesh upon which to simulate the reconstruction results. In order to reduce the computational burden without greatly reducing the image resolution, the mesh used for reconstruction is refined adaptively using the a priori image displayed in Figure 3. First, the reconstruction domain is uniformly discretized. Then, the uniform mesh is refined for the areas with large variations of the pixel values in the a priori image. Whether the mesh needs to be refined is judged by following formula
where X denotes the pixel value in the prior image, D denotes the variation of pixel values in the triangle, and E denotes the expectation operator. The corresponding triangle will be refined when the variation is larger than the assumed threshold. In such a way, the adaptively refined mesh can be generated. As shown in Figure 4, the adaptively refined mesh has 122 nodes and 212 elements. The optical properties of the target and the background with regards to Figure 2 are given in Table 1.
The reconstruction results of μ axf for the single target case from the traditional method that without using CS are summarized in Figure 5. The tomographic reconstruction with 30 measurements is shown in Figure 5(a), and the reconstruction with 18 measurements is depicted in Figure 5(b). It is seen that the quality of image reconstruction for the location of target as well as the contrast can be improved with the increasing number of the measurements, which may lead to higher computational requirements.The reconstructed absorption distributions for the single target case using the different algorithms are displayed in Figure 6, where Figure 6(a) shows the traditional reconstruction result with 30 measurements, and Figure 6(b) depicts the reconstruction result based on the proposed method with 15 measurements. As can be seen, the proposed algorithm can reconstruct the target with enhanced contrast and more accurate shape.
The accuracy of the reconstructions is analyzed quantitatively by computing the mean square error (MSE) from the distribution of the optical properties for the reconstructed data sets.
where N represents the total number of nodes in the domain. The superscript cal denotes the reconstructed results, and act denotes the actual distribution of the optical properties.
The quantitative performance of reconstruction for the one-target phantom in terms of the two metrics as computation time and MSE is given in Table 2 to further evaluate the reconstruction quality. Examining this table, we observe that our proposed algorithm can decrease the computation time of reconstruction as compared with the traditional method. Furthermore, Table 2 suggests that the proposed method can provide relatively high reconstruction quality with low computational requirements.
Dual fluorescent targets
The simulated phantom with dual fluorescent targets of different shapes is illustrated in Figure 7. The source–detector configurations of the two-target phantom are the same as the single target case.
Figure 8 depicts the prior image, which is introduced to refine the mesh for reconstruction. Thus, the adaptively refined mesh has 148 nodes and 264 elements as shown in Figure 9. Table 3 outlines the optical properties of the targets and the background with regards to Figure 7.
In Figure 10, the reconstructed images of μ axf for the dual targets case with 30 measurements (see Figure 10(a)) and 18 measurements (see Figure 10(b)) are shown. The reconstructed results in Figure 10 are obtained based on the traditional method that without using CS. We can also observe that the higher accuracy of locations of the dual targets and the contrast from the reconstruction result can be achieved with the increasing measurements. On the other hand, the computational burden of reconstruction will become greater with more measurements.Figure 11 depict the reconstructed absorption distributions for the dual targets case using the different algorithms. The traditional reconstruction result with 30 measurements and that based on our method with 15 measurements are shown in Figure 11(a) and (b), respectively. We can observe that the shape and contour of the reconstructed targets with the proposed algorithm are better than those with the traditional method. Note also that the contrast can be improved with the proposed algorithm.
Table 4 lists the performance of reconstructions in terms of the computation time and MSE for the phantom with dual targets to assess the reconstruction quality. As can clearly be seen, the computation time and MSE of the reconstruction with our algorithm are less than those with the traditional method. In other words, for the dual targets case, our proposed algorithm can also improve the reconstruction speed with high accuracy.
In this work, we developed an innovative method based on CS for image reconstruction of FMT. A scheme of grouped sources is incorporated in the reconstruction process. In comparison to traditional reconstruction approach, the CS-based reconstruction algorithm has demonstrated significant advantages as faster speed and high accuracy. Furthermore, the cost and the amount of measurements for image reconstruction can be reduced with the CS-based reconstruction algorithm. This can be expected to have a significant impact on the clinical applications of FMT.
Fluorescence molecular tomography
Radiative transfer equation
Finite element method
Mean square error.
Balas C: Review of biomedical optical imaging—a powerful, non-invasive, non-ionizing technology for improving in vivo diagnosis. Meas Sci Technol 2009, 20: 1–12.
Zhang G, Cao X, Zhang B, Liu F, Luo J, Bai J: MAP estimation with structural priors for fluorescence molecular tomography. Phys Med Biol 2013, 58: 351–372. 10.1088/0031-9155/58/2/351
Ntziachristos V: Fluorescence molecular imaging. Annu Rev Biomed Eng 2006, 8: 1–33. 10.1146/annurev.bioeng.8.061505.095831
Ntziachristos V, Ripoll J, Wang L, Weissleder R: Looking and listening to light: the evolution of whole-body photonic imaging. Nat Biotechnol 2005, 23: 313–320. 10.1038/nbt1074
Darne C, Lu Y, Sevick-Muraca EM: Small animal fluorescence and bioluminescence tomography: a review of approaches, algorithms and technology update. Phys Med Biol 2014, 59: R1-R64. 10.1088/0031-9155/59/1/R1
Hassan M, Klaunberg BA: Biomedical applications of fluorescence imaging in vivo. Comp Med 2004, 54: 635–644.
Milstein AB, Webb KJ, Bouman CA: Estimation of kinetic model parameters in fluorescence optical diffusion tomography. J Opt Soc Am A 2005, 22: 1357–1368. 10.1364/JOSAA.22.001357
Ntziachristos V, Weissleder R: Experimental three-dimensional fluorescence reconstruction of diffuse media by use of a normalized Born approximation. Opt Lett 2001, 26: 893–895. 10.1364/OL.26.000893
Donoho DL: Compressed sensing. IEEE Trans Inform Theory 2006, 52: 1289–1306.
Candes EJ, Romberg JK, Tao T: Stable signal recovery from incomplete and inaccurate measurements. Commun Pur Appl Math 2006, 59: 1207–1223. 10.1002/cpa.20124
Arridge SR, Schweiger M: A general framework for iterative reconstruction algorithms in optical tomography using a finite element method. Comput Radiol Imaging 1998, 110: 45–70. IMA Volumes in Mathematics and Its Applications, Springer-VerlagBorgers C, Natterer F
Zhai Y, Cummer SA: Fast tomographic reconstruction strategy for diffuse optical tomography. Opt Express 2009, 17: 5285–5297. 10.1364/OE.17.005285
Arridge SR, Dehghani H, Schweiger M, Okada E: The finite element model for the propagation of light in scattering media: A direct method for domains with nonscattering regions. Med Phys 2000, 27: 252–264. 10.1118/1.598868
Martelli F, Bianco SD, Ninni PD: Perturbative forward solver software for small localized fluorophores in tissue. Biomed Opt Express 2012, 3: 26–36. 10.1364/BOE.3.000026
Fedele F, Eppstein MJ, Laible JP, Godavarty A, Sevick-Muraca EM: Fluorescence photon migration by the boundary element method. J Comput Phys 2005, 210: 1–24. 10.1016/j.jcp.2005.03.031
Chaudhari AJ, Darvas F, Bading JR, Moats RA, Conti PS, Smith DJ, Cherry SR, Leahy RM: Hyperspectral and multispectral bioluminescence optical tomography for small animal imaging. Phys Med Biol 2005, 50: 5421–5441. 10.1088/0031-9155/50/23/001
Arridge SR, Schotland JC: Optical tomography: forward and inverse problems. Inverse Probl 2009, 25: 1–59.
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/R01
Donoho DL: For most large underdetermined systems of linear equations the minimal l1-norm solution is also the sparsest solution. Commun Pur Appl Math 2006, 59: 797–829. 10.1002/cpa.20132
Kak AC, Slaney M: Principles of Computerized Tomographic Imaging. New York: IEEE Press; 1999.
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.815072
Lustig M, Donoho DL, Pauly JM: Sparse MRI: the application of compressed sensing for rapid MR imaging. Magn Reson Med 2007, 58: 1182–1195. 10.1002/mrm.21391
This work was supported by Natural Science Foundation of Jiangsu Province, China under Grant No. BK20130324, Specialized Research Fund for the Doctoral Program of Higher Education (SRFDP) under Grant No. 20123201120009, and Natural Science Foundation of the Jiangsu Higher Education Institutions of China under Grant No. 12KJB510029.
The authors declare that they have no competing interests.
WZ conceived the study, implemented the algorithm, and drafted the manuscript. XYP helped to analyze the results, and revised the manuscript. Both authors read and approved the final manuscript.