Poisson noise obscures hypometabolic lesions in PET.

The technology of fluoro-deoxyglucose positron emission tomography (PET) has drastically increased our ability to visualize the metabolic process of numerous neurological diseases. The relationship between the methodological noise sources inherent to PET technology and the resulting noise in the reconstructed image is complex. In this study, we use Monte Carlo simulations to examine the effect of Poisson noise in the PET signal on the noise in reconstructed space for two pervasive reconstruction algorithms: the historical filtered back-projection (FBP) and the more modern expectation maximization (EM). We confirm previous observations that the image reconstructed with the FBP biases all intensity values toward the mean, likely due to spatial spreading of high intensity voxels. However, we demonstrate that in both algorithms the variance from high intensity voxels spreads to low intensity voxels and obliterates their signal to noise ratio. This finding has profound impacts on the clinical interpretation of hypometabolic lesions. Our results suggest that PET is relatively insensitive when it comes to detecting and quantifying changes in hypometabolic tissue. Further, the images reconstructed with EM visually match the original images more closely, but more detailed analysis reveals as much as a 40 percent decrease in the signal to noise ratio for high intensity voxels relative to the FBP. This suggests that even though the apparent spatial resolution of EM outperforms FBP, the signal to noise ratio of the intensity of each voxel may be higher in the FBP. Therefore, EM may be most appropriate for manual visualization of pathology, but FBP should be used when analyzing quantitative markers of the PET signal. This suggestion that different reconstruction algorithms should be used for quantification versus visualization represents a major paradigm shift in the analysis and interpretation of PET images.

Introduction

Positron emission tomographic (PET) images play a major role in the treatment and management of a growing number of maladies. In most cases, the interpretation of these images relies on the detection of high intensity lesions by quantifying the relative distribution of a radioactively decaying tracer. This tracer is most commonly fluoro-deoxyglucose (FDG), which allows PET to quantify the relative glucose metabolism in tissues. Hypermetabolic lesions indicate the presence of inflammation, malignancy, and/or major functional changes. The observation of these changes has been critical to the characterization and clinical management of central nervous system cancers, paraneoplastic syndrome, Huntington’s, and, when scanned during ictus, epilepsy [1,2]. In some cases, PET is used to guide resective curative neurosurgery [1,3].

The detection of hypometabolic lesions is equally clinically relevant. PET has been effectively used to characterize Parkinson’s disease, Alzheimer’s disease, interictal epilepsy, cortical dysplasia, tuberosclerosis, and even mood disorders. In these cases, the hypometabolic lesions indicate functional abnormalities or portend the location of future atrophic lesions [1-12]. In cortical dysplasia, tuberous sclerosis and central nervous system infections with multiple structural abnormalities frequently exist, but only a small subset of these lesions generates epileptic seizures. The co-localization these structural abnormalities visualized in MRI with hypometabolic lesions observed using PET can be effectively used to determine which of these structural abnormalities is generating the seizures [1]. When these co-localized lesions are resected, 86 percent of patients achieve favorable outcomes compared to 30 percent to 76 percent without co-localization [13-15].

PET, however, may be biased against the detection of these hypometabolic lesions. The technology of PET relies on the emission of positrons from radioactively decaying isotopes. The number of positrons that are emitted from each volumetric pixel, or voxel, is Poisson distributed. In a Poisson distribution, the variance of a sample is equal to its mean. Consequently, the variance of positron count increases as overall signal increases. Due to the fact that these images are reconstructed based on projections, this noise could potentially spread to nearby voxels [16,17]. For hypermetabolic lesions, this would result in the lesion dominating the signal in the reconstructed images. Unfortunately, the same signal interaction can allow surrounding normal tissue to mask hypometabolic tissue.

This potential bias against the detection of hypometabolic lesions may seem to be an issue of resolution. Modern reconstruction techniques such as ordered subset expectation maximization (OSEM) substantially increase the resolution of reconstructed images relative to the canonical filtered back-projection (FBP) [18,19]. One of the major hurdles to resolution in FBP is the streaking caused by high intensity voxels. In X-ray computed tomographic (CT) imaging, this streaking is regularly caused by bone artifacts. It is also present, albeit to a lesser degree, in PET [20,21]. The OSEM algorithm substantially decreases the effect of these streaks and thereby increases image resolution [22,23]. This improvement is visually apparent even to the untrained observer and has resulted in the pervasive adoption of the OSEM algorithm for CT and PET reconstruction.

