A modified OSEM algorithm for PET reconstruction using wavelet processing

Computer Methods and Programs in Biomedicine (2005) 80, 236—245

A modified OSEM algorithm for PET reconstruction using wavelet processing Nam-Yong Leea Yong Choib,∗ a

School of Computer Aided Science, Inje University, Korea Department of Nuclear Medicine, Samsung Medical Center, Sungkyunkwan University School of Medicine, 50 Ilwon-Dong, Kangnam-Gu, Seoul 135-710, Korea

b

Received 22 April 2004 ; received in revised form 26 August 2005; accepted 23 September 2005

KEYWORDS PET reconstruction; Ordered subsets; Wavelet shrinkage

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. © 2005 Elsevier Ireland Ltd. All rights reserved.

1. 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 ∗

Corresponding author. E-mail address: [email protected] (Y. Choi).

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

0169-2607/$ — see front matter. © 2005 Elsevier Ireland Ltd. All rights reserved. doi:10.1016/j.cmpb.2005.09.004

A modified OSEM algorithm for PET reconstruction using wavelet processing 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 backprojection (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—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

237

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—16]. Wavelet methods are also used in direct image reconstructions [17—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—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,

238

N.-Y. Lee, Y. Choi

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—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 denoising 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.

2. OSEM algorithm The purpose of the image reconstruction is to recover an image f from its projections Rf, where R is the Radon transform defined by
Rf(
, u) =

∞

−∞




∞

−∞

f(x, y)ı(u − x cos


− y sin
) dx dy.

(1)

Here, ı is the Dirac delta function. In practice, we have only discrete and noisy projections (yd )d instead of the true projections Rf(
, u) with continuous variables. The EM algorithm begins with the independent Poisson distribution assumption on projections (yd )d , i.e., Prob(yd = k) = e−d

kd , k!

k = 0, 1, . . . ,

where
d = Rf(
, u) d
du, Jd

where Jd is the unit region in the domain of projections. The boldface in (yd )d indicates that (yd )d are regarded as random variables. The ML-EM approach is based on the maximization of the probability of (yd )d as a function of (fb )b given observed samples (yd )d under positivity constraints fb ≥ 0, where fb are discretizedsamples for the original image f, for instance, f = b fb Ib , where Ib is the unit region in the domain of images and I is the indicator function, i.e., I (x) = 1 if x ∈ I, 0 otherwise. In the resulting EM iteration, the new image f n+1 is updated from the current image f n by    yd n+1 n fb =fb Pb,d n Pb,c , nd = Pb,d fbn , d c d

b

(2) where the transition matrix P is computed by Pb,d = J RIb (
, u) d
du. d The described EM algorithm had two major drawbacks: huge computing requirement and noisy images generated in EM iterations. The first problem has almost been overcome by a large reduction in the required iteration number for acceptable convergence (see, e.g., [3,4]) and rapid improvement in computing facilities in the last decade. Thus, even though OSEM are much slower than conventional Fourier transform-based ones such as FBP, in some extent, the computing time required for OSEM has become no longer a big problem. On the other hand, the second problem still needs to be improved. Noisy images are mainly due to the approximation error, the noise transferred from observed projections, and the numerical computation error accumulated from iterations, etc. To reduce noises in reconstructed images we must use some kind of smoothing to noisy projections or reconstructed images to suppress noise magnification and finish the iteration as fast as possible to prevent the numerical computation error from

A modified OSEM algorithm for PET reconstruction using wavelet processing accumulating. For these purposes, several EM-like iterative methods based on penalized least squares and Bayesian approaches have been introduced. Even though theoretical bases of penalized least squares and Bayesian methods are different, their resulting iterative algorithms have steps that give some inter-smoothing effect. Based on this result and to get the smoothing effect directly, as mentioned before, GS is commonly used in practice. We write the following modified OSEM iterations, which apply inter GS on projections in iterations and post GS after finishing iterations. Here, N and M is the main iteration and sub-iteration numbers, respectively. OSEM with GS

where Z is the set of all integers, is called an orthonormal wavelet basis for L2 (R), the set of all square integrable functions defined on the whole real line R, if for any f ∈ L2 (R), f=



f,

k,j  k,j .

k∈Z j∈Z

Associated with is a scaling function , from which one can generate functions k,j (x) = 2k/2 (2k x − j). With k,j and k,j , we have    f, m,n m,n = f, k,j  k,j n∈Z

k0 ≤k
+

fb0,0

(1) Let = 1 for all b. (2) For n = 0, 1, . . . , N − 1 (a) For subsets m , m = 0, 1, . . . , M − 1, (i) Project: n,m = Pb,d fbn,m , for d ∈ m d

˜ f = GS2 (f

N,0

f, k0 ,l k0 ,l

for any integers k0 , m, k0 < m and f ∈ L2 (R). For two dimensional orthonormal wavelets on L2 (R2 ), by using tensor products of  and , we make four functions as

b

(b) f n+1,0 = f n,M (3) Post GS

 l∈Z

(ii) Inter GS:  ˜ n,m = GS1 (n,m ) (iii) Backproject and update:    yd fbn,m+1 = fbn,m Pb,d n,m Pb,c  ˜d c∈ d∈m

239

m

)

Here, GS1 and GS2 mean one and two dimensional Gaussian smoothing, respectively. The smoothing degree can be measured by full width half maximum (FWHM) of the Gaussian kernel. To simplify our presentation, we consider the case with the same FWHM in inter GS and post GS. The described GS with proper FWHM provides better reconstructions than not only OSEM without smoothing but also FBP. However, they cannot avoid possible degraded image resolution due to the intrinsic nature inherited in the convolution-based smoothing process.

3. Wavelet processing In this section, we shall briefly review the basic wavelet theory and the wavelet-related smoothing methods. More details can be found in [13,14,28], etc. We begin with the one dimensional case. The term ‘‘wavelets’’ refers to a collection of basis elements k,j (x) = 2k/2 (2k x − j), which are generated by translations and dilations from a single function . The collection of function { k,j }k∈Z,j∈Z ,

(x, y) = (x) (y) 1 (x, y)

= (x) (y)

2 (x, y)

=

(x)(y)

3 (x, y)

=

(x) (y).

We define k,j = 2k (2k · −j) for k ∈ Z and j ∈ Z2 , i . Then for any f ∈ L2 (R2 ) we and similarly for k,j have     i i f, m,n m,n = f, k,j  k,j k0 ≤k
n∈Z2

+



f, k0 ,l k0 ,l

l∈Z2

for any integers k0 , m, k0 < m. From now on, we shall use the following notations Cnk (f) = f, k,n  k

and

Dnk,i (f) = f,

i k,n .

The scaling Cn0 (f) and wavelet coefficients Dnk,i (f) in level k = k0 , k0 + 1, . . . , m − 1 can be obtainable from scaling coefficients Cnm (f) in level m by successively using the so-called fast wavelet transform (FWT). For details, see, e.g., [28]. When we are concerned with a digital image with pixel values (fn )n , n = (n1 , n2 ), 0 ≤ n1 , n2 < 2m − 1 (to simplify our presentation, we consider only square shape images), we assume that there exist a square integrable function f defined on  = [0, 1] × [0, 1] in R2 and fn = Cnm (f). From now on, we do not distinguish the image with the corresponding square integrable function with the

240

N.-Y. Lee, Y. Choi

described sense. In this case, we do not consider (i) all shifts j ∈ Z2 , but only those shifts for which k,j intersects  nontrivially. With these assumption we can compute wavelet coefficients Djk,i (f) of the image f by using FWT and recover the image f from its wavelet coefficients by using inverse FWT. The FWT has some advantages over FFT in approximation rate, computation speed, timefrequency localization, etc. Fast approximation means that most of the energy in the image f tends to be concentrated on a few wavelet coefficients Dnk,i (f). A similar phenomenon holds for FFT too, but its approximation rate is much slower than that of FWT. On the other hand, the white noise tends to spread to all wavelet coefficients. Using this property, several researchers [12—15] applied wavelet de-noising methods to the white noise removal problems. Suppose a noisy image y = f + z, where f is a true image, is contaminated by a white noise z. The basic structure of most of the wavelet de-noising methods follows the next three steps: • Apply FWT to noisy image y to get its wavelet coefficients Dnk,i (y). • Use soft or hard thresholding to Dnk,i (y). • Apply inverse FWT to the resulting coefficients from the previous step. Here, soft thresholding S (x) and hard thresholding H (x) are defined by    x −  if x > , S (x) = 0 if |x| ≤ ,   x +  if x < −, and

H (x) =

x if |x| > , 0 if |x| ≤ .

This wavelet de-noising method, which is called wavelet shrinkage (WS), provides asymptotically near-optimal results for white noise removal if a well-chosen thresholding value  is employed, and gives better noise reduction than the classical methods such as GS do. Moreover, WS employs a spatially adaptive noise removal scheme [12,14]. WS removes small wavelet coefficients by regarding them as noise and leaves large ones as signal. This is the main advantage of WS over GS; GS can be viewed as a low-pass windowing in the frequency domain. Thus GS considers any low-frequency structure to be signal, and any high-frequency structure to be noise, no matter how large the high-frequency structure might be. This is not acceptable when we

wish to reconstruct images that can be more meaningfully characterized by their discontinuities (e.g., strong edges and small extent), since such information lies in the high-frequency domain.

4. Proposed method As mentioned in Section 1, our main concern in this paper is how to suppress noises transferred from observed projections so that they are not magnified in each OSEM iteration with a fixed small main iteration number such as 2 or 3. The measurement error in observed projections comes mainly from the intrinsic nature of radioactive decay, scattering, and the error related to attenuation corrections, etc. We use the Poisson distribution for the error related to the radioactive decay the white noise model for errors related to scattering, attenuation correction, and other sources. The noise in images in OSEM iterations, however, neither follows the Poisson distribution nor white noise model. There have been some research on the noise properties of OSEM images [21—23]. In these theoretical and related Monte Carlo studies, it has been shown that noise in images from OSEM iterations can be model as log-normal distribution and the correlation between pixels is significantly high. Since the log-normal distribution has much slowly decaying tail as compared with the normal distribution, there is more possibility of observing large noise terms in images from OSEM iterations than corrupted images by additive white noise. In this paper, we call those large noise terms as outliers, in a sense that they are not small enough to be considered as samples from normal distribution. Therefore, in this paper, outliers not only include point source like noise terms but also noise terms with relatively large intensities. Most of the smoothing algorithms do not remove well white noise and outliers simultaneously. The low-pass filtering such as GS does not remove outliers well; the energy on an outlier is spread over a wide frequency range so that the usual low-pass filter cannot remove outliers effectively. On the other hand, the median filter is well known to be good at removing outliers. A considerably long mask is, however, required for the median filter to suppress white noise effectively, and this might result in an over-smoothed image. The standard WS does not work well on outlier removal, either; most of the energy of an outlier is concentrated on a few wavelet coefficients that are too large to be removed by the standard WS without causing some degradation.

A modified OSEM algorithm for PET reconstruction using wavelet processing To overcome the described weakness of low-pass filtering, median filtering, and standard WS, we suggest to use a hybrid method of the standard WS and a robust WS (RWS) for smoothing of images in OSEM iterations. The RWS was proposed by Bruce et al. [29] and successfully used in estimations of frequency response function [30] and ultrasound pulses [16]. To explain RWS, let us consider the noise removal problem y = f + z + O, where the true image f is contaminated by white noise z and outliers O. Let Cnm = Cnm (y) be the pixel values of the observed noisy image y. Starting from Cnm , for each level k (k ≤ m) Cnk is decomposed into three components: ˜nk ), • Outlier coefficients Onk given by Onk = Sε (Cnk − C ˜nk is where Sε is the soft shrinkage operator and C k a median filtered coefficients of Cn • New scaling coefficients Cnk−1 obtained by applying the low-pass FWT to the outlier-free coefficients Cnk − Onk • New wavelet coefficients Dnk−1 obtained by applying the high-pass FWT to the outlier-free coefficients Cnk − Onk The outlier selection step is just used in the first few wavelet decomposition levels. The remaining wavelet decomposition is performed in the standard manner. The proposed method, HWS, uses a hybrid of WS and RWS. In summary, we write the proposed method as follows: OSEM with HWS (1) Let fb0,0 = 1 for all b and ˜ f 0,0 = f 0,0 . (2) For n = 0, 1, . . . , N − 1 (a) For subsets m , m = 0, 1, . . . , M − 1, (i) Project: n,m = Pb,d ˜ fbn,m , for d ∈ m d b

(ii) Backproject and update:    yd fbn,m+1 =˜ fbn,m Pb,d n,m Pb,c d c∈ d∈m

m

(iii) Wavelet processing:  n,m+1 ) if m is odd,   WS(f ˜ f n,m+1 = RWS(f n,m+1 ) if m=M/2−2,   n,m+1 f otherwise. f n,M (b) ˜ f n+1,0 = ˜ (3) ˜ f =˜ f N,0 The proposed method does not have a post processing, thus, the reconstructed image ˜ f at the final stage is ˜ f N,0 . Here, WS(f) and RWS(f) indicate the resulting images after applying WS and RWS to f, respectively. As mentioned in Section 1, RWS can accidentally removes point source-like structures. Even for such

241

case, the remaining M/2 sub-iterations, which do not have RWS, can sufficiently restores those point source-like structures from observed projections y, since y are not changed in iterations and contain all information about the true image.

5. Experiments We conducted OSEM-based PET image reconstruction experiments with HWS and GS. Our main conclusion is that the proposed wavelet method, HWS, outperforms convolution-based smoothing methods, GS. For simulation studies, we used monkey autoradiography software phantom of size 128 × 128. We generate five noisy projection data sets of size 336 × 281 by the Poisson distribution with total counts of data to be 600,000, 800,000, 1,000,000, 1,200,000, and 1,400,000. We also used real projection data sets obtained from a GE Advance tomograph (Milwaukee, WI) in two-dimensional (septaextended) mode [31,32]. The dimension of projection data from this tomograph is also 336 × 281. All experiments were performed with M = 16 as the number of ordered subsets. In ordered subsets, angle indices were set to be equally distributed to each subset. The subset order was chosen such that the angular distance between successive subsets was maximized, i.e., 0 = {
0 ,
1 , . . . ,
20 }, 1 = {
168 ,
169 , . . . ,
188 }, 2 = {
21 ,
22 , . . . ,
41 }, . . ., 15 = {
315 ,
316 , . . . ,
335 }, where
n = n/336. The performance of each smoothing method was evaluated by estimating the following figure of merits under various image considerations. (1) Image root-mean-square (RMS) Image RMS was employed to find the global differences between a reference software phantom and the images reconstructed using OSEM iterations with HWS and GS. The RMS between the original image f = (fb ) and the reconstructed image ˜ f = (˜ fb ) is defined by  |fb − ˜ f b |2 . RMS = b  2 b fb (2) Spatial resolution A point source located at the center of the field of view of the PET scanner was imaged. A profile of the source fitted with a standard Gaussian function and the spatial resolution was taken as the FWHM of the fitted profile curve. (3) Coefficient of variation (CV) CV was measured using a uniform cylinder (20 cm in diameter) containing uniform activity

242

N.-Y. Lee, Y. Choi of 3 mCi of F-18. It is defined as

R CV = , R

22k is the total number of wavelet coefficients Dnk,i for fixed k and i. To estimate , we used

where R and R denotes the standard deviation and mean of the region of interest (ROI), respectively.

ˆ =

As a choice of wavelets, we used piecewise quintic biorthogonal wavelets that are denoted by ˜ , 6,8 ˜ , 6,8 , and 6,8 in [28]. The wavelet de6,8  noising method (WS or RWS) based on the standard wavelet basis often exhibits certain artifacts, due to the lack of translation invariance. For details, see, e.g., [33]. To avoid those artifacts, whenever we applied wavelet de-noising methods, we used the so-called redundant wavelet basis [33] instead of the standard one. In RWS, the size of the median filter is determined to just the five nearest neighboring pixels (up, down, left, right, center), to prevent outliers from leaking into coarser levels. The depth, at which median filters are operated, is set to be 1, and the wavelet decomposition depth is set to be 2. Therefore, RWS behaves like the standard WS for white noise removal in the first two levels, but in contrast with the standard WS, it effectively reduces outliers by an additional median filtering in the first wavelet decomposition level. The proposed method requires shrinkage parameters in WS or RWS steps. We suggest to use the soft wavelet shrinkage with  2 log 22k . (3) This is motivated by the so-called universal parameter suggested by Donoho and Johnstone [12]. Here,

√

2 mediann,i (|Dnm−1,i |)

as the estimated standard deviation. Recall that this standard deviation to be estimated is not the one of the white noise in projection data set, but in the images generated from OSEM iterations. The parameter ε in RWS is designed to distinguish outliers from white noise. We suggest to use ε = 1.96ˆ

that is sufficiently large to have a 95% confidence and sufficiently small to locate most of outliers. For GS, we conducted simulations with various FWHMs from 1.0 to 4.0 with the increment 1.0 both in pixel unit. Table 1 shows RMS values of HWS and GS. It also shows results by WS (HWS without the robust wavelet shrinkage) and RWS (HWS without the standard wavelet shrinkage). As we can see in Table 1, HWS provided slightly smaller RMS values than WS, RWS, and GS for five different photon counts. All reconstructed images in Figs. 1–3 and the results in Table 2 were obtained with N = 3 as the main iteration number. Fig. 1 shows software phantoms, monkey brain (a), and reconstructed ones by OSEM iterations with HWS (b), WS (c), RWS (d), GS with FWHM = 1.0 (e), FWHM = 2.0 (f), FWHM = 3.0 (g), and FWHM = 4.0 (h) for the noisy projection data set generated by the Poisson distribution with 1,000,000 photon counts. For monkey brain software phantom, images reconstructed by HWS demonstrated the activity distribution patterns as those by WS, RWS,

Fig. 1 (a) Original monkey brain software phantom image. Reconstructed image (b) by HWS, (c) by WS, (d) by RWS, (e) by GS with FWHM = 1.0, (f) by GS with FWHM = 2.0, (g) by GS with FWHM = 3.0, and (h) by GS with FWHM = 4.0. All reconstructed images are computed from the noisy projection data set generated by the Poisson distribution with 1,000,000 photon counts.

A modified OSEM algorithm for PET reconstruction using wavelet processing

Fig. 2 (d).

243

Hoffman brain images reconstructed by HWS (a) and GS with FWHM = 1.0 (b), FWHM=2.0 (c), and FWHM=3.0

Fig. 3 Human FDG-PET brain images reconstructed by HWS (a) and GS with FWHM1.0 (b), FWHM = 2.0 (c), and FWHM = 3.0 (d). Table 1

RMS values (after multiplying 10,000) and corresponding parameters are shown N=3

N=2

Counts

0.6

0.8

1.0

1.2

1.4

0.6

0.8

1.0

1.2

1.4

HWS WS RWS

472 611 552

418 542 449

400 442 417

368 419 378

356 412 356

550 665 591

492 602 507

477 512 474

441 487 441

431 450 431

GS: FWHM = 1.0 FWHM = 2.0 FWHM = 3.0 FWHM = 4.0

820 478 595 772

654 421 583 780

571 414 568 788

517 403 556 778

452 400 543 771

708 553 674 842

591 504 667 856

535 496 653 863

510 493 646 855

468 492 634 849

The total photon counts for Poisson distribution are written in million units.

and GS, while preserving sharp edge structures and removing more noises. From Table 1 and Fig. 1, it is obvious that the hybrid wavelet method, HWS, is better than WS or RWS. From now on, we omit the results by WS or RWS. Fig. 2 shows Hoffman brain phantom images reconstructed by HWS (a) and GS (b, c, and d). For GS, we showed images reconstructed with FWHM = 1.0 (b), 2.0 (c), and 3.0 (d). Through this experiment, we believe that HWS removes as much noises as GS does, while preserving more of the fine structure of the image. Table 2 shows spatial resolution and noise removal test results by HWS and GS. Unlike GS, where the smoothing effect is totally controlled by the FWHM of the Gaussian kernel used in the convolution-based smoothing, the smoothing effect of HWS is dependent on the given projection data, since the shrinkage parameter in HWS is deter-

mined through noise intensity estimation, which is certainly dependent on the given data itself. To have a fair comparison, we used the same wavelet shrinkage parameters used in Hoffman image reconstructed by HWS in Fig. 2 to spatial resolution and noise removal tests. Through these results, we can conclude that HWS provides sharper spatial resolutions and better noise removal simultaneously as

Table 2 Spatial resolution and noise removal as well as corresponding parameters are shown

FGS FSR CV

HWS

GS

N/A 2.4 0.023

1.0 2.7 0.035

2.0 2.5 0.032

3.0 3.2 0.027

4.0 4.2 0.024

For all experiments, N = 3 was used as the main iteration number. Here, we used FGS and FSR for FWHM of the Gaussian kernel in GS and that in spatial resolution test, respectively.

244 compared with GS for a wide range of FWHMs of Gaussian kernels. Fig. 3 compares human FDG-PET brain images reconstructed by HWS and GS. It demonstrates similar results to those obtained by simulation with reference data and physical phantom studies.

N.-Y. Lee, Y. Choi tional Nuclear Technology Program and by a grant of the Korea Health 21 R&D Project (02-PJ3-PG6EV06-0002), Ministry of Health & Welfare, Republic of Korea.

References 6. Discussion and conclusion The application of the wavelet framework to the problem of PET image reconstruction by OSEM iterations was demonstrated conceptually and experimentally. The proposed method, HWS, is based on a hybrid algorithm of the standard WS and the robust WS. Several criteria are used to compare the performance of HWS with those of GS. Through these simulation studies, we conclude that the HWS smoothing method gives stable image reconstruction in incorporation with OSEM iterations and provides better resolution than WS, RWS, and GS, while preserving a comparable noise removal. The RWS de-noising method is adapted as an alternative tool to overcome the limitation of the standard WS, the convolution based smoothing, and median filtering in removing the white noise and outliers simultaneously. This method uses a wavelet-related median filtering to suppress outliers while dramatically reducing the white noise with the standard WS. The overuse of RWS in incorporation with OSEM iterations, however, may results in blocky artifacts in reconstructed images. In most of the cases, those artifacts do not harm the image RMS. In visual quality, however, they give some coarseness in reconstructed images. Because of this reason, even though it is very subjective, we suggest to use RWS just once per one main iteration and not to be used near at final iteration step. In HWS, we used the shrinkage parameter motivated by Donoho and Johnstone’s suggestion as the universal one. We believe that there is a close relationship between optimal shrinkage parameters and the smoothness of the original image: for smooth images, larger shrinkage parameters are preferred, while smaller ones provide better results for unsmooth images. Finding an optimal decision rule in selecting the amount of wavelet shrinkage based on the smoothness of the image to be reconstructed will be our future research program.

Acknowledgments This study is supported by Korea Institute of Science & Technology Evaluation and Planning (KISTEP) and Ministry of Science & Technology through their Na-

[1] F. Natterer, F. W¨ ubbeling, Mathematical methods in image reconstruction, SIAM, Philadelphia, 2001. [2] L.A. Shepp, Y. Vardi, Maximum likelihood reconstruction for emission tomography, IEEE Trans. Med. Imaging 1 (2) (1982) 113—122. [3] H.M. Hudson, R.S. Larkin, Accelerated image reconstruction using ordered subsets of projection data, IEEE Trans. Med. Imaging 13 (4) (1994) 601—609. [4] J. Browne, A.R. De Pierro, A row-action alternative to the EM algorithm for maximizing likelihoods in emission tomography, IEEE Trans. Med. Imaging 15 (5) (1996) 687—699. [5] L. Kaufman, Maximum likelihood, least squares, and penalized least squares for PET, IEEE Trans. Med. Imaging 13 (2) (1993) 200—214. [6] J.A. Fessler, Penalized weighted least-squares image reconstruction for positron emission tomography, IEEE Trans. Med. Imaging 13 (2) (1994) 290—300. [7] T. Hebert, R. Leahy, A generalized EM algorithm for 3-D Bayesian reconstruction from Poisson data using Gibbs priors, IEEE Trans. Med. Imaging 8 (2) (1989) 194—202. [8] K. Lange, Convergence of EM image reconstruction algorithms with Gibbs smoothing, IEEE Trans. Med. Imaging 9 (4) (1990) 439—446. [9] P.J. Green, Bayesian reconstruction from emission tomography data using a modified EM algorithm, IEEE Trans. Med. Imaging 9 (1) (1990) 84—93. [10] A.R. De Pierro, A modified expectation maximization algorithm for penalized likelihood estimation in emission tomography, IEEE Trans. Med. Imaging 14 (1) (1995) 132— 137. [11] A.R. De Pierro, M.E.B. Yamagishi, Fast EM-like methods for maximum a posteriori estimates in emission tomography, IEEE Trans. Med. Imaging 20 (4) (2001) 280—288. [12] D.L. Donoho, I.M. Johnstone, Ideal spatial adaptation by wavelet shrinkage, Ann. Stat. 81 (1994) 425—455. [13] D.L. Donoho, I.M. Johnstone, Adapting to unknown smoothness via wavelet shrinkage, J. Am. Stat. Assoc. 90 (1995) 1200—1224. [14] A. Chambolle, R.A. DeVore, N.-Y. Lee, B.J. Lucier, Nonlinear wavelet image processing: variational problems, compression, and noise removal through wavelet shrinkage, IEEE Trans. Image Process. 7 (3) (1998) 319—335. [15] R.A. DeVore, B.J. Lucier, Fast wavelet techniques for nearoptimal image processing, IEEE Military Communication Conference, 1992, pp. 1129—1135. [16] O. Michailovich, D. Adam, Robust estimation of ultrasound pulses using outlier-resistant de-noising, IEEE Trans. Med. Imaging 22 (3) (2003) 368—381. [17] E.D. Kolaczyk, A wavelet shrinkage approach to tomographic image reconstruction, J. Am. Stat. Assoc. 91 (1996) 1079—1090. [18] M. Bhatia, W.C. Karl, A.S. Willsky, Tomographic reconstruction and estimation based on multiscale natural-pixel bases, IEEE Trans. Image Process. 6 (1997) 463—478. [19] N.-Y. Lee, B.J. Lucier, Wavelet methods for inverting the Radon transform with noisy data, IEEE Trans. Image Process. 10 (1) (2001) 79—94.

A modified OSEM algorithm for PET reconstruction using wavelet processing [20] Y. Choi, J.-Y. Koo, N.-Y. Lee, Image reconstruction using the wavelet transform for positron emission tomography, IEEE Trans. Med. Imaging 20 (11) (2001) 1188—1193. [21] D.W. Willson, B.M.W. Tsui, H.H. Barrett, Noise properties of the EM algorithm II: Monte Carlo simulations, Phys. Med. Biol. 39 (5) (1994) 847—871. [22] H.H. Barrett, D.W. Willson, B.M.W. Tsui, Noise properties of the EM algorithm I: Theory, Phys. Med. Biol. 39 (5) (1994) 833—846. [23] E.J. Soares, C.L. Byrne, S.J. Glick, Noise characterization of block-iterative reconstruction algorithms: I. Theory, IEEE Trans. Med. Imaging 19 (4) (2000) 261—270. [24] M.V. Ranganath, A.P. Dhawan, A. Mullani, A multigrid expectation maximization reconstruction algorithm for positron emission tomography, IEEE Trans. Med. Imaging 7 (4) (1988) 273—278. [25] A. Raheja, T.F. Doniere, A.P. Dhawan, Multiresolution expectation maximization reconstruction algorithm for positron emission tomography using wavelet processing, IEEE Trans. Nucl. Sci. 46 (3) (1999) 594—602. [26] B.A. Mair, R.B. Carroll, J.M. Anderson, Filter banks and the EM algorithm, IEEE Medical Imaging Conference, 3, 1996, pp. 1747—1751.

245

[27] F. Anscombe, The transformation of Poisson, binomial and negative binomial data, Biometrika 35 (1948) 246—254. [28] I. Daubechies, Ten Lectures on Wavelets, SIAM, Philadelphia, 1992. [29] A.G. Bruce, D.L. Donoho, H.Y. Gao, R.D. Martin, Denoising and robust non-linear wavelet analysis, Proceedings of SPIE Conference, 1994, pp. 325—336. [30] Y.Y. Kim, J.-C. Hong, N.-Y. Lee, Frequency response function estimation via a robust wavelet de-noising method, J. Sound Vibr. 244 (4) (2001) 635—649. [31] T.R. DeGrado, T.G. Turkington, J.J. Williams, C.W. Stearns, J.M. Hoffman, R.E. Coleman, Performance characteristics of a whole-body PET scanner, J. Nucl. Med. 35 (1994) 1398— 1406. [32] J.R. Lee, Y. Choi, Y.S. Choe, K.H. Lee, S.E. Kim, S.A. Shin, B.-T. Kim, Performance measurements of positron emission tomography: an investigation using General Electric AdvanceTM Korean, J. Nucl. Med. 30 (1996) 548— 559. [33] R. Coifman, D.L. Donoho, Translation invariant de-noising, in: A. Antoniadis, G. Oppenheim (Eds.), Wavelet and Statistics Lecture Notes in Statistics, Springer-Verlag, 1995, pp. 125—150.