About this Journal Submit a Manuscript Table of Contents
Journal of Applied Mathematics
Volume 2013 (2013), Article ID 732178, 14 pages
Research Article

A Novel Adaptive Probabilistic Nonlinear Denoising Approach for Enhancing PET Data Sinogram

1Department of Electronics and Informatics (ETRO-IRIS), Vrije Universiteit Brussel, Pleinlaan 2, 1050 Brussels, Belgium
2Faculty of Information Technology and Computer Engineering, Palestine Polytechnic University, Hebron, Palestine
3Interuniversity Microelectronics Centre (IMEC), Leuven, Belgium

Received 23 March 2013; Revised 4 May 2013; Accepted 4 May 2013

Academic Editor: Hang Joon Jo

Copyright © 2013 Musa Alrefaya and Hichem Sahli. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.


We propose filtering the PET sinograms with a constraint curvature motion diffusion. The edge-stopping function is computed in terms of edge probability under the assumption of contamination by Poisson noise. We show that the Chi-square is the appropriate prior for finding the edge probability in the sinogram noise-free gradient. Since the sinogram noise is uncorrelated and follows a Poisson distribution, we then propose an adaptive probabilistic diffusivity function where the edge probability is computed at each pixel. The filter is applied on the 2D sinogram prereconstruction. The PET images are reconstructed using the Ordered Subset Expectation Maximization (OSEM). We demonstrate through simulations with images contaminated by Poisson noise that the performance of the proposed method substantially surpasses that of recently published methods, both visually and in terms of statistical measures.

1. Introduction

Positron Emission Tomography (PET) is an in vivo nuclear medicine imaging method that provides functional information of the body tissues. The PET image results from reconstructing very noisy, low resolution raw data, that is, the sinogram, in which important features are shaped as curved structures. Enhancing the PET image spurred a wide range of denoising models and algorithms. Some methodologies focus on enhancing the reconstructed PET image directly, where others prefer enhancing the sinogram prior to reconstruction.

Although developing an appropriate denoising method for filtering the PET images has been widely studied in the last two decades, this problem is still receiving high attention from researchers trying to improve the diagnosis accuracy of this image. Existing methods may suffer drawbacks such as the careful selection of a high number of parameters, smoothing of the important features boundaries, or prohibitive computation. Recently, nonlinear diffusion techniques have been investigated for PET images. Many researchers did explore the application of the well-known Perona and Malik anisotropic diffusion [1] in combination with diverse diffusivity functions on PET images [26], as well as on sinograms [79].

In [5], authors introduced a postreconstruction adaptive nonlinear diffusion (Perona and Malik) filter based on varying the diffusion level according to a local estimation of the image noise. Applying the nonlinear anisotropic filtering method on the whole body, PET image and the sinogram were presented in [4, 8, 10], respectively. Results showed that the anisotropic diffusion filtering algorithm helps improve the quantitative aspect of PET images.

In the study of [9], combining the anisotropic diffusion method with the coherence enhancing diffusion method for filtering the sinogram as a preprocessing step was proposed. However, the considered cascading approach is time consuming, and the results are highly dependent on the parameters selection criteria. Zhu et al. [11] built the diffusivity function using fuzzy rules that were expressed in a linguistic form. The mean curvature and the Gauss curvature diffusion algorithms for filtering the PET images during reconstruction were investigated in [12]. An anatomically driven anisotropic diffusion filter (ADADF) as a penalized maximum likelihood expectation maximization optimization framework was proposed in [2]. This filter enhances a single-photon emission computed tomography images using anatomical information from magnetic resonance imaging as a priori knowledge about the activity distribution. Authors in [3] proposed a reconstruction algorithm for PET with thin plate prior combined with a forward-and-backward diffusion algorithm. The main drawback of the above filter, with respect to sinogram images is that the diffusion produces important oscillations in the gradient. This leads to a poorly smoothed image [11, 12]. Moreover, the adopted diffusivity functions do not consider the special properties of the sinogram, which are high level of Poisson noise and curved-shape features.

Happonen and Koskinen [13] proposed filtering the sinogram in a new domain, that is, stackgram domain, where the signal along the sinusoidal trajectories can be filtered separately. They applied this method using the Gaussian and nonlinear filters. Radon transform is applied for inverse mapping of the filtered data from the stackgram domain to the sinogram domain. The above described method is impractical for noise reduction of the medical images such as PET, since the corresponding 3D stackgram requires an enormous space of computer memory. Furthermore, it is a complex method not suitable for clinical applications where timely reconstructien of the PET image is a very important issue.

