This article has Open Peer Review reports available.

# Homotopic non-local regularized reconstruction from sparse positron emission tomography measurements

- Alexander Wong
^{1}Email author, - Chenyi Liu
^{2}, - Xiao Yu Wang
^{1}, - Paul Fieguth
^{1}and - Hongxia Bie
^{2}

**15**:10

https://doi.org/10.1186/s12880-015-0052-5

© Wong et al.; licensee BioMed Central. 2015

**Received: **28 May 2014

**Accepted: **25 February 2015

**Published: **18 March 2015

## Abstract

### Background

Positron emission tomography scanners collect measurements of a patient’s *in vivo* radiotracer distribution. The system detects pairs of gamma rays emitted indirectly by a positron-emitting radionuclide (tracer), which is introduced into the body on a biologically active molecule, and the tomograms must be reconstructed from projections. The reconstruction of tomograms from the acquired PET data is an inverse problem that requires regularization. The use of tightly packed discrete detector rings, although improves signal-to-noise ratio, are often associated with high costs of positron emission tomography systems. Thus a sparse reconstruction, which would be capable of overcoming the noise effect while allowing for a reduced number of detectors, would have a great deal to offer.

### Methods

In this study, we introduce and investigate the potential of a homotopic non-local regularization reconstruction framework for effectively reconstructing positron emission tomograms from such sparse measurements.

### Results

Results obtained using the proposed approach are compared with traditional filtered back-projection as well as expectation maximization reconstruction with total variation regularization.

### Conclusions

A new reconstruction method was developed for the purpose of improving the quality of positron emission tomography reconstruction from sparse measurements. We illustrate that promising reconstruction performance can be achieved for the proposed approach even at low sampling fractions, which allows for the use of significantly fewer detectors and have the potential to reduce scanner costs.

## Background

The concept of positron emission tomography (PET) imaging revolves around the measurement of a patient’s *in vivo* radiotracer distribution. The system detects pairs of gamma rays emitted indirectly by a positron-emitting radionuclide (tracer), which is introduced into the body via biologically active molecules. As such, PET imaging can be described via a line-integral model of acquisition [1]. PET data are collected as projections in sinograms or listmode format [2]. The raw data collected by a PET scanner are a list of ‘coincidence events’ representing near-simultaneous detection, each coincidence event represents a line in space connecting the two detectors along which the positron emission occurred. While one can employ the inverse Radon transform to recover tomograms from acquired PET data, such an approach can be unstable, particularly when dealing with noise-contaminated data [3]. Therefore, the filtered back projection algorithm, which is a stabilized and discretized version of the inverse Radon transform, is commonly used in practice [3].

There are two common approaches for tomographic reconstruction: i) filtered back projection (FBP) [4], and ii) iterative expectation-maximization (EM) [5-7]. FBP has been widely used to reconstruct tomograms from the projections in clinical settings due to its overall simplicity and low computational complexity. However, shot noise in the raw data is prominent in the reconstructed tomograms and areas of high tracer uptake tend to form streaks across the tomogram.

Iterative algorithms for reconstructing tomograms from acquired PET data take into account the statistical nature of the acquired projection data and incorporate the physical model into reconstruction. These algorithms compute an estimate of the likely distribution of annihilation events that led to the measured data, based on statistical principles [8]. The advantage is a better noise profile and resistance to the artifacts common with FBP, but are typically more computationally intensive.

Recently, wavelet-based methods [3,8-10] have also been proposed for solving this inverse problem in PET reconstruction. Wavelet-based methods have been applied in research literature in the post-processing stage to suppress artifacts in the reconstructed tomograms [3]. Another attractive approach to reconstructing tomograms from PET data using wavelets is to incorporate wavelets into the reconstruction process itself [9], where an *L*_{1}-regularization constraint is used to enforce sparsity in the wavelet domain. For example, wavelet regularization using exponential-spline wavelets has been shown to produce good reconstruction results [10].

