Abstract
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 [2–6], as well as on sinograms [7–9].
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.
(a) probabilistic Diffusion along the Edge
(b) Diffusion accross 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.
(a) Noise-free histogram
(b) Noisy gradient histogram
(c) Chi-square and Laplacian prior for the noise free histogram
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.
(a)
(b)
(c)
(d)
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.
(a) AMPCM
(b) PSS
(c) P-M
(d) NAF
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).
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.
(a) AMPCM
(b) AD-AMPCM
(c) PSS
(d) AD-PSS
(e) P-M
(f) AD-PM
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 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.
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.
Acknowledgment
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.