Most of the above studies considered a global image noise level. Local or varying noise level has been used in [14], where a nonlinear diffusion method for filtering MR images with varying noise levels was presented. The authors assumed that the MR image can be modeled as a piecewise constant (slowly varying) function and corrupted by additive zero-mean Gaussian noise. Pizurica et al. [15] proposed a wavelet domain denoising method for subband-adaptive and spatially adaptive image denoising approach. The latter approach is based on the estimation of the probability that a given coefficient contains a significant noise-free component called “signal of interest.” The authors of [15] found that the spatially adaptive version of their proposed method yielded better results than the existing spatially adaptive ones.

In this work, based on the following PET sinogram characteristics: (i)the important features in the sinogram are curved structures with high contrast values and these represent the regions of interest in the reconstructed PET image, for example, tumor, (ii)the weak edges in the sinogram are the edges that contain low contrast values, (iii)the noise in the sinogram is a priori identified as a Poisson noise, we propose filtering the sinogram by means of a curvature constrained filter. The amount of diffusion is modulated according to a probabilistic diffusivity function that suits images contaminated with Poisson noise, being the known noise distribution of PET sinograms. Further, considering the singoram characteristics, we propose a probabilistic edge-stopping function based on Chi-square prior for the ideal sinogram gradient with a spatially adaptive algorithm for calculating the prior odd at each pixel. We show that this method is improving and denoising sinogram data which leads to enhanced reconstructed PET image.

The remainder of the paper is organized as follows, Section 2 briefly reviews the notions of curvature motion, edge affected diffusion filtering, and self-snakes. The proposed adaptive filtering scheme is introduced in Section 3. Finally, Section 4 discusses the experimental results, while the conclusions and future work are given in Section 5.

2. Geometry Driven Scale-Space Filtering

This section reviews the formulations for Mean Curvature Motion (MCM), Edge Affected Variable Conductance Diffusion (EA-VCD) and self-snakes. Also we recall the Gauge Coordinates notions. Extensive discussion can be find in [16]. Let be a scalar image defined on the spatial image domain , then the family of diffused versions of is given by where is referred to as the scale-space filter, is denoted by the scale-space image, and the scale [6]. The denoised or enhanced version of is a given that is the closest to the unknown noise-free version of . In the following we denote by .

2.1. Curvature Motion

Mean Curvature Motion (MCM) is considered as the standard curvature evolution. MCM allows diffusion solely along the level lines. In Gauge coordinates (see Section 2.4) the corresponding PDE formulation is ( is Euclidian curvature) as follows: Hence diffusion solely occurs along the -axis.

2.2. Edge Affected Variable Conductance Diffusion

Variable Conductance Filtering (VCD) is based on the diffusion with a variable conduction coefficient that controls the rate of diffusion [16]. In the case of Edge Affected-VCD (EA-VCD), the conductance coefficient is inversely proportional to the edgeness. Consequently it is commonly referred to as the edge-stopping function , in which the edgeness is typically measured by the gradient magnitude. The EA-VCD is governed by The above PDE system together with the initial condition given in (1) is completed with homogenous von Neumann boundary condition on the boundary of the image domain. Note that the Perona and Malik’s antitropic diffusion [1] is an EA-VCD.

2.3. Self-Snakes

Self-snakes are a variant of the MCM where an edge-stopping function is introduced [17]. The main goal is preventing further shrinking of the level-lines once they have reached the important image edges. For scalar images, self-snakes are governed by This equation adopts the same boundary condition as (3). Furthermore, it can be decomposed into two parts [16, 17]: The first part describes a degenerate forward diffusion along the level lines that is orthogonal to the local gradient; it allows preserving the edges. Additionally, the diffusion is limited in areas with high gradient magnitude and encouraged in smooth areas. Actually the first term is the constraint curvature motion. The second term can be viewed as a shock filter since it pushes the level-lines towards valleys of high gradient, acting as Osher’s shock filter.

2.4. Gauge Coordinates

An image can be thought of as a collection of curves with equal value, the isophotes. At extrema an isophote reduces to a point, at saddle points the isophote is self-intersecting. At the noncritical points; a Gauge coordinate system [18] is defined by two orthogonal axes and , which correspond to the directions of the tangent and normal at the isophote. The second order Gauge derivatives of the image in the and directions have the following second-order structures:

These gauge derivatives can be expressed as a product of gradients and a matrix with second-order derivatives [18]:

For the matrix equals the Hessian . For it is det . Note that the expressions are invariant with respect to the spatial coordinates. The above expression can also be obtained as follows [16]: where is given as and .