More recently, the concept of compressive sensing [11] has offered great potential for signal recovery with very high accuracy through sparse measurements acquired at very low sampling fractions. Researchers have also begun to explore the use of compressive sensing for image reconstruction in various medical modalities [12]. To date, most sparse medical image reconstruction algorithms have focused on solving the sparse reconstruction problem via *L*_{1} minimization.

The main contribution of this paper is to introduce and investigate the potential of a homotopic non-local regularization (HNLR) reconstruction framework for the reconstruction of PET tomograms from sparse measurements, with the aim to achieve high reconstruction quality. The ability to achieve strong reconstruction results at low sampling fractions can be very helpful for reducing the number of detectors needed in a PET system [7], which in effect helps can lead to reduced scanner costs and potentially improve adoption of PET in developing countries.

The rest of the paper is organized as follows. First, the underlying methodology behind the proposed use of a homotopic non-local regularization (HNLR) reconstruction framework for the reconstruction of PET tomograms from sparse measurements is described in Section “Methods”. The experimental results are presented in Section “Results and discussion”. Finally, conclusions are drawn and future work is discussed in Section “Conclusion”.

## Methods

where *P* is the set of sampled measurements, *s* is the distance of straight line *L* from the origin, *α* is the angle that the normal vector to *L* makes with the *x* axis,
is the Radon transform, and *f* is the unknown tomogram. The goal of reconstruction is to use the sampled measurements *P* (projections in the Radon transform domain through the object) to find the tomogram *f* in image space.

*L*_{
p
} minimization for sparse reconstruction

*Φ*is a measurement operator defining which of the sites are measured (non-measured projection space indices will be assigned 0 by

*Φ*). Our goal is to reconstruct

*f*(

*x*) from a sparse sampling of

*P*(

*α*,

*s*). With only few measurements made, not all sites are measured, making this problem an ill-posed inverse problem [13], and as such there may exist multiple solutions. A common method to solving this inverse problem is through the use of

*L*

_{2}minimization:

where \(\hat {f}(x)\) is the estimated signal in image space, and \(\hat {P}(\alpha,s)\) is the estimated signal in projection space. However, solving the problem in this manner results in significant artifacts such as aliasing and blurring in practical situations [12].

*f*(

*x*) can be obtained by maximizing the sparsity of the signal in the transform domain and enforcing data fidelity in projection space domain. One can formulate the reconstruction problem as a constrained

*L*

_{0}minimization problem,

where *Ψ* is the sparse transform operator. Two commonly used sparse transform operators are the finite differential transformation and the wavelet transformation [16].

*L*

_{0}problem is essentially NP hard, and as such is intractable in practice [17]. To address this problem, pioneering work done by Candes [11] and Donoho [15] showed that under certain conditions,

*L*

_{0}and

*L*

_{1}minimization can lead to the same solutions. As such, one can instead solve the

*L*

_{1}minimization as:

Theoretically, under certain conditions, solving the *L*_{1} problem can get exactly the same solution as solving the *L*_{0} problem, although a substantial increase in the number of measurements is required [12,18].

*L*

_{1}minimization framework, which is known to have an edge-preservation effect, and to account for the unavoidable noise in the measurements, the data fidelity constraint is typically replaced by a

*L*

_{2}norm constraint:

### Homotopic, non-local regularization based reconstruction

There are two fundamental limitations associated with the reconstruction framework described in Eq. 6 for the purpose of reconstructing tomograms [20]. First, the use of a TV penalty can lead to the loss of fine structure at low sampling fractions [12]. Since PET images are characterized by complex functional process variations, the TV approach may not be well-suited for reconstructing these data. Second, as already mentioned, the number of measurements required for solving the *L*_{1} problem can be noticeably higher than that for the *L*_{0} problem, which is undesirable as it leads to an increased number of detectors needed. Hence, an alternative reconstruction strategy that addresses these two problems is desired.

*L*

_{0}minimization approach without a drastic increase in the number of measurements, Trzasko et al. [12] introduced a homotopic

*L*

