- Email: [email protected]

Contents lists available at ScienceDirect

Optik journal homepage: www.elsevier.de/ijleo

Application of regularized maximum likelihood algorithm in PET image reconstruction combined with nonlocal fuzzy anisotropic diffusion Quan Zhang a,b , Yi Liu a , Huazhong Shu b , Zhiguo Gui a,∗ a b

National Key Laboratory for Electronic Measurement Technology, North University of China, Taiyuan 030051, China Laboratory of Image Science and Technology, Southeast University, Nanjing 210096, China

a r t i c l e

i n f o

Article history: Received 30 August 2012 Accepted 19 January 2013

Keywords: Positron emission tomography Anisotropic diffusion Regularized maximum likelihood algorithm Fuzzy set Nonlocal means

a b s t r a c t Although the traditional regularized maximum likelihood (RML) algorithms can obtain the high-quality reconstructed image in positron emission tomography (PET), it still suffers from over-smoothing with the increase of the iteration number. In this paper, we propose a novel regularized maximum likelihood reconstruction algorithm for PET imaging, which uses the fuzzy nonlinear anisotropic diffusion as the regularization to reduce noise and preserve edges. In addition, we introduce the nonlocal means idea which can adequately exploit the global information of the image into the reconstruction process to further improve the quality of the result image. The proposed algorithm not only absorbs the advantages of fuzzy set theory in deal with uncertain problems, but also has the good ability of anisotropic diffusion, namely protecting edges perfectly and suppressing noise. Experimental results and quantitative analysis show that the novel method has the more excellent performance for positron emission tomography imaging. © 2013 Elsevier GmbH. All rights reserved.

1. Introduction Owning to the ability of providing the valuable information of the biochemical and biological activity inside a living subject by a non-invasive way, positron emission tomography (PET) opens up the exciting and rapidly evolving ﬁeld of molecular imaging [1]. However, although the molecular imaging by PET can be applied to neuroscience research, the quantitative accuracy of PET images is limited by the imperfect system response and the reconstruction methods which are used to reconstruct the images from projection measurements [1,2]. In fact, PET reconstruction is a complicated problem, and many algorithms have been developed in the past decades. PET imaging is an ill-posed problem which means that the solution may be nonexistent, not unique or unstable, and the noise level of PET projection data is always relatively signiﬁcant. People often employ the regularized technology to control the noise propagation so as to produce a reasonable reconstruction. The normal simple quadratic membrane (QM) prior [3] has the tendency to smooth both the high frequency edge and the low frequency region during the iterative process, while edge-preserving nonquadratic [4] prior such as median root prior (MRP) tends to produce blocky artifacts. Furthermore, the applications of edge-preserving Bayesian methods are often hindered by the needed annealing procedure

∗ Corresponding author. Tel.: +86 351 3557226; fax: +86 351 3557226. E-mail addresses: [email protected], [email protected] (Z. Gui). 0030-4026/$ – see front matter © 2013 Elsevier GmbH. All rights reserved. http://dx.doi.org/10.1016/j.ijleo.2013.01.092

and some complicated parameter estimation, and this will make the solution of posterior energy function be a more difﬁcult and complicated problem [5]. Since Zedeh introduced the fuzzy set theory in his literature [6], the fuzzy technology has already been applied in several ﬁelds of image processing such as ﬁltering [7], enhancement [8]. But there have been few researches combining with the fuzzy set in positron emission tomography until recently. The iterative algorithm using fuzzy potential function proposed by Modal et al. [9] demonstrates that the method is very effective to suppress the noise without damaging the features of the reconstructed image. The nonlinear anisotropic diffusion (AD) ﬁlter, ﬁrst proposed by Perona and Malik [10], is a nonlinear ﬁlter that purports to smooth a noisy image without blurring edges. It has been widely used for image denoising, image enhancement, image segmentation [11–13], and often obtains better quality than other methods. The basic idea of anisotropic diffusion is to choose a diffusion coefﬁcient adaptively in different diffusion regions so that ﬂat regions are smoothed out and edges are reserved [10]. Chen and his colleagues [14] suggested that the diffusion coefﬁcient can be constructed by the fuzzy membership function in order to get an accurate and ﬂexible control of the noise smoothing and edge preserving. Recently, Buades et al. [15] adopted a nonlocal means (NLM) approach for image denoising. It adjusts each pixel value with a weighted average of other pixels. Illuminated by the above-mentioned ideas, we propose a novel approach for the PET imaging, which can obtain the superior reconstructed image. This proposed method combines a fuzzy nonlinear anisotropic diffusion with the nonlocal means to reconstruct