These improvements, however, only focus on decreasing the bias in signal intensity caused by surrounding tissues and ignore effects of noise. A simple theoretical proof illustrates that maximum variance of voxel intensity in image space provides an upper bound for the maximum reconstructed voxel intensity variance for FBP (see Supplementary Material) [19]. There is no analogous proof for the iterative EM algorithm, much less the OSEM algorithm. Therefore, even though the OSEM algorithm decreases bias, it has the potential to increase variance and thereby decrease the signal to noise ratio. This potential challenge has been largely ignored because, as humans, we are exquisitely capable of detecting changes in the mean but relatively weak at detecting changes in the spread. The development of PET and CT reconstruction has focused on the generation of visually interpretable images; therefore, previous literature has focused almost exclusively on trends of the mean. As more quantifiable markers of PET are developed, we believe that an in-depth treatment of the variance is critically important to achieving accurate and clinically relevant measurements.

In this paper, we use Monte Carlo simulations to characterize the statistical properties of the variance in both EM and FBP reconstructed images. We demonstrate that in both algorithms, the Poisson noise from hypermetabolic voxels obliterates the signal to noise ratio for hypometabolic lesions, resulting in a bias against the detection of hypometabolic lesions. The understanding of this effect has a profound impact on the interpretation of hypometabolic lesions on PET images.

Materials and Methods

In this Monte Carlo simulation, 10 million exemplars of reconstructed Poisson noise were measured from pixels with integer initial intensity from 1 to 100. Figure 1 illustrates one example of a noisy image and the two reconstructions with a common intensity scale. Images were sequentially realized until each intensity value had been reconstructed at least 100,000 times. This took 2,614 realizations and approximately 7 cpu-days. We constructed our simulation to match PET images of an average human brain with (2mm)3 voxels. Each simulation image was 80 voxels by 80 voxels with circle of radius 35 voxels centered on the 40th voxel in each dimension. Each voxel within the circle was given an initial intensity from a discrete uniform distribution ranging from 1 to 100. All voxels outside the circle had intensity 0. Poisson noise with parameter equal to the initial intensity of each voxel was then added. Radon projections were used to simulate the actual data collected by sensors for integer angles from 0 to 179 degrees of this noisy image. By realizing many independent images in this way, edge effects and the effect of particular configurations were minimized.

Figure 1

These circles illustrate examples in which each pixel intensity is initialized using a discrete uniform distribution with range of 1 to 100. An independent Poisson random variable with parameter equal to this intensity is then realized for each pixel. We then used the filtered back projection (FBP) and expectation maximization algorithm (EM) to reconstruct this circle based on its projection, as is done for PET images.

Images were reconstructed from the simulated sensor data using the filtered back-projection (FBP) and expectation maximization (EM) algorithms. For the FBP reconstruction, the ramp filter and linear interpolation were used and the image was padded with zeros up to 126 voxels by 126 voxels. This reconstruction exactly mimics the canonical implementation of the algorithm. For the iterative EM reconstruction, the initial image had uniform intensity 1. The canonical full form formula was used for the updates of the EM. The A matrix was formed by calculating the explicit point spread function for all integer angles from 0 to 179 degrees (see Supplementary Material for algorithmic details). The pervasive OSEM algorithm is a subset of the EM algorithm that substantially decreases the computational load of reconstruction; therefore, all results shown for the EM algorithm generalize to OSEM.

Due to the high spatial frequency in the focus of the image, 300 iterations were used for each EM reconstruction. The image did not appear qualitatively different after 20 iterations. The magnitude of variance was also observed to decrease asymptotically with iteration number (data not shown). The choice of 300 was made to maximize the potential for high spatial frequency noise that may better match the underlying data.

All simulations were conducted in MATLAB 7.14 (Mathworks), and all statistical analysis was conducted in R (see Supplementary Material: Regression of Trends in Reconstructed Voxel Intensity). Signal to noise ratio was calculated as the ratio of original intensity to the standard deviation of the reconstructed intensity. This is equivalent to a hypothesized two-fold change in original intensity.

Results

