An efficient wavelet and curvelet-based PET image denoising technique

Abstract

Positron emission tomography (PET) image denoising is a challenging task due to the presence of noise and low spatial resolution compared with other imaging techniques such as magnetic resonance imaging (MRI) and computed tomography (CT). PET image noise can hamper further processing and analysis, such as segmentation and disease screening. The wavelet transform–based techniques have often been proposed for PET image denoising to handle isotropic (smooth details) features. The curvelet transform–based PET image denoising techniques have the ability to handle multi-scale and multi-directional properties such as edges and curves (anisotropic features) as compared with wavelet transform–based denoising techniques. The wavelet denoising method is not optimal for anisotropic features, whereas the curvelet denoising method sometimes has difficulty in handling isotropic features. In order to handle the weaknesses of individual wavelet and curvelet-based methods, the present research proposes an efficient PET image denoising technique based on the combination of wavelet and curvelet transforms, along with a new adaptive threshold selection to threshold the wavelet coefficients in each subband (except last level low pass (LL) residual). The proposed threshold utilizes the advantages of adaptive threshold taken from BayesShrink along with the neighborhood window concept. The present method was tested on both simulated phantom and clinical PET datasets. Experimental results show that our method has achieved better results than the existing methods such as VisuShrink, BayesShrink, NeighShrink, ModineighShrink, curvelet, and an existing wavelet curvelet-based method with respect to different noise measurement metrics, such as mean squared error (MSE), signal-to-noise ratio (SNR), peak signal-to-noise ratio (PSNR), and image quality index (IQI). Furthermore, notable performance is achieved in the case of medical applications such as gray matter segmentation and precise tumor region identification.

Appendices

Appendix A: Theoretical preliminary

Positron emission tomography plays a vital role in investigating functional analysis in the human brain. The functional activity of PET images should be defined as the distinguishable intensity difference between the objects and their neighborhood. Intensity distribution in the PET image depends on the applied radiopharmaceutical.

1.1 A.1 Noise characteristic in PET image

High noise and poor contrast in PET images can hamper various types of feature recognition processes that are mostly based on the intensity values. Noise characteristics in PET are not well-known till now, so sometimes well-known traditional noise removal techniques cannot perform well in PET image denoising. Typically, it is assumed that noise in PET is characterized as Poisson [5, 9, 10, 44] and mixed Gaussian-Poisson [38]. The effect of Gaussian noise, Poisson noise, and mixed Gaussian-Poisson noise in PET image are shown in Fig. 21 using histograms. For Poisson and mixed Gaussian-Poisson noisy PET image, the intensity scale is distributed with positive values only which are shown in Fig. 21c and d, respectively, whereas, the intensity distribution is seen to have negative values (Fig. 21b) with Gaussian noise which is additive in nature with zero mean and unit deviation. Gaussian noise and Poisson noise are treated as an uncorrelated and correlated component, respectively. Generally, the useful signal is represented as a mean value of that signal, whereas the standard deviation denotes the signal noise.

Gaussian noise

The probability density function of a Gaussian random variable is formulated as

$$ p_{G}(i) = \frac{1}{\sigma\sqrt{2\pi}}e^{-\frac{(i-\mu)^{2}}{2\sigma^{2}}}. $$

(25)

Fig. 21

Noise effect in Clinical PET image using histograms. a Histogram of above clinical PET image. b Histogram after adding Gaussian noise (σ = 60). c Histogram after applying Poisson noise. d Histogram after applying mixed Gaussian-Poisson noise (σ = 60)

In Eq. 25, i is the intensity level, μ and σ represent mean and standard deviation, respectively. Gaussian distribution is continuous, so a continuous noise can be modeled using the Gaussian distribution.

Let X{X_i,j, i, j = 1,2,…,N} be the noise-free image and Y{Y_i,j, i,j = 1,2,…,N} be the noisy image. If the noise type is Gaussian, the noisy image is generated by Eq. 26, where G(i,j) has a normal distribution N(0, 1) and σ is the noise standard deviation.

$$ Y(i, j) = X(i, j )+\sigma G(i, j) . $$

(26)

Poisson noise

The probability density [45] function of a Poisson [46] random variable for a time interval t is formulated as

$$ P(k)=\frac{e^{-\lambda t}(\lambda t)^{k}}{k!} . $$

(27)

Here k is a Poisson random variable, e denotes Euler’s number, and the average number of events occurring within an interval is denoted by λ. For Poisson distribution, Eq. 27 is also called probability mass function. In Poisson distribution, the occurrence of each event is independent with respect to other occurrences. A key characteristic of the Poisson distribution is that the variance is the mean. The Poisson distribution is discrete in nature, so a discrete noise can be modeled by using the Poisson distribution.

Gaussian and Poisson noise can be distinguished by the relationship between the mean of pixel intensity and the amplitude of the signal noise which is shown with a simulated phantom in Fig. 22. In the Gaussian noisy image, the amplitude of the Gaussian noise is almost constant throughout the image (Fig. 22e), whereas the amplitude of the Poisson noise is not constant (Fig. 22c) throughout the image because it is proportional to the pixel intensity within its neighborhood region.