4562

Q. Zhang et al. / Optik 124 (2013) 4561–4565

an adaptive diffusion coefﬁcient according to the characteristics of the image. The proposed method can remove the noise while preserving the detail information of the image, and the reconstructed image can be greatly improved. Simulated experiments show that the new approach can acquire a good balance between lowering noise and preserving edges. 2. The proposed fuzzy anisotropic diffusion ﬁltering Recently, image processing technology based on partial differential equation (PDE) has become a research hotspot [10,16,17]. To some extent, it overcomes the defects of traditional ﬁlters of easily leading to decline of the image spatial resolution and causing the ﬂattening of sharp edges [18]. And it is a family of adaptive smoothing technology, which can adaptively keep the sharpness of images. That is, the diffusion coefﬁcient can automatically increase in ﬂat regions to remove noise, while can adaptively decrease near the edge of image to preserve the edge information. In 1990, Perona and Malik [10] ﬁrst proposed the following anisotropic diffusion equation for an image f:

⎧ ⎨ ∂f = div(c (|∇ f |)∇ f ) P−M ∂t , ⎩

(1)

where the diffusion coefﬁcient is 1 1 + (|∇ f |/K)

or cP−M (|∇ f |) = exp

(2)

2

2 |∇ f | −

K

.

(3)

And the discrete form of P–M equation can be deﬁned as follow: f

h+1

= f + div(c(|∇ f |)∇ f ), h

(4)

where ∈ [0, 1] controls the strength of diffusion process. Usually, a small value of is used to avoid destabilizing the diffusion process. The function c(•) is a nonnegative monotonically decreasing function with c(0) = 1.0and c(b) = 0 when b → ∞. Black et al. [19] thought that the P–M function can be replaced by other new functions such as Tukey’s biweight function. Furthermore, the diffusion process can be considered as a fuzzy classiﬁcation in the view of fuzzy mathematics. That is, the more pixels belonging to the ﬂat region, the stronger the diffusion effects. Therefore, the diffusion coefﬁcient can be controlled by a membership function which gives the degree of belongingness in the region to be smoothed. Accordingly, the membership function, which maps the value of pixel into a numerical value [0,1] (a fuzzy space), can represent the diffusion coefﬁcient. In this paper, we adopt a membership function described as follow: (|∇ f |) = ω exp

2 |∇ f | −

K

,

(5)

where the weight ω is a positive value that denotes the interaction degree between the center pixel and the neighborhood pixel. And it is usually considered to be inversely proportional to the distance between the center pixel and the neighborhood pixels. Here, we propose a novel membership function by introducing the nonlocal means (NLM) idea which is ﬁrst proposed by Buades et al. [15] for image denoising when calculating the weight ω. For the pixel b and the pixel j, the corresponding weight ωbj can by denoted by the following formulas: ωbj =

exp(−||f (nb ) − f (nj )||2E /h2 ) dbj

f (nb ) = {fl : l ∈ nb },

,

(6)

(7)

f (nj ) = {fl : l ∈ nj }

(8)

||f (nb ) − f (nj )||2E =

2

(fl ∈ nb − fl ∈ nj ) ,

(9)

l

dbj =

f (x, y, 0) = f0 (x, y)

cP−M (|∇ f |) =

Fig. 1. The 3D plot of the proposed nonlocal fuzzy anisotropic diffusion coefﬁcient.

1

2

(bx − jx ) + (by − jy )

2

.

(10)

where f(nb ) and f(nj ) represent the pixel values in the two comparing neighborhood nb centered on the pixel b and nj centered on the pixel j, respectively. The parameter h is a factor which controls the decay of the exponential function. (bx , by ) and (jx , jy ) are the 2D positions of the pixel b and the pixel j, respectively. The weight ωbj reﬂects the degree of similarity or connectivity between the pixel b and the pixel j in the neighborhood. Based on the above description, the proposed nonlocal fuzzy anisotropic diffusion coefﬁcient (| ∇ f|, ωbj ) and the ﬂux function (∇ f, ωbj ) can be formalized as follows: (|∇ f |, ωbj ) = ωbj exp