The two expressions of (6) can be combined linearly in a PDE setting. In this way, an image is evolved according to a weighted combination of these two invariants: with .

Equation (9) comprises a diffusion modulated by along the image edges (a smoothing term) and a diffusion adjustable by across the image edges (a sharpening term). Careful modeling of these terms allows efficiently denoising the PET sinograms, whilst keeping their interesting features.

3. The Probabilistic Curvature Motion Filter

The proposed probabilistic curvature motion filters are based on the idea of the probabilistic diffusivity function [19], where the diffusivity function is expressed as the probability that the observed gradient presents no edge of interest under a suitable marginal prior distribution for the noise-free gradient. In [19], the probabilistic diffusivity function has been defined as where the normalizing constant is set to to ensure that , the hypothesis describes the notion whether an edge element of interest is present given the considered noise and an edge element of interest is absent. For a noisy gradient model , we set with being the ideal, noise-free, gradient magnitude, and being the standard deviation of the noise . In [19] it has been demonstrated that with the prior odds defined as and the likelihood ratio Considering a Laplacian prior , we have [19]: and the parameter can be estimated as with denoting the variance of the noisy image and as defined above. Due to limited space, the reader is referred to [19] for the detailed expression of in (14).

It has to be noted that the diffusivity function (12) fits in the cluster of backward-forward diffusivities, and it has no free parameters to be set. Moreover, for the considered PET sinograms, the noise standard deviation, , in (16) is estimated as , where the image noise, , is reconstructed via wavelet decomposition of , from the two finest resolution levels coefficients, using Daubechies (4) wavelet.

3.1. Probabilistic Self-Snakes (PSS)

It can be demonstrated that the diffusion of scalar images via EA-VCD can be decomposed into (5), which can be rewritten as follows [16]: consolidating the properties of both the self-snakes and the EA-VCD into a single diffusion schema.

Considering (9) and the probabilistic diffusivity function, if we set , and the sharpening term, , we obtain the probabilistic self-snakes (PSS) [7]. By its nature, the PSS enhances the weak edges and/or features in the sinogram. The second term in (32) allows the sharpening. In this way, weak but important edges are enhanced whilst the noise is removed efficiently.

PSS proved to be very effective and flexible for the sinogram image where the high contrast regions, which represent a tumor in the reconstructed PET, should be smoothed wisely without blurring the poor edges [7]. Like EA-VCD, the main advantage of this filter is that the average gray value of the image is not altered during the diffusion process which is a significant issue in the sinogram.

3.2. Adaptive Probabilistic Diffusivity Function Based on a Chi-Square Prior for a Sinogram

The probabilistic diffusivity function (12) does not take into account the spatially varying noise levels such as the case for the sinograms. It has a global threshold parameter, which is related to the image noise standard deviation . Hence, in such formulation, if two pixels/voxels have equal gradient magnitude, they will have the same values, no matter the noise level at these pixels. In this work, the probabilistic diffusivity function is improved by considering the local statistical noise at each element. We adopt the ideas of [15] and adapt the estimator to the local spatial context in the image.

3.2.1. Spatially Adaptive Probabilistic Diffusivity Function

The most appropriate way to achieve a spatial adaptation is to estimate the prior probability of signal presence adaptively for each pixel instead of fixing it globally. This can be achieved by conditioning the hypothesis on a local spatial activity indicator such as the locally averaged magnitude or the local variance of the observed gradient.

To estimate the probability that “signal of interest” is present at the position , we consider a local spatial activity indicator, , at each position. is defined as the locally averaged gradient magnitude within a relatively small window, , with a fixed size, , around the position in the gradient image.

Starting from the prior odd, (13), we replace the ratio of “global” probabilities by a locally adaptive prior, ; that is, and are conditioned on the local spatial indicator where () is the local likelihood ratio:

Considering Bayes’ rule, the probability that an “edge of interest” is present at position , , is given as The spatially adaptive probabilistic diffusivity function can then be formulated as with We ensure that , because the minimum of is at and thus peaks at .

Intuitively, the proposed method considers an “observed gradient” at a given location as how probable this location presents useful information compared to its neighborhood, based on(1)the likelihood ratio via ,(2)a measurement from the local surrounding via ,(3)the global prior odd via .

The local spatial activity indicator is defined as where is the gradient magnitude at location . Assume that all the elements within the small window are equally distributed and conditionally independent. With these simplifications, the conditional probability of given in a window of size , which is denoted as , is given by convolutions of with itself as follows: while the conditional probability of given is given by convolutions of with itself:

Due to limited space, the reader is referred to [19] for the detailed expression of , , and .

