Research  Open  Published:
Compressedsensingbased fluorescence molecular tomographic image reconstruction with grouped sources
BioMedical Engineering OnLinevolume 13, Article number: 119 (2014)
Abstract
Background
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.
Methods
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.
Results
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.
Conclusions
The proposed algorithm can significantly improve the speed and accuracy of the reconstruction of FMT. Furthermore, our compressedsensingbased method can reduce the number of measurements.
Background
Optical molecular imaging has important applications in many fields such as molecular biology and clinical diagnosis, where diseasespecific tracers are combined with the optical methods to detect and localize the abnormalities at their molecular stage [1]. 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 highthroughput capability [2]. 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 [7].
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 [8]. 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 undersampled measurements [9, 10]. Therefore, the CSbased 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 [11]. In [12], a modelorder 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 CSbased 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.
Methods
Forward model
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 [13]. 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 [14].
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; $i=\sqrt{1}$; Φ_{ 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; ${\mu}_{\mathit{sx},m}^{\prime}\left(\mathbf{r}\right)$ are the isotropic scattering coefficients, and c _{ n } = c/n represents the velocity of light in tissue.
The numerical solutions for the excitation and emission density distributions in equations (1) and (2) are obtained using the Robin type boundary conditions, which take the form as follows:
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 [15].
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 [16]. 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 ${\mathrm{\Phi}}_{x,m}={\displaystyle \sum _{i=1}^{N}{\mathbf{\Phi}}_{\mathit{xi},\mathit{mi}}{\psi}_{i}}$ [17]. Hence, making use of finite element discretization of equations (1) and (2), we can yield the matrix equation as shown below:
Inverse problem
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 [18], 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.
Compressed sensing
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 [19]. A discrete signal x ∈ ℝ^{N} is said to be rsparse, if x contains r nonzero entries (r < <N). A discrete signal x ∈ ℝ^{N} is said to be rcompressible, 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.
Substituting equation (10) into equation (11), we have:
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 1norm.
The unknown signal x can be ultimately recovered by equation (10).
Compressedsensingbased 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 [20]. Therefore, sparsity can be preassumed 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:
According to equation (13), the sparse solution $\mathrm{\Delta}\tilde{\mathbf{x}}$ in equation (15) can be solved by:
min ${\mathrm{\Delta}\tilde{\mathbf{x}}}_{1}$ s.t.:
To implement CSbased reconstruction of FMT, the reconstruction model is formulated by:
where E is the objective function, _{2} denotes l 2norm, λ denotes the regularization parameter [21], 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 conjugategradient method [22]. After obtaining the sparse solution $\mathrm{\Delta}\tilde{\mathbf{x}}$, 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 CSbased 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 (N^{3}) + O (N^{3}). 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 $\frac{S}{2}\times \left(N+1\right)$ 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

(1)
Set x = x _{0} with x _{0} being an initial guess;
Set m = 0;

(2)
if (m%2==0) then
Compute Δy and J at x with the first group of sources;
else
Compute Δy and J at x with the second group of sources;
end if
Set m = m + 1;

(3)
Solve equation (17) with the nonlinear conjugategradient method;
Recover Δx using equation (14);

(4)
Update x with x = x + Δx;
Calculate the corresponding objective function J (x) by equation (18);

(5)
if (J(x) < δ) then
Stop and output x;
else
Return to step (2);
end if
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 signaltonoise 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 finiteelement 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 onetarget 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 twotarget 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.
Conclusions
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 CSbased 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 CSbased reconstruction algorithm. This can be expected to have a significant impact on the clinical applications of FMT.
Abbreviations
 FMT:

Fluorescence molecular tomography
 CS:

Compressed sensing
 RTE:

Radiative transfer equation
 FEM:

Finite element method
 MSE:

Mean square error.
References
 1.
Balas C: Review of biomedical optical imaging—a powerful, noninvasive, nonionizing technology for improving in vivo diagnosis. Meas Sci Technol 2009, 20: 1–12.
 2.
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/00319155/58/2/351
 3.
Ntziachristos V: Fluorescence molecular imaging. Annu Rev Biomed Eng 2006, 8: 1–33. 10.1146/annurev.bioeng.8.061505.095831
 4.
Ntziachristos V, Ripoll J, Wang L, Weissleder R: Looking and listening to light: the evolution of wholebody photonic imaging. Nat Biotechnol 2005, 23: 313–320. 10.1038/nbt1074
 5.
Darne C, Lu Y, SevickMuraca EM: Small animal fluorescence and bioluminescence tomography: a review of approaches, algorithms and technology update. Phys Med Biol 2014, 59: R1R64. 10.1088/00319155/59/1/R1
 6.
Hassan M, Klaunberg BA: Biomedical applications of fluorescence imaging in vivo. Comp Med 2004, 54: 635–644.
 7.
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
 8.
Ntziachristos V, Weissleder R: Experimental threedimensional fluorescence reconstruction of diffuse media by use of a normalized Born approximation. Opt Lett 2001, 26: 893–895. 10.1364/OL.26.000893
 9.
Donoho DL: Compressed sensing. IEEE Trans Inform Theory 2006, 52: 1289–1306.
 10.
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
 11.
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, SpringerVerlagBorgers C, Natterer F
 12.
Zhai Y, Cummer SA: Fast tomographic reconstruction strategy for diffuse optical tomography. Opt Express 2009, 17: 5285–5297. 10.1364/OE.17.005285
 13.
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
 14.
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
 15.
Fedele F, Eppstein MJ, Laible JP, Godavarty A, SevickMuraca EM: Fluorescence photon migration by the boundary element method. J Comput Phys 2005, 210: 1–24. 10.1016/j.jcp.2005.03.031
 16.
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/00319155/50/23/001
 17.
Arridge SR, Schotland JC: Optical tomography: forward and inverse problems. Inverse Probl 2009, 25: 1–59.
 18.
Gibson AP, Hebden JC, Arridge SR: Recent advances in diffuse optical imaging. Phys Med Biol 2005, 50: R1R43. 10.1088/00319155/50/4/R01
 19.
Donoho DL: For most large underdetermined systems of linear equations the minimal l1norm solution is also the sparsest solution. Commun Pur Appl Math 2006, 59: 797–829. 10.1002/cpa.20132
 20.
Kak AC, Slaney M: Principles of Computerized Tomographic Imaging. New York: IEEE Press; 1999.
 21.
Roy R, Godavarty A, SevickMuraca EM: Fluorescenceenhanced optical tomography using referenced measurements of heterogeneous media. IEEE Trans Med Imaging 2003, 22: 824–836. 10.1109/TMI.2003.815072
 22.
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
Acknowledgements
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.
Author information
Additional information
Competing interests
The authors declare that they have no competing interests.
Authors’ contributions
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.
Authors’ original submitted files for images
Rights and permissions
About this article
Received
Accepted
Published
DOI
Keywords
 Fluorescence molecular tomography
 Inverse problem
 Grouped sources