2 |∇ f | −

K

,

(∇ f, ωbj ) = (|∇ f |, ωbj ) · ∇ f = ωbj exp

(11)

2 |∇ f | −

K

· ∇f

(12)

In (11) and (12), for a noise image, how to determine the gradient threshold K is very critical for the performance of the ﬁlter. Perona et al. [10] proposed using Canny’s noise estimator to get K. Torkamani-Azar et al. [20] suggested using the mean of the absolute gradient as K. In emission tomography, the gradient generated by noise is comparable to or even larger than the edge gradients. Under such a condition, choosing a proper K is very difﬁcult [20]. In order to overcome this weakness, Jian Ling [21] proposed a new ﬁltering method which combines with anisotropic diffusion and median ﬁlter, forming the anisotropic median diffusion ﬁltering method to remove noise. Here, we also propose a novel approach combined with the median ﬁlter, in which the anisotropic diffusion smoothes the region with small gradients and protects the area with large gradients, while the median ﬁlter effectively eliminates the large gradient generated by the noise according to the fact that the median ﬁlter does not affect the gradient generated by the edges. The new nonlocal fuzzy anisotropic diffusion is able to systematically use the global self-prediction conﬁguration information in the image. Fig. 1 presents the 3D plot of the novel nonlocal fuzzy anisotropic diffusion coefﬁcient in (11). When either f is large or ωbj is small at some pixels, the coefﬁcient (| ∇ f|, ωbj ) is small. Then, the diffusion process stops at these pixels, and edges are preserved. In other words, the diffusion process has little effect on those pixels with large difference even though they are very close to the central point in distance.

Q. Zhang et al. / Optik 124 (2013) 4561–4565

4563

Fig. 3. Noisy sinogram and modiﬁed Sheeo–Logan phantom: (a) projection sinogram; (b) modiﬁed Sheep–Logan phantom.

Fig. 2. The 3D plot of ﬂux function.

Fig. 2 shows the 3D plot of the ﬂux function in (12). From the ﬁgure, we observe that when f is large or ωbj is small, the value of the ﬂux function descend rapidly to zero, which indicates that the diffusion process will stop and the feature of the image will be preserved. In contrast, when the difference is small, the diffusion process will have a signiﬁcant effect on the image to remove the noise. 3. Regularized maximum likelihood algorithm based on the fuzzy anisotropic diffusion

we propose a new algorithm, named as RML-NLMFuzzyAD, which incorporates the novel nonlocal fuzzy anisotropic diffusion (NLMFuzzyAD) described in (5)–(10) and a median ﬁlter. Consequently, the result image update procedure becomes fjk,h

I

fjk+1,h =

I

f j

= fjk+1,h +

k+1,h+1

J

aij yi

a f k,h l=1 il l

a i=1 ij i=1

,

k+1,h k+1,h ((|∇ fj,b |, ωbj )∇ fj,b ),

(19)

(20)

b ∈ j

In PET, the measurements follow independent Poisson random distribution. yi is the number of coincidences which are counted by the ith detector during the data collection. f represents the image vector and the element of f denotes the activity of image. A = {aij } is the system matrix. Therefore, we can get the following formulas:

where the diffusion coefﬁcient can be described by (5)–(10) and Median(•) represents the median ﬁltering operation.

yi ∝ Poisson(y¯ i (f )),

4. Simulation and results

i = 1, ..., I

(13)

fjk+1,h+1 = Median(f j

k+1,h+1

),

(21)

J

y¯ i (f ) =

aij fj

(14)

j

P(y|f ) =

I y y¯ (f ) i i

i

exp(−y¯ i (f )),

yi !

(15)

where I is the number of detector pairs, J is the number of the objective image pixels, and P(y|f) is the probability of the detected measurement vector y with image intensity f. The corresponding log-likelihood can be described as follow: L(f ) = log P(y|f ) =

I

J

(yi log(

i=1

j=1

aij fj ) −

J

aij fj )

(16)

j=1

According to the maximum likelihood algorithm, the reconstructed image can be obtained by maximizing the log-likelihood function L(f), i.e., f = argmaxL(f )