Figure 1 illustrates the original probabilistic diffusivity function and adaptive probabilistic diffusivity function. The adaptive function has lower values, as shown in Figure 1(a), which allows better preservation of important edges than the original function. The negative peaks of the proposed function occur usually at larger gradient magnitude than the original one, which indicates stronger edge enhancement capability. This property indicates a quick smoothing of the nearly uniform areas while maintaining and enhancing the strong edges.

Figure 1: The Adaptive probabilistic diffusivity function versus the original probabilistic diffusivity function . (a) Diffusivity functions along the edge (). (b) Diffusivity functions across the edge ().
3.2.2. A Chi-Square Prior for Noise-Free Sinogram Gradient

The noise in the sinogram is a priori identified as a Poisson noise [12, 20]. The magnitude of Poisson noise varies across the image, as it depends on the image intensity. This makes removing such noise very difficult. In the original probabilistic diffusivity function (Section 3), the Laplacian prior was imposed for the ideal image gradient that is contaminated with Gaussian noise. In order to take into account the sinogram’s noise distribution in the filtering scheme, in this section we aim at redefining the diffusivity function for handling Poisson noise. This can be accomplished by finding an appropriate prior for the ideal noise-free sinogram gradient.

In the following, we argue that we can represent Poisson distribution by a Gaussian distribution as . Afterwards, we demonstrate that the gradient magnitude of the noise-free sinogram follows a Chi-square distribution, and finally we reformulate the probabilistic diffusivity function based on a Chi-square prior for the noise-free sinogram gradient.

In the literature, several studies demonstrated that the Poisson distribution (of probability mass , with being the number of occurrences of an event and being the expected number of occurrences during a given interval) approaches a Gaussian density function in the case of high number of counts [21, 22]. Moreover, Miller et al. [21] showed that the Gaussian approximation is surprisingly accurate, even for a fairly small number of counts. To illustrate this, we use the logarithmic function to simplify the proof: Using Stirling’s formula (for large as we are assuming here) , we have

Assuming that the sinogram gradient is approximated by absolute difference of neighboring pixel values on a 2-connected grid, we demonstrate that the gradient of Poisson random variables follows a Skellam distribution. Then, we show that the Skellam distribution can be approximated as a Gaussian distribution.

Let and be two statistically independent adjacent pixels in the observed sinogram be with a Gaussian distribution and , respectively. The distribution of the difference, , of two statistically independent random variables and , each having Poisson distributions with different expected values and , is denoted as the Skellam distribution [23] and can be given as where is the modified Bessel function of the first kind.

The difference between two Poisson variables has the following properties: (i) and (ii) . Considering these properties, the cross-correlation, and the delta function, the approximated distribution of the sinogram gradient can be given as .

Based on the assumption that the sinogram gradient follows the distribution, we can show that the distribution of this gradient leads a Chi-square distribution as follows:

The Chi-square distribution is defined by the following probability density function: where denotes the Gamma function and is a positive integer that specifies the number of degrees of freedom. For the noise gradient model , the Chi-square with degrees of freedom (2 degrees because we are dealing with 2D images) is given as

Based on the above, the prior odd, (15), can be reformulated considering Chi-square prior instead of Laplacian prior and the noise . Note that the Chi-square with 2 degrees of freedom is almost a special case of Laplacian prior with a rate parameter (scale parameter) . This parameter determines the “scale” or statistical dispersion of the probability distribution. It indicates that the distribution of the ideal gradient is independent of a rate parameter since the Poisson noise is a pixel dependent (i.e., depends on the number of counts). Therefore, it is more natural to use Chi-square prior for estimating the ideal gradient of Poisson data rather than Laplacian prior with a parameter, , which is based on the variances of the image and noise gradients as indicated by (16).

In Figure 2, we illustrate the estimation of the Laplacian and Chi-square priors for the noise-free gradient. Figures 2(a) and 2(b) show the histograms of the gradient magnitudes for the noise-free and noisy images, respectively. The noise-free gradient histogram is typically sharply peaked at zero since the noise-free images typically contain large portions of relatively uniform regions that produce negligible gradient values. Sharp edges and textured regions produce some relatively large gradients, building in this way long tails of the gradient histogram. From the noisy histogram of Figure 2(b), we estimate the parameter of the prior for the noise-free sinogram gradient. The results are illustrated in Figure 2(c), which shows the estimated Laplacian (dotted curve) and Chi-square (continues curve) priors in comparison to the ideal, noise-free histogram.