A detailed statistical analysis of the reconstructed images reveals important trends. Figure 2 illustrates the probability density of the reconstructed values with respect to their initial intensity. In this figure, all densities above 0.1 are rounded down to 0.1 to facilitate comparisons between the distributions. In the right panel, the Poisson nature of the original image is evident: the spread increases linearly with respect to the original intensity. The probability densities of the reconstructed intensities are markedly different from that of the original image. For both algorithms, the variance is much more homogenous and more extreme values shrink toward the mean. In order to formalize these observations, we fit statistical models to these trends. All intervals below reflect 95 percent confidence intervals.

Figure 2

This figure illustrates the probability distribution of the reconstructed voxel intensity for each of the reconstruction algorithms. For comparison, the right panel illustrates the original probability distribution before reconstruction.

First, we address the observation that all intensities shrink back toward the mean, albeit less so for EM than FBP (Figure 3). This regression back to the mean appears to be linear for FBP reconstructed voxel intensities (FBP-RVI) and quadratic for EM reconstructed voxel intensities (EM-RVI). This quadratic trend results in fitting high intensity voxels more closely compared to low intensity voxels. We fit statistical models to quantify and compare these trends across reconstruction algorithms. The FBP reconstructed voxel intensities (FBP-RVI) regressed back to this mean linearly with slope of -0.64 and intercept of 31 units (-0.0642 to -0.0634 and 31.56 to 31.63). The EM reconstructed voxel intensities (EM-RVI) regressed quadratically back to the mean with acceleration of 0.0018 units-1 (0.00176 to 0.00179). After controlling for this quadratic term, the EM-RVI had a slope of 0.038 closer to zero and a 7 units smaller intercept than the FBP-RVI (-7.44 to -7.37 and 0.037 to 0.040). The F statistic of this composite model was 3.7 million with 4 and 195 degrees of freedom, resulting in a model-wide p value of less than 10-16. There was no evidence that the residuals deviated from the assumption of independent identically distributed Gaussians. Even though the EM algorithm converges quadratically to the maximum likelihood solution [22,23], calculating more iterations does not significantly change any of these fitted parameters.

Figure 3

This figure illustrates the magnitude of the reconstructed intensity bias of each of the algorithms. The line thickness represents the standard error for each point. This standard error is small due to the large sample size. The FBP is indicated by cyan and the EM is indicated by green.

The focus of this report is the signal to noise ratio of reconstructed voxels. Controlling for the biases addressed above, the signal to noise ratio to detect a hypothesized two-fold change in intensity was substantially reduced for EM reconstructions compared to FBP reconstructions (Figure 4). The maximum signal to noise ratio for FBP-RVI was 60 percent larger than that of the maximum for EM-RVI. For both algorithms, this original intensity dependent increase in the signal to noise ratio with respect to original intensity reflects similar trends in variance as seen in the bias. The FBP-RVI variance increases linearly with intercept of 73 units2 and slope of 0.013 units (73.36 to 73.82 and 0.009 to 0.017). The EM-RVI variance increased quadratically with acceleration of 0.004 (0.0040 to 0.0044). After controlling for this quadratic term, the EM-RVI variance had a 13 units2 smaller intercept and a 1 unit larger slope (-13.6 to -12.8 and 1.04 to 1.08). This means that the EM performs slightly better for extremely low intensity voxels, but variance in EM-RVI and FBP-RVI also quickly increases as original intensity increases. The F statistic of this composite model was 280,000 with 4 and 195 degrees of freedom, resulting in a model-wide p value of less than 10-16. There was no evidence that the residuals deviated from the assumption of independent identically distributed Gaussians.

Figure 4

This figure illustrates the magnitude of the signal to noise ratio of the reconstructed intensity each of the algorithms. Signal to noise ratio was calculated as the original intensity divided by the standard deviation of the biased reconstructed intensity. This corresponds to a hypothesized two-fold change in intensity. The line thickness represents the standard error for each point. This standard error is small due to the large sample size. The FBP is indicated by cyan, and the EM is indicated by green.

Discussion