(17)

f ≥0

To solve the optimization problem (17), Shepp and Vardi [22] proposed the EM algorithm, i.e., maximum likelihood expectation maximization (MLEM) algorithm. And the iterative formula can be described as follow: fjk+1 =

fjk

I

I

a i=1 ij i=1

aij yi

J

a fk l=1 il l

(18)

With the increment of the number of iteration, the reconstructed image obtained by MLEM suffers noise artifacts. To address the problem, the usual method is to introduce the regularization during the iterative process. In addition, another way suggested by Yan [16] is to employ a proper ﬁlter such as anisotropic diffusion ﬁlter to remove noise in the iterative process. In this paper,

In experiments, we compared the reconstructed results obtained by using the new algorithm with those reconstructed by RML-AD, MLEM, RML-MRP and RML-QM methods to demonstrate the validity of the proposed algorithm. We simulated each sinogram with total count amount to 1 × 106 . Each simulated data is Poisson distributed and all assumed to 128 radial bins with 128 angular views evenly spanned over 180◦ . Fig. 3(a) shows the simulated sinogram and all of these were performed on the modiﬁed Sheep–Logan phantom with the size of 128 × 128 pixels shown in Fig. 3(b). RML-AD combines with the traditional anisotropic diffusion and the diffusion coefﬁcient is cP–M described in (3). The hyper parameter ˇ for regularized maximum likelihoods, threshold parameter k and for algorithms were manually chosen to give the best stable reconstructed image in terms of the tradeoff between noise reduction and edge preservation. The value of ˇ was set to 0.5, 0.000055 for RML-MRP and RML-QM, respectively. We set to 0.1 in RML-NLMFuzzyAD, and 0.22110 in RML-AD. k in RML-NLMFuzzyAD and RML-AD was set to 2.6559 and 1.989, respectively. The size of neighborhoods in RMLNLMFuzzyAD was 3 × 3. Suitable parameter h was set to 20 by hand, and the AD iterative number was set to 15 according to visual errors. The result images reconstructed by different algorithms are shown in Fig. 4. According to Fig. 4, we can observe that both RMLNLMFuzzyAD and RML-AD have the well performance in removing noise, but the former does perverse detail edges accurately, especially the thin edges and protect lots of detail information, while the latter does not. In addition, there are still artifacts observed in the reconstructed image using RML-AD. Comparing with RML-AD, the image reconstructed by RML-NLMFuzzyAD is similar to that of the standard phantom. Furthermore, the edge-preserving performance of the proposed algorithm can be comparable to the well-known MRP algorithm and overcomes its shortcoming that MRP is not capable of suppressing noise and artifact. All of these demonstrate

4564

Q. Zhang et al. / Optik 124 (2013) 4561–4565

Fig. 6. The zoomed images of original image and the zoomed images of the reconstructed images in Fig. 4: (a) the zoomed images of original image; (b) RMLNLMFuzzyAD; (c) RML-AD; (d) RML-MRP; (e) RML-QM; (f) MLEM.

Fig. 4. Reconstructed images by different algorithms: (a) original Image; (b) RMLNLMFuzzyAD; (c) RML-AD; (d) MLEM; (e) RML-MRP; (f) RML-QM.

that the proposed approach is very valid on suppressing noise and preserving edges of image, and it can produce high-quality images. Proﬁles along the 60th row line through the original image and the reconstructions produced by the different methods are illustrated in Fig. 5. This line passes through the edges and different activity regions. We evaluated the performance according to the intensity recovery and edge-preserving capabilities of different algorithms. As depicted in Fig. 5, RML-NLMFuzzyAD can recover image intensity effectively. However, RML-AD does not do better than the proposed algorithm, for it does not perverse the structure