Figure 2: Chi-square prior versus Laplacian prior. (a) Histogram of the noise-free gradient. (b) Histogram of the noisy gradient magnitude. (c) Laplacian and Chi-square priors for the noisy-free gradient, estimated from the noisy data.
3.2.3. Algorithm: The Adaptive Probabilistic Curvature Motion Filter Based on Chi-Square Prior

In summary the general equation of the proposed Adaptive Modified Probabilistic Curvature Motion (AMPCM) filter is given by Let , as given by (21), and the sharpening term ; the overall algorithm is given as shown in Algorithm 1.

Algorithm 1. Adaptive Modified Probabilistic Curvature Motion Filter Algorithm(i)Build the Probabilistic Diffusivity Function:(a) compute the prior odd based on the Chi-Square prior: (ii)Build the spatially adaptive diffusivity function  , for each pixel  : (a) compute the local spatial activity indicator (b) compute the likelihood ratio for each window: (c) compute the normalizing constant : (d) compute as defined by (20) (iii)compute the diffusivity function  (iv)Filter the sinogram,  , based on (8) recursion, for , (number of iterations 1):

4. Experiments and Discussion

For the analysis of the proposed denoising approaches, we use a simulated thorax PET phantom, containing three hot regions of interest (tumors) of activities 1.18, 1.8, and 1.57. realizations (noisy sinograms) with added noise of coincident events have been generated. Each sinogram has a size of pixels and its spacing is , with 128 detectors and 128 projection angles. Figure 3(a) illustrates the noise-free sinogram, and Figure 3(b) illustrates a noisy realization.

Figure 3: Experimental dataset. (a) Noise-free simulated sinogram. (b) Noisy sinogram. ((c)-(d)) Corresponding reconstructions.

For reconstructing the PET images, we adopt the Ordered Subset Expectation Maximization (OSEM) reconstruction algorithm [24]. One way to represent the imaging system is with the following linear relationship: where is the set of observations (sinogram), is the known system model (projection matrix), is the unknown image, and is the noise. The goal of reconstruction is to use the data values (projections through the unknown object) to find the image . Under the Poisson assumption, counts of all the detectors around the object are considered as independent Poisson variables. OSEM groups projection data into an ordered sequence of subsets and progressively processes every subset of projections in each iteration process [24]. The OSEM method results are highly affected by the number of iterations and subsets used. The iterations should be ended before the noise becomes too dominant and not too early to avoid losing important information. In our experiments, the parameters of the OSEM algorithm are set as follows: we use subsets and run them for iterations. Such PET reconstructions are illustrated in Figures 3(c) and 3(d) for the noise-free and noisy sinograms, respectively.

For the experimental assessment of the proposed diffusion filtering, we use two sets of evaluation measures. The first set stems from measuring the quality of the filtering techniques on the sinogram, whilst the second set originates from validating the quality of the PET reconstruction, after filtering the sinogram observations. As ground-truth information, the former uses the noise-free sinogram, while the latter needs prior identification of the important areas by a medical professional.

A fundamental issue with scale-spaces induced by diffusion processes, as the ones proposed in this paper, is the automatic selection of the most salient scale. For PET sinogram denoising application, we use an earlier proposed optimal scale selection approach [25], where the maximum correlation method has been adopted: with being the so-called outlier noise estimated using wavelet-based noise estimation. is the zeroth scale; thus and represents the initial amount of noise.

4.1. Sinogram Denoising Evaluation

In this work, we are interested in comparing the proposed AMPCM filter with filtering approaches from the literature: the Probabilistic self-snakes (PSS) [7], Perona and Malik filter (P-M) [1], and the Noise-Adaptive Nonlinear Diffusion Filtering (NAF) [14].

The Lorentzian diffusivity function is used for applying the filter. This function was proposed by Perona and Malik [1] and widely used in the literature. The contrast parameter is calculated based on the value given by a percentage of the cumulative histogram of the gradient magnitude [16]. The same diffusivity constant is used for all filters ().

Figure 4 illustrates the optimal enhanced scale of the sinograms for all the considered filtering methods. The filtered sinogram by AMPCM and PSS shows better enhanced results especially the curvy shape features.

Figure 4: Optimal enhanced scale () of the noisy sinogram of Figure 3(b).  

For the qualitative assessment of the denoised sinogram, , with respect to the noise-free image , we adopt the following measures [25]:DQ1 The Peak Signal to Noise Ratio (PSNR) is a statistical measure of error, used to determine the quality of the filtered images. It represents the ratio of a signal power to the noise power corrupting. Obviously, one sees that the higher the PSNR, the better the quality: DQ2 The correlation () between the noise-free and the filtered image. The higher this correlation is, the better the quality is: DQ3 The calculated variance of the noise (NV) describes the remaining noise level. Therefore, it should be as small as possible: We are interested in comparing the maximum of each measure for different filtering approaches.