_{0}minimization framework, which uses an increasing approximation framework to get close to the

*L*

_{0}minimization problem:

where *ρ* is the homotopic approximation of the *L*_{0} norm, which approaches the *L*_{0} norm as *σ* approaches 0. Experimental results showed that this strategy is able to approach the capabilities of the *L*_{0} minimization approach [21], thus addressing the problem pertaining to the number of measurements required.

*L*

_{0}minimization framework [25] for the purpose of tomogram reconstruction from sparse PET measurements, as it is well-suited for handling fine structural details. The proposed homotopic, non-local regularization framework for sparse PET reconstruction is formulated as follows:

*η*denotes the homotopic non-local regularization functional:

*σ*controls the approximation degree, and

*w*(

*x*,

*y*,

*σ*) is defined as

*N*(

*x*) is a neighborhood around

*x*. Essentially, the idea is to minimize the homotopic non-local regularization function, denoted by

*η*(|

*Φ*

*f*(

*x*)|,

*σ*), at decreasing values of approximation degree, denoted by

*σ*, such that

*η*approximates the

*L*

_{0}norm closer and closer. An approximate solution of Eq. 8 can be obtained using an iterative optimization approach, with decreasing values of

*σ*as the number of iterations increases. To reduce the computation complexity, the neighborhood search space is limited to a window search around the pixel to be estimated. The regularization term

*η*(|

*Ψ*

*f*(

*x*)|,

*σ*) is designed to efficiently suppress artifacts associated with incomplete measurements while preserving detail and structure, and is enforced via steepest descent [26,27], while the data fidelity term is designed to ensure that the reconstructed tomogram complies with sampled measurements while accounting for some level of noise artifacts. The data fidelity term is enforced at measured projection space indices at each iteration

*i*as follows:

where *H*(.) is the Heaviside step function and *b* is the number of measured projection space indices. The pseudocode for the proposed homotopic non-local regularized reconstruction method is shown in Algorithm 1. Note that a *L*_{2} data term is used here for simplicity in the iterative optimization realization shown in Algorithm 1, and more advanced data terms such as the Kullback-Leibler divergence [28,29], which takes better advantage of the Poisson noise statistics of PET, may be used.

### Experiments

#### Data description

To evaluate the effectiveness of the proposed method, sparse reconstructions were performed using a simulated phantom generated using ASIM [30] at 30% sampling fraction (defined here as the percentage of projection angles from which measurements are obtained). The simulated phantom consists of 4 spots of different sizes and is useful for observing the effects of reduced sampling on spatial resolution. Also, since the simulated phantom generated using ASIM is contaminated by simulated Poisson-distributed noise, it is also useful for observing the performance of the proposed method when faced with noisy data. Since the total number of projection angles used in this study is 180 projection angles, that means at a 30% sampling fraction the number of projection angles used for sparse reconstruction is 54 projection angles. Also, the object covers approximately 70% of the field of view. Furthermore, sparse reconstructions at different sampling fractions were performed using three real PET data sets hosted by Harvard Medical School [31], in accordance with Harvard Medical School ethics procedures. The PET data sets are acquired with Fluorine-18 Deoxyglucose (FDG) as the biologically active molecule. The first PET data set (PET1) is acquired from a 73 year-old man who sought medical attention due to a grand mal seizure and progressive speech difficulties, with a brain biopsy revealing grade II astrocytoma. The second PET data set (PET2) was taken from a 53 year-old man who sought medical attention due to a grand mal seizure, with a brain biopsy revealing grade IV astrocytoma. The third PET data set (PET3) is acquired from a 70 year-old man with mild Alzheimer’s disease about 9 months prior to imaging. The PET data sets are 256×256×23 voxels. The data sets were forward-projected from a previous reconstruction obtained via filtered back-projection (FBP) [2] based on measurements from 100% of the projection angles in order to obtain projection data for evaluation.

#### Implementation