Fig. 22

Noise effect in simulated Phantom. a Simulated noise-free image with eight homogeneous regions. b Poisson noisy image. c Mean and standard deviation (SD) of Poisson noisy image. d Gaussian noisy image. e Mean and standard deviation of the Gaussian noisy image. In both c and e, 3×3 neighborhood is chosen for measuring the mean and SD of each pixel

1.2 A.2 Wavelet

The wavelet transform [23, 47,48,49] corresponds to the decomposition of a quadratic integral function s(x) 𝜖L²(R) into a family of scaled and translated functions ψ_k,l(t) is shown in Eq. 28. The functions used in the wavelet transform are localized in the real and Fourier space.

$$ \psi_{k,l}(t)=k^{-1/2}\psi \left( \frac{t-l}{k}\right) . $$

(28)

The function ψ(x) is called wavelet function, which describes the band-pass behavior of the signal. The wavelet coefficient (d_k,l) is shown in Eq. 29, where ∗ refers to complex conjugate function, k ∈ $\mathbb {R}$⁺ and l ∈ $\mathbb {R}$.

$$ d_{k,l}=\frac{1}{\sqrt{k}}\displaystyle{\int}_{-\infty}^{\infty} S(x){\Psi}^{\ast}\left( \frac{x-l}{k}\right)dx . $$

(29)

Wavelet function must be orthogonal to its discrete translation that represents a mathematical condition such as in Eq. 29 called dilation equation where S is treated as a scaling factor.

1.3 A.3 Curvelet

Curvelet transform [50] provides multi-scale and multi-directional transformation with several features compared with wavelets for dealing with directional properties such as edges and curves. Candes et al. [50] decomposed an image into a series of subbands, then applied the concept of random transform and ridgelet transform [14,15,16,17, 51] on each band. The concept of continuous ridgelet function ψ_a,b,𝜃(y₁,y₂) was introduced by Cands [14, 50] and formulated as

$$ \begin{array}{@{}rcl@{}} \psi_{a,b,\theta}(y_{1},y_{2})=a^{-1/2}\psi((y_{1}\cos\theta + y_{2}\sin\theta-b)/a). \end{array} $$

(30)

This is oriented at angles 𝜃 and is constant along the lines y₁ cos 𝜃 + y₂ sin 𝜃=const. Here, a > 0 is the scale and b is the location. More details of curvelet are described in related Refs. [14,15,16,17, 52]. The presence of ridgelet in the curvelet transform may increase redundancy property. To solve these problems, several researchers [50, 51, 53] redesigned the concept of the curvelet transform to make it simple and easier to implement. The block diagram of the curvelet transform is shown in Fig. 23.

Fig. 23

Block diagram of first generation curvelet transform

Appendix B: Supplementary data

The performance of the different denoising methods for Gaussian and mixed Gaussian-Poisson noisy clinical PET images are shown in Figs. 24 and 25, respectively. The numerical results in Tables 17 and 18 show that the proposed method achieves better results than other denoising methods with respect to MSE, SNR, PSNR, and IQI. The detailed comparative performance analysis of different methods is described in 3.3. The average computation time of different denoising methods for two noisy PET images (Figs. 24 and 25) is shown in Table 19. The detailed analysis of the average computation time is discussed in Section 3.7.

Fig. 24

Performance comparison of various denoising techniques on Gaussian noisy (σ = 64) brain PET with noisy PSNR = 33 dB with 3×3 neighborhood window. a Original image. b Noisy image. c VisuShrink [19]. d BayesShrink [20]. e NeighShrink [23]. f ModineighShrink [24]. g Om et al. [25]. h Pogam et al. [7]. i Starck et al. [17]. j Proposed wavelet-curvelet

Fig. 25

Performance comparison of various denoising techniques on mixed Gaussian(σ = 80)-Poisson noisy brain PET image with noisy PSNR = 26 dB with 3×3 neighborhood window. a Original image. b noisy image. c VisuShrink [19]. d BayesShrink [20]. e NeighShrink [23]. f ModineighShrink [24]. g Om et al. [25]. h Pogam et al. [7]. i Starck et al. [17]. j Proposed wavelet-curvelet

Table 17 Metrics obtain for Gaussian noisy (σ = 64) brain PET (Fig. 24) image (noisy PSNR = 33 dB) and the denoised images using different denoising techniques with 3×3 neighborhood window

Table 18 Metrics obtain for the mixed Gaussian(σ = 80)-Poisson noisy brain PET (Fig. 25) image (noisy PSNR = 26 dB) and the denoised images using different denoising techniques with 3×3 neighborhood window

Table 19 Performance analysis between the proposed and existing methods with respect to average time requirement in seconds per PET image. Here two different PET images (Figs. 24 and 25) are used for measuring the time requirement