Table 1 lists the above measures for each of the considered filtering approaches. The best performance (maximum of the measure) is displayed in bold. As it can be seen the best performing filtering is achieved when using the AMPCM and PSS. Furthermore, we notice that for all measures, AMPCM and PSS outperform the NAF and (P-M) filters. The performance of the smoothing algorithms proposed in this paper is remarkable in discriminating a true edge from image noise (see Figure 4). We also notice the improved performance of the probabilistic curvature motion algorithms as compared to the standard diffusion algorithms. Table 2 gives the number of iterations to reach the optimal scale as defined in (40).

Table 1: Denoising quality measures.
Table 2: Optimal number of Iterations.

The absolute difference between the mean of the filtered sinograms and the noise-free one is another indication of the behavior and stability of the filtering methods. The mean results of the filtered realizations by the PSS and AMPCM filters are close to the ideal image, as shown in Figure 5. Comparing the mean results between the AMPCM and PSS with the other methods clearly demonstrates the effect of the sharpening/enhancing term which yields a better enhanced edges. On the other hand, P-M keeps edges unsmoothed, which is appearing as sharper edges in Figure 5. Figure 5 shows more noise remained, represented as larger values in the absolute difference images of P-M and NAF. Furthermore, using AMPCM and PSS, the absolute difference approaches zero (black image), indicating that the noise has been effectively and adaptively reduced from the noisy sinograms.

Figure 5: Mean of the filtered sinograms (left column). Absolute difference between the mean of the filtered sinograms and the noise-free image (right column).

The heavy noise is clearly eliminated without destroying the texture (curves) by the probabilistic curvature motion filter. From the above and other various examples, we have observed that the proposed algorithm has ability to eliminate the Poisson noise. No stair-casing has been observed, nor severe blurring is introduced. Figure 4(c) shows that P-M filter performs well for relatively low levels of noise, while it results in poor quality of images for high noise levels. However, Perona’s operator does not work well when applied to very noisy images because the noise introduces important oscillations of the gradient.

4.2. PET Reconstruction Evaluation

In this section, we use the smoothed sinograms, , for PET reconstruction via the OSEM algorithm. The reconstructed PET image is denoted as . For evaluating the effect of sinogram predenoising on the PET reconstruction, we use the contrast recovery method which aims to measure the quality of by measuring the contrast recovery of the interesting objects (ROIs) in the image. The contrast recovery is computed based on the contrast gain. The latter measures how much the ROIs (i.e., tumor) in the PET image are discriminated from the background by sharp edges and on the variance of the contrast gain for different realizations. Further, the contrast gain evaluates the stability of the applied algorithm. The contrast recovery curve is estimated using the set of regions of interest that were identified by a medical professional. This curve is widely used in the literature for evaluating the quality of the reconstructed PET images [26].

In Figure 3(c), there are white spots that represent tumors. We calculate the contrast gain for each realization ( denotes the number of realizations), and its overall variance . Let ( in our case) be the set of identified ROIs, and let be a representative background tissue area then where denotes a pixel, the scale number, and the contrast mean, . represent the mean intensities inside background, defined as .

Plotting the contrast gain in function of the variance given the smoothness parameter, which is in this work, the scale parameter yields the contrast recovery curve [26].

In order to perform, under the same conditions, the comparison of the contrast curves of the different diffusion scheme, one should identify the scales, , emanating from different diffusion schemes, having similar information content. For this, we use the scale synchronization proposed by Niessen et al. [27], being a Shannon-Wiener entropy used to compare nonlinear scale-space filters: where represents the averaging operator that maps an image onto a constant image, where each pixel equals the average feature vector.

The contrast gain and the variance were computed for 6 scales, for each filtering methods, with similar information content, selected based on Niessen et al. [27] approach. Figure 6 shows the contrast recovery curves for the considered diffusion schemes AMPCM, PSS, P-M, and CCM-Sapiro [17]. The best quality PET reconstruction is situated in the upper, that is, high contrast gain, left, that is, high stability, area of the plot. The figure shows that AMPCM has better contrast recovery of the interesting objects (ROIs) in the reconstructed PET images.

Figure 6: The contrast recovery curves for reconstructed PET by OSEM.

Figure 7 illustrates the reconstructed PET starting from the enhanced sinogram, by the AMPCM, along with the variance images of the reconstructed PET, the mean of the reconstructed PET images, and the absolute difference between the mean of the filtered PET and the noise-free image. We immediately notice that the background seems much flatter when adopting the OSEM reconstruction.