For comparison purposes, the proposed homotopic non-local regularization (HNLR) algorithm was compared with the standard filtered back-projection (FBP) algorithm [2], as well as the more advanced iterative expectation maximization algorithm with total variation regularization (TVR) algorithm [6]. The standard FBP algorithm was chosen for baseline comparisons as it is widely used in clinical settings for its speed and ease of implementation. The TVR algorithm was chosen as it is one of the more advanced reconstruction approaches available for PET reconstruction [6,7]. To allow for quantitative assessment, the signal-to-noise ratio (SNR) was computed for reconstructed data at different sampling fractions. Furthermore, a qualitative visual assessment is performed on the reconstructed data to investigate the reconstruction performance and the preservation of details achieved using the tested methods. Different sampling fractions were achieved by sampling the projection angles at evenly spaced intervals. For example, to achieve a 50*%* sampling fraction, every other projection angle from the set of total projection angles was sampled. This methodology was used for all the data sets and the comparison of all three methods. For this study, the parameters of the HNLR algorithm were set as follows as it yielded strong reconstruction performance. The neighborhood *N*(*x*) was chosen as a 9×9 neighborhood around *x*, *σ*_{1} and *σ*_{
t
} are set as 20% and 1% of the dynamic range in image space, respectively, *λ* is set as 0.8, *ε* is set as 2% of the dynamic range in projection space multiplied by the number of measured projection space indices (*b*), and *n* is set as 100 as these parameters were found empirically to provide strong reconstruction performance. Note that the parameters of the HNLR algorithm can be tuned based on the PET imaging system the algorithm is incorporated into to obtain optimal performance for clinical data obtained from the PET system.

### Results and discussion

where *f*(*x*) is original image, \(\hat {f}(x)\) is reconstructed image, and *N* is the number of pixels in each image.

## Conclusion

The potential use of homotopic non-local regularization (HNLR) reconstruction framework to maintaining high reconstruction quality while reducing sampling fraction in sparse PET reconstruction was studied. The proposed reconstruction framework is designed to address issues associated with artifacts and detail loss with sparse PET reconstruction. A comparative analysis using real PET data sets was performed to compare traditional FBP as well as TVR reconstruction with the proposed HNLR reconstruction method. Based on experimental results, we illustrate that promising reconstruction performance can be achieved using HNLR even at low sampling fractions. Future work includes a more comprehensive evaluative study using larger PET data sets to further validate the performance of the proposed HNLR framework, the investigation of different sparsifying transforms to study their potential for improving reconstruction quality, and the investigation of incorporating alternative representations such as multiscale texture representations and local phase [32,33] for enforcing regularization. Furthermore, we wish to investigate more advanced data terms such as the Kullback-Leibler divergence [28,29], which takes better advantage of the Poisson noise statistics of PET, as well as investigate the application of the HNLR framework for other medical imaging modalities such as diffusion weighted magnetic resonance imaging [34,35] and correlated diffusion imaging [36]. We also aim to investigate methods for automatically optimizing the parameters of the HNLR framework to obtain optimal reconstruction accuracy. Finally, we wish to investigate the effects of employing the HNLR framework for sparse reconstruction on medical image analysis tasks such as medical image annotation and contouring [37], as well as medical image registration [32,33,38-41].

## Declarations

### Acknowledgements

This work was supported by the Natural Sciences and Engineering Research Council (NSERC) of Canada, the Canada Research Chairs program, the Ontario Ministry of Economic Development and Innovation, and the Chinese Council Scholarship. The authors would also like to thank Dr. Keith A. Johnson for the test data.

## Authors’ Affiliations

## References