These striking results have a profound impact on the interpretation of PET images using quantitative and visual measures. We demonstrated in our simulations that PET is insensitive to all but large scale changes in hypometabolic regions. Therefore, we caution against the interpretation of hypometabolic lesions when reading PET images both visually and quantitatively. However, we confirm that EM improves the spatial resolution of reconstructed images by decreasing the bias introduced by nearby voxels when compared to the FBP, but we also illustrate that this bias correction results in a substantial decrease in the signal to noise ratio. Consequentially, even though EM reconstructed images are more consistent with our knowledge of the underlying biological structures, this increased spatial resolution comes at the cost of decreased statistical power of quantitative measures of the signal.

When interpreting PET images, our results suggest that one should focus on regions that are normally hypermetabolic relative to the surrounding tissue and caution against interpretation of changes in hypometabolism. For example, this is particularly important when interpreting images from patients with tuberous sclerosis for identification of epileptic foci. The tubers that characterize this disease can be small and distributed throughout the brain. Lee & Salamon suggests that hypometabolic lesions corresponding with structural abnormalities are candidates for epileptic foci [1]. If a structural lesion is in a hypometabolic region, however, our results suggest that there is very little power to detect metabolic abnormalities. This has the potential to increase the false negative rate for foci detection and thereby lead to patients with multifocal epilepsy being diagnosed with focal epilepsy. This misclassification can lead to patients undergoing focal surgical treatment that fails to control their seizures [13-15].

However, this does not hinder the ability of PET to recognize changes in relatively hypermetabolic tissue. It is important when reading PET images for one to consider the expected metabolism in the region of interest. If the expected metabolism is high, then most observed changes are interpretable and clinically relevant. Conversely, if the expected metabolism is low, then one should recognize that only comparatively large changes in metabolism are interpretable.

These findings also provide further motivation for the development of focused radioactive PET tracers to improve sensitivity [24]. Focused tracers target particular receptors or tissue types. For example, in Parkinson’s disease, there is increased neural death in the substantia nigra, resulting in decreased metabolism [25]. As we have shown, the power to detect these subtle, highly localized hypometabolic lesions is limited with FDG-PET. Our results suggest that the radioactive serotonin analog, 18F-DOPA, that differentially localizes to the substantia nigra in normal tissue has a stronger signal to noise ratio compared to FDG [26]. If the relative localization of this tracer is reduced, this may provide early diagnostic or more detailed prognostic information for the patient [27]. From a research perspective, this early detection could result in the development of novel pharmaceutical intervention that could slow the progression of disease. This also suggests that PET experiments will have higher signal to noise ratios if they are designed such that they focus on changes in tissue that is the target of the tracer.

The implication of these findings is particularly salient for quantitative PET analysis that has the potential to capture more subtle or distributed trends in metabolism. Conventional analysis of PET segments the brain into focused regions of interest and averages the reconstructed metabolic rate across the entirety of each region [28,29]. Although it is tempting to suggest that this averaging improves the signal to noise ratio with respect to the factors we have modeled, this is, unfortunately, not the case. Instead, the linearity of the noise spreading suggests that the variance from hypermetabolic voxels spreads across the entire projection and is not corrected by reconstruction algorithms. This suggests that the noise across a local hypometabolic region is correlated. Because of this correlation, the average then estimates the value of the signal plus the noise instead of separating the two. Therefore, our results suggest that PET is systematically insensitive to the detection of changes in hypometabolic tissue even when averaged over lower resolution regions of interest.

Our guidance to bias against interpretation of changes in hypometabolic tissue, however, is most generalizable to high resolution changes in metabolism. If these changes are widely distributed over tissue, then the number of hypermetabolic voxels that contribute to each projection decreases. Consequently, the noise is dominated only by the voxel with highest metabolism within the larger region. This is especially relevant to current analysis of epileptic foci because, due to the low spatial resolution of surgical procedures, only large magnitude, low resolution changes are clinically meaningful.

The substantial decrease in signal to noise ratio caused by the EM reconstruction suggests that while it vastly outperforms the FBP in manual interpretability and spatial resolution, EM may be not be appropriate for quantitative analysis of the PET signal. Based on our results, studies based on EM reconstructed images need 2.5 more patients or images than studies based on FBP reconstructed images to achieve the same signal to noise ratio (see Supplementary Material). Computation time for both reconstructions is relatively inexpensive; therefore, our results suggest that both methods should be employed. The EM reconstruction should be used for visual interpretation and the FBP reconstruction should be used for statistical inferences. This guidance, however, is balanced by the fact that FBP reconstructed images are more biased than EM reconstructed images. Using the mean squared error, which incorporates both bias and variance, the signal to noise ratio for the EM remains less than the FBP for the majority of voxel intensities.