Figure 7: Evaluation of the PET reconstruction using enhanced sinograms.

Experiments results show that combining the probabilistic diffusivity function with the curvature motion diffusion produces a powerful nonlinear filtering method that is appropriate for the unique characteristics of the PET images. The AMPCM filter preserves the boundaries of the curvy shape features and wisely smooths the regions of interest as well as the other regions. Figure 7 demonstrates that the edge enhancement around image data is significantly improved. The effect of edge strengthening is even more pronounced for weaker edges in PET images or in image areas affected by a high level of noise.

5. Conclusions

Adaptive probabilistic curvature motion filter (AMPCM) for enhancing PET images is developed and discussed in this work. The filter is applied on the 2D sinogram pre-reconstruction. For considering the special characteristics of the sinogram data, a Chi-square is used as a marginal prior for noise-free sinogram gradient in the diffusivity function.

The qualitative evaluation results results show that, among other diffusivity functions, the probabilistic adaptive function seems more effective in smoothing all the homogenous regions that contain high level of noise. Further, AMPCM retain in an accurate way the location of the edges that defines the shape of the represented structures in the sinogram. It preserves the boundaries of the curvy shape features and wisely smooths the regions of interest as well as the other regions. The application of the proposed diffusion scheme results in well-smoothed images and preserves the edges. It gains the advantages of the curvature motion diffusion and the shock filter. Further, it deals better with the problem of poor and discontinuous edges which are common in PET sinograms.

Applied as a preprocessing step before PET reconstruction, we demonstrated via qualitative measures (the contrast curve, the variance, and the mean images) that the adaptive probabilistic diffusion function has a great capability and stability to detect and enhance the important features in the reconstructed PET image, which make it a reliable and practical approach as it has no free parameters to be optimized. All parameters are image based and are automatically estimated and proved to give the best results.