- Valk P, Bailey D, Townsend D, Maisey M. Positron emission tomography. London: Springer; 2002.Google Scholar
- Brinks R, Busch M. Local compensation for respiratory motion in list-mode PET. Springer Proceedings in Physics: Adv Med Imaging. 2007; 114:31–6.View ArticleGoogle Scholar
- Kalifa J, Laine A, Esser P. Tomographic reconstruction with non-linear diagonal estimators. Proc SPIE. 2000; 4119:1–11.View ArticleGoogle Scholar
- Wang W, Gindi G. Noise analysis of MAP-EM algorithms for emission tomography. Phys Med Biol. 1997; 42:2215–32.View ArticlePubMedGoogle Scholar
- Wang C, Snyder W, Bilbro G, Santago P. Performance evaluation of filtered backprojection reconstruction and iterative reconstruction methods for PET images. Comput Biol Med. 1998; 28:13–24.View ArticlePubMedGoogle Scholar
- Jonsson E, Huang S, Chan T. Total variation regularization in positron emission tomography. Technical report, Los, Angeles, California: Dept Mathematics, Univ California; 1998, pp. 1–25.Google Scholar
- Valiollahzadeh S, Chang T, John W Clark J, Mawlawi OR. Image recovery in pet scanners with partial detector rings using compressive sensing. IEEE Nucl Sci Symp Med Imaging Confererence. 2012:3036–9.Google Scholar
- Zhou J, Coatrieux J, Bousse A, Shu H, Luo L. A bayesian map-em algorithm for PET image reconstruction using wavelet transform. IEEE Trans on Nucl Sci. 2007; 54:1660–9.View ArticleGoogle Scholar
- Pustelnik N, Chaux C, Pesquet JC, Sureau FC, Dusch E. Adapted convex optimization algorithm for wavelet-based dynamic PET reconstruction. In: the proceeding of the 10th, International Meeting on Fully Three-Dimensional Image Reconstruction in Radiology and Nuclear Medicine: 2009. p. 5–10.Google Scholar
- Verhaeghe J, Ville DVD, Khalidov I, Asseler YD, Lemahieu I, Unser M. Dynamic PET reconstruction using wavelet regularization with adapted basis functions. IEEE Trans Med Image. 2008; 27:943–59.View ArticleGoogle Scholar
- Candës E, Romberg J, Tao T. Robust uncertainty principles: exact signal reconstruction from highly incomplete frequency information. IEEE Trans Inf Theory. 2006; 52:489–509.View ArticleGoogle Scholar
- Trzasko J, Manduca A. Highly undersampled magnetic resonance image reconstruction via homotopic l0-minimization. IEEE Trans Med Imaging. 2009; 28:106–21.View ArticlePubMedGoogle Scholar
- Fieguth P. Statistical image processing and multidimensional modeling. New York: Springer; 2010.Google Scholar
- Hanif A, Mansoor A, Ejaz T. Iterative tomographic image reconstruction by compressive sampling. In: Image Process (ICIP): 2010. p. 4313–6.Google Scholar
- Donoho D. Compressive sensing. IEEE Trans Inf Theory. 2006; 52:1289–306.View ArticleGoogle Scholar
- Lustig M, Donoho D, Pauly J. Sparse MRI: The application of compressed sensing for rapid MR imaging. Magn Reson Med. 2007; 58:1182–95.View ArticlePubMedGoogle Scholar
- Natarajan BK. Sparse approximate solutions to linear systems. SIAM J Comput. 1995; 24:227–34.View ArticleGoogle Scholar
- Wong A, Mishra A, Clausi D, Fieguth P. Sparse reconstruction of breast MRI using homotopic l0 minimization in a regional sparsified domain. Biomed Eng IEEE Trans. 2010; 60:743–52.View ArticleGoogle Scholar
- Gilboa G, Osher S. Nonlocal operators with applications to image processing. Multiscale Model Simul. 2008; 7:1005–28.View ArticleGoogle Scholar
- Liang D, Wang H, Ying L. Sense reconstruction with nonlocal TV regularization. Proc IEEE Eng Med Biol Soc. 2009; 2009:1032–5.Google Scholar
- Trzasko J, Manduca A, Borisch E. Sparse MRI reconstruction via multiscale L0-Continuation. In: Proc IEEE Workshop Stat Signal Process: 2007. p. 176–80.Google Scholar
- Wong A, Mishra A, Fieguth P, Clausi D. A perceptually adaptive approach to image denoising using anisotropic non-local means. In: Proc Int Conference Image Process: 2008. p. 1–4.Google Scholar
- Wong A, Orchard J. An adaptive non-local means approach to exemplar-based inpainting. In: Proc Int, Conference Image Process: 2008. p. 2600–3.Google Scholar
- Orchard J, Ebrahimi M, Wong A. Efficient nonlocal-means denoising using the svd. In: Proc Int Conference Image Process: 2008. p. 1732–5.Google Scholar
- Liu C, Wong A, Fieguth P, Bizheva K, Bie H. Homotopic, non-local sparse reconstruction of optical coherence tomography (OCT) imagery. Opt Express. 1020; 20:0–11.Google Scholar
- Rosenbloom P. The method of steepest descent. Proc Symp Appl Math. 1956; 6:127–76.View ArticleGoogle Scholar
- Wang X. Method of steepest descent and its applications. IEEE Microwave Wireless Components Lett. 2008; 12:24–26.Google Scholar
- Pustelnik N, Chaux C, Pesquet JC, Comtat C. Parallel algorithm and hybrid regularization for dynamic pet reconstruction. In: IEEE Nucl Sci Symp Med Imaging Confererence: 2010. p. 2423–7.Google Scholar
- Anthoine S, Aujol J, Mélot C, Boursier Y. Some proximal methods for cbct and pet tomography. inverse problems in imaging. Inverse Probl Imaging. 2012; 6:565–98.View ArticleGoogle Scholar
- Anthoine S, Aujol J, Mélot C, Boursier Y. Asim pet simulation software. 2012. http://depts.washington.edu/asimuw/.
- Johnson K, Becker J. The whole brain atlas. 2008. www.med.harvard.edu/AANLIB/home.html.
- Wong A, Orchard J. Robust multimodal registration using local phase-coherence representations. J Signal Process Syst. 2009; 54:89–100.View ArticleGoogle Scholar
- Wong A, Clausi D, Fieguth P. Cpol: Complex phase order likelihood as a similarity measure for mr-ct registration. Med Image Anal. 2010; 14:50–7.View ArticlePubMedGoogle Scholar
- Koh D, Padhani A. Diffusion-weighted mri: a new functional clinical technique for tumour imaging. Br J Radiology. 2006; 79:633–5.View ArticleGoogle Scholar
- Shafiee M, Lui D, Haider S, Cameron A, Wong A, et al. Apparent ultra-high b-value diffusion-weighted image reconstruction via hidden conditional random fields. IEEE Trans Med Imaging. 2015; 1:1–15.View ArticleGoogle Scholar
- Wong A, Glaister J, Cameron A, Haider M. Correlated diffusion imaging. BMC Med Imaging. 2013; 13:1–7.View ArticleGoogle Scholar
- Mishra A, Wong A, Zhang W, Fieguth P, Clausi D. Improved interactive medical image segmentation using enhanced intelligent scissors (eis). In: Annu Int Conference IEEE Eng Med Biol Soc: 2008. p. 1–4.Google Scholar
- Wong A, Bishop W, Orchard J. Efficient multi-modal least-squares alignment of medical images using quasi-orientation maps. In: Proc Int Conference Image Process, Comput Vision, Pattern Recognit: 2006. p. 74–80.Google Scholar
- Wong A, Bishop W. Efficient least squares fusion of mri and ct images using a phase congruency model. Pattern Recognit Lett. 2008; 29:173–80.View ArticleGoogle Scholar
- Roche A, Malandain G, Pennec X, Ayache N. The correlation ratio as a new similarity measure for multimodal image registration. Lecture Notes Comput Sci: Med Image Comput Comput-Assisted Intervention. 1998; 1496:1115–24.View ArticleGoogle Scholar
- Wong A, Fieguth P. Fast phase-based registration of multimodal image data. Signal Process. 2009; 89:724–37.View ArticleGoogle Scholar

## Copyright

This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly credited. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.