of the edges accurately, especially the thin edges. Our proposed method does not have this shortcoming and we can obtain this in Fig. 5. In order to further test the performance in preserving detail edges capability of the new algorithm, we compared the zoomed regions of interest of the reconstructed image using RMLNLMFuzzyAD with those using RML-AD, RML-QM, and RML-MRP, shown in Fig. 6. The comparisons show that the proposed algorithm (seen in Fig. 6(b)) can be comparable to RML-AD (seen in Fig. 6(c)) in removing noise while RML-MRP does not. From the zoomed images, although the RML-AD reconstruction is characterized by less noise, it does not accurately preserve the detail edges of image and the edges are changed. From Fig. 6(d), we see that RML-MRP can preserve the edges but has apparent step artifacts. Fig. 6(e) and (f) gives the reconstructions using RML-QM and MLEM algorithm, respectively. It is obvious that there is so much noise in the reconstructed images, especially in the one using MLEM. Thus, we obtain that noise effect and staircase artifacts will not be observed in the corresponding reconstructed image by using the proposed method. And it may be conﬁrmed that the proposed method precedes the other algorithms in edge-preserving and intensity recovery. In addition, to evaluate quantitatively the performance of the proposed method, we further computed the signal to noise ratio (SNR) and the root mean square error (RMSE) of the algorithms. (1) Signal-to-noise ratio (SNR) SNR indicates is deﬁned as

⎧ ⎫ ⎨ J (fˆj − f¯ )2 ⎬ j=1 SNR = 10 log 10 , ⎩ J ˆ ¯ 2⎭ j=1

where the numerator and denominator stand for the signal power and the noise power, respectively. fˆ , f¯ and f denote the reconstructed image, the mean of all pixels in fˆ , the corresponding true phantom image, respectively. A bigger SNR value indicates a better performance of the corresponding reconstruction result. (2) Root mean square error (RMSE). This item measures how the reconstructed image is close to the phantom image. That is

RMSE =

Fig. 5. Proﬁles along a line through the original image and the result images produced by different algorithms.

(fj − fj )

J j=1

(fˆj − f¯j ) J2

2

,

where notation fˆj and notation fj denote the value of pixel j of the reconstructed image and that of the original image, respectively. J is the number of pixels of the image. The best algorithm will give the minimum value of root mean square error. Fig. 7 shows the curves of SNR versus the iteration. From Fig. 7, we can get the result that the SNR produced by the proposed method is distinctly higher than those produced by other algorithms. Therefore, the RML-NLMFuzzyAD algorithm signiﬁcantly improves the quality of reconstruction in terms of SNR. That is,

Q. Zhang et al. / Optik 124 (2013) 4561–4565

4565

accurately preserving details in the objective PET images. Visually, the image using the proposed algorithm shows signiﬁcant improvement in terms of image quality. In a word, the new algorithm absorbs the advantage of PDE, fuzzy theory, nonlocal means and median ﬁlter, which assure the proposed algorithm can not only suppress noise effectively, but also preserve the information of edges. Simulated experiments demonstrate that this approach can produce a high quality image which is similar to the original phantom. Acknowledgements This work was supported by the National Natural Science Foundations of China under Grant No. 61071192 and No. 61271357, the Natural Science Foundation of Shanxi Province, China, under Grant No. 2009011020-2, and Program for the Top Young Academic Leaders of Higher Learning Institutions of Shanxi. The authors would like to thank the anonymous reviewers for their constructive comments and suggestions.

Fig. 7. The plot of SNR versus iterations.

References

Fig. 8. The plot of RMSE versus iterations. Table 1 SNRs and RMSEs for the reconstructed images in Fig.4.

SNRs RMSEs

RML-QM

RML-MRP

RML-AD

RML-NLMFuzzyAD

12.1001 13.6086

13.4407 13.3439

13.0948 12.8415

15.8520 9.7289

the RML-NLMFuzzyAD algorithm is superior to the RML algorithms with MRP and traditional AD. In addition, from Fig. 8 which gives the trend curves of RMSE versus the iteration number, we can observe that the RMSE value of the proposed RML-NLMFuzzyAD is obviously smaller than those of other algorithms after 30 iterations, which shows RML-NLMFuzzyAD outperforms other methods. From these quantitative results, it is obvious that the use of RML-NLMFuzzyAD iterative method allows us to obtain better reconstructed images. For further detailed descriptions, we also list the SNRs and RMSEs in Table 1. The SNRs and RMSEs comparisons indicate the proposed reconstruction method can produce superior image than other reconstruction methods. 5. Conclusion In this work, we introduce a novel nonlocal fuzzy anisotropic diffusion into the PET image reconstruction to get a better reconstructed image. At the same time, we combine with median ﬁlter to reduce the impact of the gradient threshold on the reconstructed images. We named the new algorithm as RML-NLMFuzzyAD. From the above experiments and analysis, we can conclude that our proposed method is capable of effectively suppressing noise and