The authors would like to thank Professor M. Defrise, Division of Nuclear Medicine at UZ-VUB, for providing them the PET phantom images, the discussions, and feedback on the Poisson noise and Chi-square prior, as well as on the evaluation of the PET reconstruction results.


  1. P. Perona and J. Malik, “Scale-space and edge detection using anisotropic diffusion,” IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 12, no. 7, pp. 629–639, 1990. View at Publisher · View at Google Scholar · View at Scopus
  2. D. Kazantsev, S. R. Arridge, S. Pedemonte et al., “An anatomically driven anisotropic diffusion filtering method for 3D SPECT reconstruction,” Physics in Medicine and Biology, vol. 57, no. 12, pp. 3793–3810, 2012. View at Publisher · View at Google Scholar
  3. Z. Quan and L. Yi, “Image reconstruction algorithm for positron emission tomography with Thin Plate prior combined with an anisotropic diffusion filter,” Journal of Clinical Rehabilitative Tissue Engineering Research, vol. 15, no. 52, 2011.
  4. O. Demirkaya, “Post-reconstruction filtering of positron emission tomography whole-body emission images and attenuation maps using nonlinear diffusion filtering,” Academic Radiology, vol. 11, no. 10, pp. 1105–1114, 2004. View at Publisher · View at Google Scholar · View at Scopus
  5. D. R. Padfield and R. Manjeshwar, “Adaptive conductance filtering for spatially varying noisein PET images,” Progress in Biomedical Optics and Imaging, vol. 7, no. 3, 2006.
  6. J. Weickert, Anisotropic Diffusion in Image Processing, European Consortium for Mathematics in Industry, B. G. Teubner, Stuttgart, Germany, 1998. View at Zentralblatt MATH · View at MathSciNet
  7. M. Alrefaya, H. Sahli, I. Vanhamel, and D. Hao, “A nonlinear probabilistic curvature motion filter for positron emission tomography images,” in Scale Space and Variational Methodsin Computer Vision, vol. 5567 of Lecture Notes in Computer Science, pp. 212–2223, 2009.
  8. O. Demirkaya, “Anisotropic diffusion filtering of PET attenuation data to improve emission images,” Physics in Medicine and Biology, vol. 47, no. 20, pp. N271–N278, 2002. View at Publisher · View at Google Scholar · View at Scopus
  9. W. Wang, “Anisotropic diffusion filtering for reconstruction of poisson noisy sinograms,” Journal of Communication and Computer, vol. 2, no. 11, pp. 16–23, 2005.
  10. C. Negoita, R. A. Renaut, and A. Santos, “Nonlinear anisotropic diffusion in positron emission tomography kinetic imaging,” in SIAM Conference on Imaging Science, Salt Lake City, Utah, USA, 2004.
  11. H. Zhu, H. Shu, J. Zhou, C. Toumoulin, and L. Luo, “Image reconstruction for positron emission tomography using fuzzy nonlinear anisotropic diffusion penalty,” Medical and Biological Engineering and Computing, vol. 44, no. 11, pp. 983–997, 2006. View at Publisher · View at Google Scholar · View at Scopus
  12. H. Zhu, H. Shu, J. Zhou, X. Bao, and L. Luo, “Bayesian algorithms for PET image reconstruction with mean curvature and Gauss curvature diffusion regularizations,” Computers in Biology and Medicine, vol. 37, no. 6, pp. 793–804, 2007. View at Publisher · View at Google Scholar · View at Scopus
  13. A. P. Happonen and M. O. Koskinen, “Experimental investigation of angular stackgram filtering for noise reduction of SPECT projection data: study with linear and nonlinear filters,” International Journal of Biomedical Imaging, vol. 2007, Article ID 38516, 2007. View at Publisher · View at Google Scholar · View at Scopus
  14. A. A. Samsonov and C. R. Johnson, “Noise-adaptive nonlinear diffusion filtering of MR images with spatially varying noise levels,” Magnetic Resonance in Medicine, vol. 52, no. 4, pp. 798–806, 2004. View at Publisher · View at Google Scholar · View at Scopus
  15. A. Pizurica, P. Scheunders, and W. Philips, “Multiresolution multispectral image denoisingbased on probability of presence of features of interest,” in Proceedings of Advanced Concepts for Intelligent Vision Systems (Acivs '04), 2004.
  16. I. Vanhamel, Vector valued nonlinear diffusion and its application to image segmentation [Ph.D. thesis], Vrije Universiteit Brussel, Faculty of Engineering Sciences, Electronics and Informatics (ETRO), Brussels, Belgium, 2006.
  17. G. Sapiro, Geometric Partial Differential Equations and Image Analysis, Cambridge University Press, Cambridge, UK, 2001. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet
  18. A. Kuijper, “Geometrical PDEs based on second-order derivatives of gauge coordinates in image processing,” Image and Vision Computing, vol. 27, no. 8, pp. 1023–1034, 2009. View at Publisher · View at Google Scholar · View at Scopus
  19. A. Pižurica, I. Vanhamel, H. Sahli, W. Philips, and A. Katartzis, “A Bayesian formulation of edge-stopping functions in nonlinear diffusion,” IEEE Signal Processing Letters, vol. 13, no. 8, pp. 501–504, 2006. View at Publisher · View at Google Scholar · View at Scopus
  20. A. Teymurazyan, T. Riauka, H. S. Jans, and D. Robinson, “Properties of noise in positron emission tomography images reconstructed with filtered-back projection and row-action maximum likelihood algorithm,” Journal of Digital Imaging, vol. 26, no. 3, pp. 447–456, 2013. View at Publisher · View at Google Scholar
  21. G. Miller, H. F. Martz, T. T. Little, and R. Guilmette, “Using exact poisson likelihood functions in Bayesian interpretation of counting measurements,” Health Physics, vol. 83, no. 4, pp. 512–518, 2002. View at Scopus
  22. H. A. Gersch, “Simple formula for the distortions in a Gaussian representation of a Poisson distribution,” American Journal of Physics, vol. 44, no. 9, pp. 885–886, 1976. View at Publisher · View at Google Scholar
  23. J. G. Skellam, “The frequency distribution of the difference between two Poisson variates belonging to different populations,” Journal of the Royal Statistical Society A, vol. 109, no. 3, p. 296, 1946. View at MathSciNet
  24. X. Lei, H. Chen, D. Yao, and G. Luo, “An improved ordered subsets expectation maximization reconstruction,” in Advances in Natural Computation, vol. 4221, pp. 968–971, 2006.
  25. I. Vanhamel, C. Mihai, H. Sahli, A. Katartzis, and I. Pratikakis, “Scale selection for compact scale-space representation of vector-valued images,” International Journal of Computer Vision, vol. 84, no. 2, pp. 194–204, 2009. View at Publisher · View at Google Scholar · View at Scopus
  26. C. Comtat, P. E. Kinahan, J. A. Fessler et al., “Clinically feasible reconstruction of 3D whole-body PET/CT data using blurred anatomical labels,” Physics in Medicine and Biology, vol. 47, no. 1, pp. 1–20, 2002. View at Publisher · View at Google Scholar · View at Scopus