It is unclear if alternate reconstruction algorithms not studied here address the limited signal to noise ratio of hypometabolic voxels. Our recommendation that manual image analysis should be modulated by the expected metabolism can be incorporated into reconstruction algorithms by using a Bayesian prior as was done in [30]. This prior would decrease the variance of the reconstructed image, but it may introduce further bias against detecting hypometabolic lesions that are incommensurate with the prior. Additionally, in each of the projections from which images are reconstructed, the magnitude of the noise from hypermetabolic voxels likely obscures any signal changes in hypometabolic voxels, as we illustrate and discuss here. The improved recovery of this and other information is a major motivation for the continued efforts to improve reconstruction techniques.

One could naturally suggest acquiring multiple PET images from the same patient to better quantify the noise distribution, but this practice is limited by expense. Statistically, one can expect that collecting multiple samples will increase the signal to noise ratio by a multiplicative factor of the square root of the number of samples. For hypometabolic voxels, however, the signal to noise ratio is so low that hundreds of PET images would be insufficient to reveal relatively large changes. Each PET, however, has substantial cost in physician, scanner, and patient time and resources. Simply splitting each scanning session into smaller time windows also does not solve the problem because the spatial resolution is a function of the total number of positrons observed [20,21]. The acquisition and incorporation of three-dimensional projections of the PET data could, however, increase the number of informative samples per image to improve the signal to noise ratio [31,32]. For extremely hypometabolic lesions, however, we expect that the inherent signal to noise ratio is so low that the signal is not recoverable through these techniques.

The knowledge that the noise in reconstructed space is likely heteroscedastic can be incorporated into statistical models of the signal in two theoretically equivalent ways. From a frequentist perspective, this can be done by relaxing the assumption that the residuals of the model are identically distributed. Instead, the variance of the residuals can be modeled as a linear or quadratic function of signal strength. By modeling this source of noise, therefore removing its contribution to the standard error of the model, we expect that the fit of the model would increase [33-35]. From a Bayesian perspective, one could introduce a prior that linearly or quadratically de-weights the contribution of hypometabolic regions. This de-weighting may also help ill-posed models, like those used in machine learning, to reduce their propensity to over fit the data by incorporating additional knowledge. As shown recently by Chu et al., this incorporation of additional biological and physical information may result in improved predictive performance [36].

These simulated results can be extended to address the signal to noise ratio in specific regions of interest. In particular, this approach of simulating the Poisson noise can be used to determine the sensitivity of FDG-PET to detect differences in numerous regions of interest in the brain. This could be used to give a more detailed explanation of the power of PET to describe high resolution metabolic changes. This could lead to an improved interpretability of smaller magnitude changes that indicate subtle phenomena. In particular, these subtle changes could be used in the aging population to predict which patients will progress to Alzheimer’s disease, as is currently being actively addressed using genotypic and MRI-based measures [37].

In order to maximize the interpretability of our results and reduce computational complexity, we made a number of simplifying assumptions. The regular shape and voxel intensities were chosen to improve the interpretability of our results. This is illustrated in the linear and quadratic trends in Figures 3 and 4. We do not expect that the discrete uniform distribution of voxels in the original image space to have any other effect on reconstructed intensities. This simplification resulted in a deeper understanding of the forces generating our results. The results can be easily generalized to the interpretation of a diverse set of hypometabolic lesions on cranial PET. The diversity of lesion location and type is prohibitively large to address in a single publication. Similarly, our results focus on the effect of Poisson noise from radiologic decay in FBP and EM reconstructed PET images. We hypothesize that similar limitations will be present in images reconstructed with other algorithms but more investigation is needed to explore this hypothesis.

We also ignore the effects of attenuation, randoms, scatter, dead time, detector normalization, scan length, decay, interpolation and the specific reconstruction filter because their inclusion does not influence our conclusions. These factors either uniformly increase the variance of the reconstructed intensities or exaggerate the contribution of hypermetabolic voxels to the total variance of the image. Even though these are two dimensional images, the concepts are readily generalizable to three dimensions (see Supplementary Material).