[1] J.A. Fessler, Penalized weighted least-squares image reconstruction for positron emission tomography, IEEE Trans. Med. Imaging 13 (2) (1994) 290–300. [2] J. Zhou, L.M. Luo, Sequential weighted least squares algorithm for PET image reconstruction, Digital Signal Process. 16 (2006) 735–745. [3] P.J. Green, Bayesian reconstruction from emission tomography data using a modiﬁed EM algorithm, IEEE Trans. Med. Imag. 1 (9) (1990) 84–93. [4] W. Chlewicki, F. Hermansen, S.B. Hansen, Noise reduction and convergence of Bayesian algorithms with blobs based on the Huber function and median root prior, Phys. Med. Biol. 49 (2004) 4717–4730. [5] S.J. Lee, A. Rangarajan, G. Gindi, Bayesian image reconstruction in SPECT using higher order mechanical models as priors, IEEE Trans. Med. Imag. 14 (4) (1995) 669–680. [6] L.A. Zadeh, Knowledge representation in fuzzy logic, IEEE Trans. Knowledge Data Eng. 1 (1) (1989) 89–100. [7] V.S. VanDe, M. Nachtegael, D.V.D. Weken, Noise reduction by fuzzy image ﬁltering, IEEE Trans. Fuzzy Syst. 11 (2003) 429–435. [8] F. Russo, An image enhancement techniques combining sharpening and noise reduction, IEEE Trans. Med. Imag. 2 (1982) 113–122. [9] P.P. Modal, K. Rajan, Iterative image reconstruction for emission tomography using fuzzy potential, in: Nuclear Science Symposium Conference Record, IEEE, Roma, Italy, 2004, pp. 3616–3619. [10] P. Perona, J. Malik, Scale space and edge detection using anisotropic diffusion, IEEE Trans. Pattern Anal. Mach. Intell. 12 (7) (1990) 629–639. [11] J. Maeda, T. Iizawa, T. Ishizaka, C. Ishikawa, Y. Suzuki, Segmentation of natural images using anisotropic diffusion and linking of boundary edges, Pattern Recogn. 31 (12) (1998) 1993–1999. [12] O. Demirkaya, Anisotropic diffusion ﬁltering of PET attenuation data to improve emission images, Phys. Med. Biol. 47 (20) (2002) 271–278. [13] Y.L. You, M. Kaveh, Fourth-order partial differential equations for noise removal, IEEE Trans. Image Process. 9 (2000) 1723–1730. [14] R.C. Chen, P.T. Yu, Fuzzy selection ﬁlters for image restoration with neural learning, IEEE Trans. Signal Process. 47 (1999) 1446–1450. [15] A. Buades, B. Coll, J.M. Morel, A review of image denoising algorithms with a new one, Multiscale Model. Simul. 4 (2) (2005) 490–530. [16] J.H. Yan, Investigation of Positron Emission Tomography Image Reconstruction, Huazhong University of Science & Technology, Wuhan, 2007. [17] N. Paragios, Variational methods and partial differential equations in cardiac image analysis, in: International Symposium on Biomedical Imaging: From Nano to Macro, IEEE, Arlinton, 2004, pp. 17–20. [18] S.M. Chao, D.M. Tsai, An improved anisotropic diffusion model for detail and edge preserving smoothing, Pattern Reorg. 31 (2010) 2012–2023. [19] M.J. Black, G. Sapiro, D.H. Marimont, D. Heeger, Robust anisotropic diffusion, IEEE Trans. Image Process. 7 (3) (1998) 421–432. [20] F. Torkamani-Azar, K.E. Tait, Image recovery using the anisotropic diffusion equation, IEEE Trans. Image Process. 5 (1996) 1573–1578. [21] J. Ling, A.C. Bovik, Smoothing low-SNR molecular images via anisotropic median diffusion, IEEE Trans. Med. Imaging 4 (21) (2002) 377–384. [22] L.A. Shepp, Y. Vardi, Maximum likelihood restoration for emission tomography, IEEE Trans. Med. Imaging 1 (1982) 113–122.