Open Access
Add to BibSonomy
Abstract
We present a reconstruction method involving maximum-likelihood expectation maximization (MLEM) to model Poisson noise as applied to fluorescence molecular tomography (FMT). MLEM is initialized with the output from a sparse reconstruction-based approach, which performs truncated singular value decomposition-based preconditioning followed by fast iterative shrinkage-thresholding algorithm (FISTA) to enforce sparsity. The motivation for this approach is that sparsity information could be accounted for within the initialization, while MLEM would accurately model Poisson noise in the FMT system. Simulation experiments show the proposed method significantly improves images qualitatively and quantitatively. The method results in over 20 times faster convergence compared to uniformly initialized MLEM and improves robustness to noise compared to pure sparse reconstruction. We also theoretically justify the ability of the proposed approach to reduce noise in the background region compared to pure sparse reconstruction. Overall, these results provide strong evidence to model Poisson noise in FMT reconstruction and for application of the proposed reconstruction framework to FMT imaging.
© 2018 Optical Society of America under the terms of the OSA Open Access Publishing Agreement
1. Introduction
Fluorescence molecular tomography (FMT) is finding several applications in 3D visualization and quantification of the distribution of molecular target within biological tissue [1]. In particular, FMT has received substantial interest in small animal imaging for applications such as studying tumor physiology and for pharmaceutical research [2, 3]. In FMT imaging, fluorescence molecules are first injected into biological tissue. External illumination sources are used to excite the fluorescence molecules. The photons emitted by the excited fluorescence molecules are collected by detectors at the tissue surface. The objective in FMT is to use these surface measurements to reconstruct the 3D distribution of fluorescence molecules within the tissue.
The reconstruction problem in FMT is known to be highly ill-posed, and is sensitive to noise and modeling errors such as discretization [4,5]. Over the past two decades, various reconstruction methods for FMT have been proposed [6]. Tikhonov regularization is a popular regularization applied to FMT reconstruction problem. The regularized problem can be solved iteratively with methods such as Newton method and algebraic reconstruction technique (ART) [6, 7]. However, such regularization tends to over-smooth the reconstructed images, leading to loss of localized features during reconstruction [8]. More recently, reconstruction methods that exploit sparsity of the fluorescence distribution have been studied [5, 9–11]. In these methods, ℓ0 or ℓ1 regularization on the fluorescence distribution is applied to enforce sparsity while performing the reconstruction. These regularization problems can be solved with methods such as greedy algorithms and iterative thresholding methods [12].
In FMT systems, often the detector system is charged-couple device (CCD)-based camera or photon multiplier tube (PMT). In these systems, the noise is described by a Poisson distribution [13–15]. For this noise distribution, MLEM-based reconstruction techniques have yielded reliable results, especially in nuclear medicine imaging [16–20]. The MLEM technique has several advantages, such as accurately modeling the Poisson noise distribution in the acquired data, constraining the activity values to be non-negative without the need for a specific regularizer, and ensuring the conservation of the total number of photons across multiple iterations. In optical tomography, several studies have applied MLEM for reconstruction in bioluminescence tomography [21–23]. In [24], MLEM has also been applied for FMT reconstruction. However, the MLEM technique typically suffers from slow convergence for optical tomography modalities, with thousands of iterations and large amount of time per iteration being required [23, 25, 26]. This makes MLEM a time-consuming method and thus not very practical [21, 24]. As a result, MLEM has not been widely used for image reconstruction in optical tomography.
The performance of MLEM is influenced by different factors. An important factor being the initial estimate provided to the algorithm. Conventionally, MLEM starts with a uniform initial estimate, as we explain later. However, different initializations for MLEM yield different reconstruction results [27, 28]. In this work, we studied the use of sparse reconstruction to initialize the MLEM approach. The overall motivation for this approach is that the sparse reconstruction method would account for the sparsity of the fluorescence distribution, while the MLEM would accurately model the Poisson noise in the FMT system. However, this combined approach is also able to exploit several inherent advantages of these two techniques, as we describe below. Our method yields reliable and improved results in comparison to pure sparse reconstruction as well as uniformly initialized MLEM methods. Preliminary versions of this work have been presented previously [29, 30]. We begin by describing our method in the next section.
2. Methods
2.1. The forward model and reconstruction problem in FMT
The forward model in FMT is described by a pair of coupled equations. The first equation describes the propagation of excitation photons from source at location rs to location r in the medium and the second one describes the propagation of emitted fluorescence photons from location r to detector at location rd, where rs, r and rd are three-dimensional vectors. These coupled equations are given by:
and where ϕex(r) and ϕem(rd) are the excitation light field at location r and emission light field at detector location rd, respectively, gex(rs, r) is the Green’s function of excitation light at location r due to a source at location rs, gem(r, rd) denotes Green’s function of emission light detected by detector at location rd due to the fluorescence source at location r, x(r) is the fluorescence yield at location r, and Ω denotes object support,. If we discretize Ω into N voxels, we obtain the linear matrix equation for the forward model: where is the sensitivity matrix of the system, Φ is an M × 1 vector denoting detector measurements, x is an N × 1 vector representing unknown fluorescence yield, Ns and Nd are number of sources and detectors, respectively, and M = Ns × Nd is the total number of measurements. Due to the limited number of sources and detectors, typically M < N in FMT.Modeling the measurement noise denoted by the M-dimensional vector n, Eq. (3) becomes:
In FMT, the data collected by the detectors is corrupted by Poisson noise [13]. The reconstruction problem in FMT is to reconstruct x given sensitivity matrix G and detector measurements Φ. In the next section, we derive the MLEM-based reconstruction technique that models this noise distribution accurately.2.2. Modeling Poisson noise in the reconstruction
The likelihood function for Poisson distributed data is:
where (Gx)m and ϕm denote the mth element of the vector (Gx) and Φ, respectively. Taking the logarithm of the likelihood function yields: The first order derivative of the log-likelihood function is given by Setting yields Multiplying both side with x and replacing x with a sequence of estimates yields the fixed-point iteration: where . This is referred to the MLEM technique [27]. The MLEM technique described in Eq. (9) provides an iterative way to maximize the Poisson likelihood function.The MLEM iteration starts from an initial estimate , and the results of this technique can be influenced by its initial estimate [27]. Typically, the initial estimate is uniform, where all the elements in are assumed to be a constant [31, 32]. However, with this estimate, MLEM updates all the voxels in every iteration, increasing the computational requirements. In Eq. (9), note that will always be zero if due to the multiplicative nature of the technique. Thus, the zero elements can be excluded from during MLEM iteration. Matrix G used for MLEM iteration can be formulated with columns corresponding to non-zero elements in . This reduces the size of matrices in the reconstruction problem and accelerates the computation speed. In this context, in many FMT applications, fluorescence molecules tend to concentrate in a small target region. Thus, if we could exploit this property, we could generate a sparse initial estimate, which allows us to accelerate the MLEM technique. Such a technique would inherently exploit the sparsity-based prior information in FMT as well as model the Poisson noise in FMT accurately. Inspired by this, we developed a sparse reconstruction method and used the output from this method as the initial estimate for MLEM. In the next section, we describe the method we used to obtain sparse initial estimate of MLEM.
2.3. Sparse reconstruction and preconditioning of sensitivity matrix
To provide the sparse initial estimate for MLEM, the following minimization problem can be formulated based on Eq. (4):
While directly solving this problem is computationally complex, Eq. (10) can be approximately solved with greedy algorithms or convex relaxation techniques [12]. The theory of compressed sensing (CS) provides the conditions under which such approximate solvers are valid. Further, approaches based on singular value decomposition (SVD) can be applied to the sensitivity matrix to improve sparse reconstruction in FMT [33–36]. This technique is known as preconditioning of sensitivity matrix. Here, we follow truncated singular value decomposition (TSVD) described in [34] as the preconditioning method. First, expressing the matrix G in terms of its singular vectors and singular values using SVD, Eq. (4) becomes: where U and V are M × M and N × N unitary matrices where the columns are left-singular vectors and right-singular vectors, respectively, and Σ is a diagonal matrix where the diagonal elements are the singular values. By multiplying both sides of Eq. (11) with Σ−1UT, we could potentially use VT as the new sensitivity matrix. However, since the reconstruction problem in FMT is highly ill-posed, the inversion of small singular values contained in Σ will cause large noise amplification. To address this issue, we keep only the K largest singular values of matrix Σ and discard the rest, before performing the inversion of Σ. The corresponding columns in U and V are also discarded. This process is referred to as truncation. Then Eq. (11) becomes where the size of Ut, Σt and Vt are M × K, K × K and N × K, respectively. Since small singular values are discarded, usually K < M. Applying to both sides of Eq. (12) yields Denoting y = MΦ, and n′ = Mn, Eq. (13) can be written as We now solve Eq.(14) as a sparse reconstruction problem. More specifically, we implemented convex relaxation technique in this work. Our objective is to minimize the ℓ1 norm of the vector x. Thus the sparse reconstruction problem is posed as We applied the fast iterative shrinkage-thresholding algorithm (FISTA) for solving the minimization problem in Eq. (15) [37]. The output with this method is then input to the MLEM technique as the initial estimate. Note that results from sparse reconstruction might contain negative elements. As we explained previously, MLEM constrains the activity values to be non-negative. To enable this, the negative elements in are set to zero.2.4. Experiments
To validate the proposed method, different simulation experiments were conducted. Three different reconstruction methods were implemented for comparison, namely, (a) a pure sparsity-based reconstruction method that used TSVD in conjunction with FISTA, (b) the MLEM method with uniform initial estimate of the image (more specifically, the initial activity values in all the voxels was set to unity) and (c) the MLEM method with an initialization that was obtained using the method described in (a). We will refer to these methods as pure sparsity-based reconstruction method, uniformly initialized MLEM and sparsity-initialized MLEM, respectively.
In the first set of experiments, a 5 × 5 × 5 cm3 cubic phantom was considered, as shown in Fig. 1(a). The phantom was discretized into 20 × 20 × 20 voxels. The absorption coefficient of the phantom was set to µa = 0.05 cm−1 and the reduced scattering coefficient was set to . 20 sources and 144 detectors were placed on the side surfaces. This configuration generated 2880 measurements. Two cylindric fluorescence bars with radius of 0.375 cm and length of 2.5 cm each were inserted into the phantom. The fluorescence intensity in these bars was set to unity. The cross section of the phantom at y = 2.5 cm is shown in Fig. 1(b). The Green’s function in the forward model of FMT was computed using Monte Carlo method, where a large number of photons were simulated to generate approximately noiseless measurements [38]. The measurements were then scaled to different levels and corresponding Poisson noise was applied using a Poisson distributed pseudo random number generator. This yielded detector measurements with different signal-to-noise ratio (SNR) values.
Fig. 1 (a) The experimental setup of cube phantom. (b) Cross section at y = 2.5 cm of the simulated phantom.
To study the effect of MLEM iteration number on reconstruction performance, 1000 iterations were performed for MLEM with different initializations with the SNR initially set to 18dB, and the truncation number K set to 760. The region of interest (ROI) corresponded to the region occupied by the fluorescence bars. The rest of the region was defined as background. For quantitative study, different figures of merit were computed. Specifically, we computed absolute bias in the estimated uptake in the ROI and the background, spatial variance within the pixels in the ROI and the background, and the root mean square error (RMSE) for the entire image. The mean of the fluorescence uptake within the ROI, denoted by θROI, is defined as
where r denotes the rth voxel in the ROI, and NR is the number of voxels in the ROI. Similarly, the background mean, denoted by θB, is defined as where b denote the bth voxel in the background region, and NB is the number of voxels in the background. Then the ROI absolute bias, denoted by bROI, was computed as: where k denotes the kth noise realization, denotes the true uptake in the kth voxel in the ROI, and R is the total number of noise realizations. The background absolute bias, denoted by bB, was computed as: where denotes the true uptake in the kth voxel in the background. We also computed the spatial variance within the pixels in the ROI (denoted by ) and in the background (denoted by ) as follows: The RMSE over the entire 3D image was computed as below: where the subscript k denotes the kth noise realization and the subscript i denotes the ith voxel. In this and all the other experiments in this paper, 100 noise realizations were used to compute the various figures of merit. To study the sensitivity of our method to noise, experiments were conducted with SNR ranging from 5 dB to 40 dB, with step size of 5 dB.In the second set of experiments, we conducted simulation studies with a digital mouse phantom [39]. Three fluorescence targets were placed in the mouse brain. Two of them had a radius of 0.8 mm and the third had a radius of 1.2 mm. The optical properties of the mouse head are listed in Table 1. The whole brain was discretized into 2942 voxels. 48 sources and 51 detectors were placed at the surface of the mouse head, as shown in Fig. 2(a). The cross section of the phantom at z = 16mm is shown in Fig. 2(b).
Table 1. Optical properties of digital mouse phantom [40]
Fig. 2 (a) The experimental setup of mouse phantom. (b) Cross section of digital mouse phantom at z = 16mm.
First, 1000 iterations were performed for MLEM with uniform and sparse initialization to study the effect of iteration number on MLEM performance. The SNR was set to 18 dB. The truncation number K was set to 120. Next, quantitative performance of pure sparse reconstruction, sparsity-initialized MLEM and uniformly initialized MLEM methods at different noise levels were evaluated. The SNR value ranged from 5 dB to 40 dB, with step size of 5 dB.
The selection of truncation number plays an important role in the quality of the reconstructed image acquired from sparse reconstruction [35, 36]. For this reason, we also studied the effect of truncation number on reconstruction results of pure sparse reconstruction method and the proposed sparsity-initialized MLEM method. To compare our proposed method and pure sparse reconstruction method, we conducted experiments with different truncation number K. 450 iterations were used for MLEM with sparse initial estimate. For quantitative study, RMSE was computed as a function of the truncation number. The experiments were conducted for two noise levels, namely SNR=40 dB and SNR=20 dB.
3. Results
3.1. Uniform cube phantom
Figure 3 shows cross sections reconstructed by MLEM with different iteration numbers. For sparsity-initialized MLEM, iteration number n = 0 corresponds to the case of pure sparse reconstruction. The fluorescence intensity in all figures were normalized to the range of [0, 1]. We determined the computation time for 1000 iterations of the MLEM algorithm with different initializations. The result is shown in Table 2. Note that these results do not include the time for generating the initial estimate for MLEM. The result is shown in Table 2. It can be observed that sparsity-initialized MLEM is about 8 times faster than uniformly initialized MLEM.
Fig. 3 Cross sections at y = 2.5 cm reconstructed by MLEM with different iteration number n for the cube phantom. SNR=18dB. The reconstructed images are from MLEM with uniform initial estimate for the top row and MLEM with sparse initial estimate for the bottom row.
Table 2. Computation time required by MLEM for 1000 iterations
Figure 4 shows the quantitative results as a function of iteration number. We observe from these plots that sparsity-initialized MLEM converges at a lower number of iterations. Sparsity-initialized MLEM has lower ROI bias, background bias, background spatial variance, and image RMSE. We also computed the variance of the mean ROI and the mean background uptakes, and found that these were much lower (less than 1%) compared to the bias. Thus, we do not show these results here. From Fig. 4(f), we notice that sparsity-initialized MLEM reached its lowest RMSE after only 50 iterations, but for uniformly initialized MLEM, the lowest RMSE was obtained after 800 iterations. Based on this result, we chose 50 iterations for sparsity-initialized MLEM and 800 iterations for uniformly initialized MLEM for the first set of experiments with different SNR values. The plots of quantitative results for the different reconstruction methods at different SNR values are shown in Fig. 5. we again observe that sparsity-initialized MLEM leads to lower ROI bias, background bias, background spatial variance, and image RMSE for all noise levels.
Fig. 4 Quantitative results of different reconstruction methods as functions of iteration number for cube phantom. (a) Plot of ROI bias vs. number of iterations. (b) Plot of ROI spatial variance vs. number of iterations. (c) Plot of ROI spatial variance vs. ROI bias. (d) Plot of background bias vs. number of iterations. (e) Plot of background variance vs. number of iterations. (f) Plot of RMSE vs. number of iterations.
Fig. 5 Quantitative results of different reconstruction methods as functions of SNR for cube phantom. (a) Plot of ROI bias vs. SNR. (b) Plot of ROI variance vs. SNR. (c) Plot of background bias vs. SNR. (d) Plot of background variance vs. SNR. (e) Plot of RMSE vs. SNR.
3.2. Digital mouse phantom
Figure 6 shows quantitative performance of different reconstruction methods as a function of iteration number. We observe that sparsity-initialized MLEM achieves lower ROI bias, background bias and RMSE. It was observed that for the sparsity-initialized MLEM and uniformly initialized MLEM, 450 and 900 iterations yielded the minimum RMSE. Thus, these values were chosen for the two methods for subsequent experiments. Quantitative performance of different reconstruction methods at different noise levels is shown in Fig. 7. The sparsity-initialized MLEM method shows better performance for ROI bias, background bias and RMSE compared to the other two methods.
Fig. 6 Quantitative results of different reconstruction methods as functions of iteration number for digital mouse phantom. (a) Plot of ROI bias vs. number of iterations. (b) Plot of ROI spatial variance vs. number of iterations. (c) Plot of ROI spatial variance vs. ROI bias. (d) Plot of background bias vs. number of iterations. (e) Plot of background variance vs. number of iterations. (f) Plot of RMSE vs. number of iterations.
Fig. 7 Quantitative results of different reconstruction methods as functions of SNR for digital mouse phantom. (a) Plot of ROI bias vs. SNR. (b) Plot of ROI variance vs. SNR. (c) Plot of background bias vs. SNR. (d) Plot of background variance vs. SNR. (e) Plot of RMSE vs. SNR.
Figure 8 shows the cross sections reconstructed by pure sparse reconstruction and MLEM with sparse initial estimate for different truncation number. From Fig. 8, we notice that for small truncation number, pure sparse reconstruction generates blurry images. As truncation number increases, the resolution improves, but the background noise also increases due to the amplification of noise during preconditioning. For truncation number larger than 550, the signal is totally overwhelmed by the noise. As a comparison, the proposed method is able to largely reduce the background noise as truncation number increases. The RMSE as a function of truncation number is plotted in Fig. 9. The sparsity-initialized MLEM leads to lower RMSE for both noise levels.
Fig. 8 Cross sections of fluorescence target reconstructed with pure sparse reconstruction method for the top row and the proposed method for the bottom row with different truncation number K for digital mouse phantom for SNR=40dB.
Fig. 9 Plot of RMSE vs. truncation number for pure sparse reconstruction method and the proposed reconstruction method for different noise levels. (a) Plot of RMSE vs. truncation number for SNR=40 dB. (b) Plot of RMSE vs. truncation number for SNR=20 dB.
4. Discussion
In this paper, we have proposed an MLEM-based technique to reconstruct the fluorescence distribution from FMT data. In our framework, the initial estimate for the MLEM algorithm is derived from a sparse reconstruction method. Often an uniform initial estimate is used with MLEM-based techniques, but here we observe that a sparsity-initialized technique yields several advantages compared to uniformly initialized MLEM. First, sparsity-initialized MLEM has faster convergence speed. From Table 2, Fig. 3, Fig. 4(a) and Fig. 6(a), we observe that sparse initial estimate speeds up the convergence by both shortening the computation time for each iteration and requiring fewer iterations for convergence. In addition, sparsity-initialized MLEM also provides improved quantitative performance in ROI bias, background bias, ROI spatial variance, RMSE, and bias-variance trade-off compared to uniformly initialized MLEM, as shown in Figs. 4–7. Further, while results in both the cube phantom and the digital mouse phantom experiments indicate that the proposed method leads to higher ROI spatial variance compared to uniformly initialized MLEM for the same number of iterations, Fig. 4(c) and Fig. 6(c) show that the proposed method still provides better bias-variance trade-off compared to MLEM with uniform initial estimate. Further, sparsity-initialized MLEM often requires fewer iterations, which enables it to provide lower ROI spatial variance compared to uniformly initialized MLEM, as we observe from Fig. 5(b) and Fig. 7(b).
We also observe that sparsity-initialized MLEM provides advantages over pure sparse reconstruction method. From Fig. 8, we notice that sparsity-initialized MLEM is less sensitive to the choice of truncation number. For pure sparse reconstruction method, when the truncation number is small, the reconstructed image is blurry. As truncation number increases, image resolution is improved, but the noise in the background region is also increased due to the noise amplification during preconditioning. On the other hand, for small truncation number, sparsity-initialized MLEM is able to improve the resolution compared to pure sparse reconstruction method. For large truncation numbers, sparsity-initialized MLEM reduces noise in the background. These properties make MLEM with sparse initial estimate more robust to the choice of truncation number compared to pure sparse reconstruction method. The plots of RMSE vs. truncation number in Fig. 9 also demonstrate this point. Sparsity-initialized MLEM also improves quantitative performance of reconstructed images compared to pure sparse reconstruction method. Fig. 5 and Fig. 7 indicate this for different SNR values. Apart from improved background bias and spatial variance due to the reduction of background noise, sparsity-initialized MLEM also reduces the ROI bias compared to pure sparse reconstruction method, especially at low SNR values. At low SNR values, small truncation number is preferred to avoid noise amplification, which results in only a small number of measurements used for reconstruction. Small truncation number not only generates blurry images, as we discussed previously, but also causes severe bias in the reconstruction results. For example, for SNR= 5 dB, we observe that pure sparse reconstruction generated 71% ROI bias in the cube phantom experiments and 57% ROI bias in the digital mouse phantom experiments. As a comparison, MLEM uses the original system matrix and detector measurements for reconstruction, enabling it to compensate for the bias in the image, which reduced the ROI bias to 40% in the cube phantom experiments and 33% in the digital mouse phantom experiments.
We have observed, for example in Fig. 3, that sparsity-initialized MLEM is able to suppress the noise in the background region that is presented in the sparse initial estimate. To explain this observation, here we provide a theoretical justification. For a set of detector measurements denoted by Φ, consider two reconstructed images x1 and x2, where x1 is the image with noise in the background region (referred to as background noise), and x2 is the image that does not contain this background noise, as shown in Fig. 10(a) and (b), respectively. We denote the background noise as ϵ = x1 − x2, where ϵn ≥ 0 for all n. Before we proceed further, we introduce the concept of Kullback-Leibler (KL) distance. This distance measures how two probability distributions diverge from another. It is known that MLEM attempts to find an estimate that minimizes the KL distance between the measured data Φ and the data predicted by an estimate Gx. Thus, our objective is to assess whether the KL distance of x2 is less than x1, which would explain why MLEM would yield a solution x2 in comparison to x1. For x1, the KL distance is:
For x2, the KL distance is: Then the difference is: We denote . If ϕm ≤ (Gx2)m, fm ≤ 0 since fm is monotonically decreasing for (Gϵ)m ≤ 0 and fm(0) = 0. If ϕm > (Gx2)m, the plot of fm is shown in Fig. 10(c). To estimate the zeros of fm, we use second-order Taylor expansion to approximate fm, which yields: Let fm = 0, we get and where the inequality comes from the fact that ϕm > (Gx2)m. Also, note that the function fm has its maxima at (Gϵ)m = ϕm − (Gx2)m. Thus, if the detector response to the noise spot has a similar pattern as ϕm − (Gx2)m, fm will be close to its maximum for most detector indices m. This leads to a higher probability that ∆DKL > 0, which means MLEM tends to update towards noisy image. This is the case for the noise close to ROI. On the other hand, for noise spot in the background region, the detector response to noise spot will have a very different pattern compared to ϕm − (Gx2)m, as shown in Fig. 10(d). For detector index m where ϕm − (Gx2)m > 0, (Gϵ)m is either close to 0 or too large. This yields fm either close to zero or have a negative value. In this case, there is a higher probability that ∆DKL < 0, implying that MLEM tends to update towards results without the noise spot.Fig. 10 (a) Image with background noise. The noise spot in the background is marked with red circle. (b) Image without background noise. (c) Plot of fm. (d)Plot of (Gϵ)m and 2(ϕ − Gx2)m
The noise model in FMT is often assumed to be Gaussian [5,8–10,34,36]. In very few cases is the Poisson noise model applied [14]. Gaussian noise model is a good approximation when SNR is high, i.e. sufficient number of photons are detected. However, in some applications, the SNR value might be low, such as in brain imaging [41], dynamic FMT [42] and early-photon FMT [43]. Our results demonstrate that incorporating the Poisson noise model is especially valuable in these scenarios. More specifically, the pure sparse reconstruction method was formulated based on Gaussian noise model, while the proposed method incorporated both the sparsity information and Poisson noise model. We observe that the performance of the proposed method improves in comparison to the pure sparse reconstruction method as the SNR value decreases, and the proposed method is substantially more reliable at low SNR values. This shows the importance of accurately modeling Poisson noise for applications of FMT when insufficient number of photons are detected. We also point out that there might be some instances where the data is not a pure Poisson distribution. This could occur for example in case of CMOS detectors where the noise is a combination of Poisson and Gaussian distributions and is effected by electronic gain [44]. However, even in those cases, MLEM provides a convenient way to account for non-negativity constraints and enforces photon conservation. These factors may help improve the quality of reconstruction in comparison to mere use of sparsity-based methods. Thus, exploring the performance of this framework with systems that have non-Poisson noise distribution would be another interesting area of future study.
In this work, we only considered the case where the background uptake of fluorescence distribution is zero. While this is a common assumption in FMT studies [5, 8–10, 34–36], it is possible that the background uptake is not zero. Exploring the performance of the proposed method for this task would be an important future direction. The proposed method has been validated with extensive simulation experiments. Evaluating the performance of the method with physical phantom and in vivo animal experiments is another important direction of research. Finally, we used the MC-based method to model photon propagation to obtain the Green’s function in this work. However, there have been several analytical methods proposed for modeling light transport [45–50]. These methods can also be used to obtain an expression for the Green’s function. Analytical methods offer the advantage that they might be less sensitive to photon noise. Thus, implementing this reconstruction method using the analytical approaches is another area of further study.
5. Conclusion
We have presented a reconstruction framework for FMT involving sparsity-initialized MLEM. Simulation experiments on cube phantom and digital mouse phantoms demonstrate that the proposed method yields improved qualitative and quantitative performance compared to uniformly initialized MLEM as well as pure sparse reconstruction techniques. Further, compared to uniformly initialized MLEM, the proposed method is faster to execute, overcoming another barrier to application of MLEM technique for optical tomography. Moreover, compared to pure sparse reconstruction, the proposed method is more robust to noise amplification. We have also provided theoretical justification for the ability of the proposed method to reduce noise in the background region. Overall, this paper provides strong evidence that the proposed sparsity initialized MLEM-based reconstruction framework is feasible and advantageous for reconstruction in FMT imaging systems.
Funding
NIH BRAIN Initiative (Award R24 MH106083).
Acknowledgments
The authors thank Drs. Eric Frey and Jin Kang for helpful discussions.
Disclosures
Dean Wong acknowledges contract work with Lilly, Lundbeck, Intracellular, Five Eleven Pharma, Roche and Dart pharmaceuticals.
References and links
1. V. Ntziachristos, “Fluorescence molecular imaging,” Annu. Rev. Biomed. Eng. 8,1–33 (2006). [CrossRef] [PubMed]
2. N. C Biswal, A. Aguirre, Y. Xu, S. Zanganeh, Q. Zhu, C. Pavlik, M. B Smith, L. T Kuhn, and K. P Claffey, “Imaging tumor hypoxia by near-infrared fluorescence tomography,” J. Biomed. Opt. 16(6), 066009 (2011). [CrossRef] [PubMed]
3. F. Stuker, J. Ripoll, and M. Rudin, “Fluorescence molecular tomography: principles and potential for pharmaceutical research,” Pharmaceutics 3(2), 229–274 (2011). [CrossRef] [PubMed]
4. L. Zhou and B. Yazici, “Discretization error analysis and adaptive meshing algorithms for fluorescence diffuse optical tomography in the presence of measurement noise,” IEEE Trans. Image Process. 20(4), 1049–1111 (2011).
5. L. Zhao, H. Yang, W. Cong, G. Wang, and X. Intes, “Lp regularization for early gate fluorescence molecular tomography,” Opt. Lett. 39(14), 4156–4159 (2014). [CrossRef] [PubMed]
6. S. R. Arridge and J. C. Schotland, “Optical tomography: forward and inverse problems,” Inverse Probl. 25(12), 123010 (2009). [CrossRef]
7. V. Ntziachristos and R. Weissleder, “Experimental three-dimensional fluorescence reconstruction of diffuse media by use of a normalized Born approximation,” Opt. Lett. 26(12), 893–895 (2001). [CrossRef]
8. J. Ye, C. Chi, Z. Xue, P. Wu, Y. An, H. Xu, S. Zhang, and J. Tian, “Fast and robust reconstruction for fluorescence molecular tomography via a sparsity adaptive subspace pursuit method,” Biomed. Opt. Express 5(2), 387–406 (2014). [CrossRef] [PubMed]
9. D. Han, J. Tian, S. Zhu, J. Feng, C. Qin, B. Zhang, and X. Yang, “A fast reconstruction algorithm for fluorescence molecular tomography with sparsity regularization,” Opt. Express 18(8), 8630–8646 (2010). [CrossRef]
10. J. Dutta, S. Ahn, C. Li, S. R. Cherry, and R. M. Leahy, “Joint ℓ1 and total variation regularization for fluorescence molecular tomography,” Phys. Med. Biol. 57(6), 1459 (2012). [CrossRef] [PubMed]
11. O. Lee, J. Kim, Y. Bresler, and J. Ye, “Compressive diffuse optical tomography: noniterative exact reconstruction using joint sparsity,” IEEE Trans. Med. Imag. 30(5), 1129–1142 (2011). [CrossRef]
12. A. M. Bruckstein, D. L. Donoho, and M. Elad, “From sparse solutions of systems of equations to sparse modeling of signals and images,” SIAM Rev. 51(1), 34–81 (2009). [CrossRef]
13. V. Ntziachristos and R. Weissleder, “Charge-coupled-device based scanner for tomography of fluorescent near-infrared probes in turbid media,” Med. Phys. 29(5), 803–809 (2002). [CrossRef] [PubMed]
14. L. Adhikari, D. Zhu, C. Li, and R. F. Marcia, “Nonconvex reconstruction for low-dimensional fluorescence molecular tomographic poisson observations,” 2015 IEEE International Conference on Image Processing (ICIP), 2404–2408 (2015).
15. Q. Pian, R. Yao, L. Zhao, and X. Intes, “Hyperspectral time-resolved wide-field fluorescence molecular tomography based on structured light and single-pixel detection,” Opt. Lett. 40(3), 431–434 (2015). [CrossRef] [PubMed]
16. Y. Vardi, L. Shepp, and L. Kaufman, “A statistical model for positron emission tomography,” J. Am. Stat. Assoc. 80(389), 8–20 (1985). [CrossRef]
17. K. Lange, M. Bahn, and R. Little, “A theoretical study of some maximum likelihood algorithms for emission and transmission tomography,” IEEE Trans. Med. Imag. 6(2), 106–114 (1987). [CrossRef]
18. J. Llacer, E. Veklerov, K. J. Coakley, E. J. Hoffman, and J. Nunez, “Statistical analysis of maximum likelihood estimator images of human brain FDG PET studies,” IEEE Trans. Med. Imag. 12(2), 215–231 (1993). [CrossRef]
19. R. M. Lewitt and S. Matej, “Overview of methods for image reconstruction from projections in emission computed tomography,” Proc. IEEE 91(10), 1588–1611 (2003). [CrossRef]
20. A. K. Jha, E. Clarkson, M. A. Kupinski, and H. H. Barrett, “Joint reconstruction of activity and attenuation map using LM SPECT emission data,” Medical Imaging 2013: Physics of Medical Imaging 8668, 86681W (2013). [CrossRef]
21. G. Alexandrakis, F. R. Rannou, and A. F. Chatziioannou, “Tomographic bioluminescence imaging by use of a combined optical-PET (OPET) system: a computer simulation feasibility study,” Phys. Med. Biol. 50(17), 4225 (2005). [CrossRef] [PubMed]
22. N. V. Slavine, M. A. Lewis, E. Richer, and P. P. Antich, “Iterative reconstruction method for light emitting sources based on the diffusion equation,” Med. Phys. 33(1), 61–68 (2006). [CrossRef] [PubMed]
23. M. Jiang, T. Zhou, J. Cheng, W. Cong, and G. Wang, “Image reconstruction for bioluminescence tomography from partial measurement,” Opt. Express 15(18), 11095–11116 (2007). [CrossRef] [PubMed]
24. L. Cao and J. Peter, “Bayesian reconstruction strategy of fluorescence-mediated tomography using an integrated SPECT-CT-OT system,” Phys. Med. Biol. 55(9), 2693 (2010). [CrossRef] [PubMed]
25. J. Qi and R. M. Leahy, “Iterative reconstruction techniques in emission computed tomography,” Phys. Med. Biol. 51(15), R541 (2006). [CrossRef] [PubMed]
26. S. Ahn, A. J. Chaudhari, F. Darvas, C. A. Bouman, and R. M. Leahy, “Fast iterative image reconstruction methods for fully 3D multispectral bioluminescence tomography,” Phys. Med. Biol. 53(14), 3921 (2008). [CrossRef] [PubMed]
27. H. H. Barrett and K. J. Myers, Foundations of Image Science, 1 ed. (Wiley, 2004).
28. D. Ma, P. Wolf, A. V. Clough, and T. G. Schmidt, “The performance of MLEM for dynamic imaging from simulated few-view, multi-pinhole SPECT,” IEEE Trans. Nucl. Sci. 60(1), 115–123 (2013). [CrossRef]
29. Y. Zhu, A. K. Jha, J. K. Dreyer, H. N. D. Le, J. U. Kang, P. E. Roland, D. F. Wong, and A. Rahmim, “A three-step reconstruction method for fluorescence molecular tomography based on compressive sensing,” Proc. SPIE 10059, 1005911 (2017). [CrossRef]
30. Y. Zhu, A. K. Jha, D. Wong, and A. Rahmim, “Improved sparse reconstruction for fluorescence molecular tomography with poisson noise modeling,” Biophotonics Congress: Biomedical Optics Congress 2018, JTu3A.51 (2018).
31. G. Kontaxakis and L. G Strauss, “Maximum likelihood algorithms for image reconstruction in Positron Emission Tomography,” Radionuclides Oncol. 8, 73–106 (1998)
32. J. H. Chang, J. MM Anderson, and J. R Votaw, “Regularized image reconstruction algorithms for Positron Emission Tomography,” IEEE Trans. Med. Imag. 23(9), 1165–1175 (2004). [CrossRef]
33. A. Jin, B. Yazici, A. Ale, and V. Ntziachristos, “Preconditioning of the fluorescence diffuse optical tomography sensing matrix based on compressive sensing,” Opt. Lett. 37(20), 4326–4328 (2012). [CrossRef] [PubMed]
34. J. Shi, X. Cao, F. Liu, B. Zhang, J. Luo, and J. Bai, “Greedy reconstruction algorithm for fluorescence molecular tomography by means of truncated singular value decomposition conversion,” J. Opt. Soc. Am. A 30(3), 437–447 (2013). [CrossRef]
35. A. Jin, B. Yazici, and V. Ntziachristos, “Light illumination and detection patterns for fluorescence diffuse optical tomography based on compressive sensing,” IEEE Trans. Med. Imag. 23(6), 2609–2624 (2014). [CrossRef]
36. R. Yao, Q. Pian, and X. Intes, “Wide-field fluorescence molecular tomography with compressive sensing based preconditioning,” Biomed. Opt. Express 6(12), 4887–4898 (2015). [CrossRef] [PubMed]
37. A. Beck and M. Teboulle, “A fast iterative shrinkage-thresholding algorithm for linear inverse problems,” SIAM J. Imaging Sci. 2(1), 183–202 (2009). [CrossRef]
38. Q. Fang and D. A. Boas, “Monte carlo simulation of photon migration in 3D turbid media accelerated by graphics processing units,” Opt. Express 17(22), 20178–20190 (2009). [CrossRef] [PubMed]
39. W. P. Segars, B. M. W. Tsui, E. C. Frey, G. A. Johnson, and S. S. Berr, “Development of a 4-D digital mouse phantom for molecular imaging research,” Mol. Imag. Biol. 6(3), 149–159 (2004). [CrossRef]
40. G. Strangman, M. A. Franceschini, and D. A. Boas, “Factors affecting the accuracy of near-infrared spectroscopy concentration calculations for focal changes in oxygenation parameters,” Neuroimage 18(4), 865–879 (2003). [CrossRef] [PubMed]
41. SB Raymond, ATN Kumar, DA Boas, and BJ Bacskai, “Optimal parameters for near infrared fluorescence imaging of amyloid plaques in alzheimer’s disease mouse models,” Phys. Med. Biol. 54(20), 6201 (2009). [CrossRef] [PubMed]
42. X. Liu, F. Liu, Y. Zhang, and J. Bai, “Unmixing dynamic fluorescence diffuse optical tomography images with independent component analysis,” IEEE Trans. Med. Imag. 30(9), 1591–1604 (2011). [CrossRef]
43. F. Leblond, H. Dehghani, D. Kepshire, and B. W Pogue, “Early-photon fluorescence tomography: spatial resolution improvements and noise stability considerations,” J. Opt. Soc. Am. A 26(6), 1444–1457 (2009). [CrossRef]
44. A. K. Jha, “Retrieving information from scattered photons in medical imaging,” Ph.D dissertation, The University of Arizona, 2013.
45. O. Lehtikangas, T. Tarvainen, and A. D. Kim, “Modeling boundary measurements of scattered light using the corrected diffusion approximation,” Biomed. Opt. Express 3(3), 552–571 (2012). [CrossRef] [PubMed]
46. T. Tarvainen, M. Vauhkonen, V. Kolehmainen, S. R Arridge, and J. P Kaipio, “Coupled radiative transfer equation and diffusion approximation model for photon migration in turbid medium with low-scattering and non-scattering regions,” Phys. Med. Biol. 50(20), 4913 (2005). [CrossRef] [PubMed]
47. PS. Mohan, T. Tarvainen, M. Schweiger, A. Pulkkinen, and S. R Arridge, “Variable order spherical harmonic expansion scheme for the radiative transport equation using finite elements,” J. Comput. Phys. 230(19), 7364–7383 (2011). [CrossRef]
48. A. K. Jha, M. A. Kupinski, T. Masumura, E. Clarkson, A. V. Maslov, and H. H. Barrett, “Simulating photon-transport in uniform media using the radiative transport equation: a study using the Neumann-series approach,” J. Opt. Soc. Am. A 29(8), 1741–1757 (2012). [CrossRef]
49. A. K Jha, M. A Kupinski, H. H Barrett, E. Clarkson, and J. H Hartman, “Three-dimensional Neumann-series approach to model light transport in nonuniform media,” J. Opt. Soc. Am. A 29(9), 1885–1899 (2012). [CrossRef]
50. A. K. Jha, Y. Zhu, S. Arridge, D. F. Wong, and A. Rahmim, “Incorporating reflection boundary conditions in the Neumann series radiative transport equation: application to photon propagation and reconstruction in diffuse optical imaging,” Biomed. Opt. Express 9(4), 1389–1407 (2018). [CrossRef] [PubMed]
