A modified OSEM algorithm for PET reconstruction using wavelet processing

Summary

Ordered subset expectation–maximization (OSEM) method in positron emission tomography (PET) has been very popular recently. It is an iterative algorithm and provides images with superior noise characteristics compared to conventional filtered backprojection (FBP) algorithms. Due to the lack of smoothness in images in OSEM iterations, however, some type of inter-smoothing is required. For this purpose, the smoothing based on the convolution with the Gaussian kernel has been used in clinical PET practices. In this paper, we incorporated a robust wavelet de-noising method into OSEM iterations as an inter-smoothing tool. The proposed wavelet method is based on a hybrid use of the standard wavelet shrinkage and the robust wavelet shrinkage to have edge preserving and robust de-noising simultaneously. The performances of the proposed method were compared with those of the smoothing methods based on the convolution with Gaussian kernel using software phantoms, physical phantoms, and human PET studies. The results demonstrated that the proposed wavelet method provided better spatial resolution characteristic than the smoothing methods based on the Gaussian convolution, while having comparable performance in noise removal.

Introduction

Positron emission tomography (PET) provides in vivo images of biological processes and has been increasingly utilized in clinical and investigative applications. PET images are obtained by reconstructing a function f, which represents the distribution of radioactivity in a body cross section, from its indirect data y, where y is the measured Radon transform $Rf$ (for definition, see Section 2). The difficulty in inverting the Radon transform $R$is that it is a smoothing transform, i.e., applying $R$ to an image f results in $Rf$ that has, roughly speaking, one-half derivatives more smoothness than the original image in two-dimension. Thus, $R^{-1}$, the inverse Radon transform, has properties like one-half of a derivative, i.e., it is unbounded. The difference between the measured data y and the true Radon transform $Rf$, makes the problem more difficult.

PET images can be reconstructed using either direct or iterative methods. All direct reconstruction methods can be cast in terms of the Fourier Projection Theorem (see Theorem 2.1 in [1]) of the Radon transform. Among them, the filtered back-projection (FBP) method is most commonly used. It provides fast reconstruction and reasonably good accuracy for the cases when noise level is low and the sufficient number of projections are available. The FBP method uses a low-pass filter to reduce the noise effect, which is usually dominant at high frequencies. Thus, the noise reduction obtained by FBP comes at the expense of degraded image resolution, since high frequency features in the image are also smoothed. Furthermore, the numerical implementation (FFT: Fast Fourier Transform) of the Fourier transform with continuous variables often causes aliasing artifacts in reconstructed images.

The major advantage of iterative methods over direct ones is the option to allow statistical noise models to be included as well as incorporation of prior knowledge. The maximum likelihood (ML) approach using the expectation–maximization (EM) algorithm [2] was introduced to deal with noise due to the low count rate and number of physical factors causing singularities and systemic biases frequently observed in the images reconstructed by FBP. The ML-EM reconstruction approach was further improved in terms of computation time by developing block iterative algorithms [3], [4].

PET images reconstructed by OSEM, however, still demonstrate poor noise characteristics. To overcome this problem, several EM-like iterative methods based on penalized least squares [5], [6] and Bayesian approaches [7], [8], [9], [10], [11] have been introduced. Even though penalized least squares and Bayesian methods are derived from different bases, both of them, incorporating with OSEM algorithm, give an additional inter-smoothing effect. In practice, the convolution with Gaussian kernel (this will be called as Gaussian smoothing (GS) from now on) is often applied to projections or images generated from iterations. Therefore, just like in FBP, the noise reduction obtained by inter-smoothing in OSEM with GS comes at the expense of degraded image resolution.

There has been great interest in the use of wavelet bases in signal/image processing applications. One of the key features of wavelet bases is the so-called energy concentration, which means that most of the energy in the image is concentrated on a few wavelet coefficients based on the smoothness of the image, while noises tend to spread to all wavelet coefficients. Using this property, several researchers apply wavelet shrinkage (WS) methods ([12], [13]), which shrink noisy wavelet coefficients towards zero, to noise removal problems, [14], [15], [16].

Wavelet methods are also used in direct image reconstructions [17], [18], [19], [20]. These wavelet-related methods provide better image reconstruction over FBP. As being a direct method, however, those wavelet methods still use the Fourier transform in computing either the wavelet coefficients or the backprojection part. Thus, those wavelet-based direct image reconstruction methods cannot be exempted from the aliasing effects caused by difference between the continuous Fourier transform and its discrete version, FFT. Obviously, those aliasing effects are more severe when there is an insufficient number of projections.

Due to an intrinsic nature of iterative method, it is rather difficult to give a simple statistical noise model on images generated in OSEM iterations. In fact, it has been shown that the noise in OSEM iterations can be modeled as log-normal distribution and the correlation between pixels is significant [21], [22], [23]. Thus, the direct application of WS, which is known to be good at removing the white noise (noise term in each pixel identically and independently follows the normal distribution with mean 0), to images in OSEM iterations is not a good approach: the log-normal distribution has much slowly decaying tail as compared with the normal distribution, and hence it can produce more frequently large noise terms. Obviously, those noise terms are not easily removed by WS with small shrinkage parameter. On the other hand, WS with large shrinkage parameter gives over-smoothed images.

In this paper, to remove the noise in images from OSEM iterations effectively, we suggest to use the standard WS and the robust WS (RWS) in a mixed manner. Roughly speaking, WS removes effectively noise terms with small intensities in images from OSEM iterations, whether they are independently distributed or correlated, while RWS does well noise terms with large intensities. By its definition, RWS potentially treats point source-like structure as noise and removes it. Therefore, RWS is not to be used too often and in the final iteration step. We suggested to use WS and RWS selectively in OSEM iterations, but fewer times for RWS than WS. For details, see Section 4. From now on, this hybrid method will be called as HWS.

We show through reconstruction experiments with software phantom data and real PET sinograms that the proposed method, HWS, provides better image reconstruction than GS. To support our claim, we compared the root mean square errors, the spatial resolution, and noise removal performance of the proposed wavelet method with those of Gaussian smoothing. In experiments with software phantom data, we generated the Poisson distribution to simulate emission related errors.

Other wavelet-related EM-like iterative reconstruction methods can be found in [24], [25], [26] , etc. In [24], [25] , authors used the EM algorithm on the multigrid resolution given by wavelets. By doing so, they reduced the required number of iterations for convergence and thus prevented noise contamination due to numerical errors from accumulating. As compared with their approach, our method used 2 or 3 as the fixed main iteration number. Those iteration numbers are commonly used in clinical PET practices for human studies. In this setting, our main concern is, obviously, not in reducing the iteration number required for convergence, but how to effectively suppress noise transferred from observed projections in each OSEM iteration.

Mair et. al. [26] applied Donoho and Johnstone’s SureShrink wavelet de-noising [13] to the observed projections and then used the standard EM algorithm, where the Anscombe transform [27] was used to stabilize the variance of Poisson projection data. They also discussed the method in which the wavelet de-noising step is inserted in each EM iteration. As compared with this, we used the wavelet de-noising method only to images, not to projections. By doing so, we can have an adaptive de-noising scheme, which is not easily obtainable by any operations on projections. We also used RWS in a very few OSEM iterations, which gives big contribution in removing noise.

This paper is organized as follows. In Section 2 , we briefly review the OSEM algorithm. Section 3 reviews the basic wavelet theory. In Section 4 , we present a wavelet-based OSEM iterations for PET image reconstruction. In Section 5 , we compare the performance of the proposed method with that of Gaussian smoothing in several criteria. The discussion and conclusion are given in Section 6.