Mathematical Problems in Engineering

Volume 2014 (2014), Article ID 491239, 7 pages

http://dx.doi.org/10.1155/2014/491239

Research Article

Penalized Maximum Likelihood Algorithm for Positron Emission Tomography by Using Anisotropic Median-Diffusion

Qian He1,2 and Lihong Huang1,3

¹College of Mathematics and Econometrics, Hunan University, Changsha 410082, China

²College of Information Science and Engineering, Hunan City University, Yiyang 413000, China

³Department of Information and Technology, Hunan Women’s University, Changsha 410004, China

Received 1 October 2013; Revised 9 December 2013; Accepted 24 December 2013; Published 20 January 2014

Academic Editor: Gelan Yang

Copyright © 2014 Qian He and Lihong Huang. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

Nowadays, positron emission tomography (PET) is widely used in engineering. In this paper, a novel penalized maximum likelihood (PML) algorithm is presented for improving the quality of PET images. The proposed algorithm fuses an anisotropic median-diffusion (AMD) filter to the maximum-likelihood expectation-maximization (MLEM) algorithm. The fusing algorithm shows its positive effect on image reconstruction and denoising. Experimental results present that the proposed method denoises and reconstructs images with high quality. Furthermore, by comparing with other classical reconstructing algorithms, this novel algorithm shows better performance in the edge preservation.

1. Introduction

PET technology, which has been widely used in neurology, oncology, and new medicine exploitation, is one of the advanced and noninvasive diagnostic techniques in modern nuclear medical. In order to obtain a high quality reconstructed image from clinical projection data with strong noise, an excellent image reconstruction algorithm is indispensable.

The MLEM algorithm is a classic method in PET image reconstruction when the measured data follows Poisson distribution [1]. One problem of this algorithm is the ill-posed problem, which represents that the reconstructed images cannot remove the noise of projection data [2]. Today, an ill-posed image reconstruction problem, such as MLEM, can be transformed into a well-posed one through the use of regularization term. The reconstructed results should be not only content with measured data to some extent but also be consistent with additional regularization term that is independent of those data at the same time. That is usually called PML or Bayesian algorithm. Numerous PML algorithms have been proposed in the past decades [3–10].

Thereinto, Green proposed a Bayesian algorithm, known as the one-step-late (OSL) algorithm [6]. The key of this algorithm is to find an appropriate energy function, which is defined by Gibbs probability distribution. Unfortunately, the selection of the energy function is difficult. The median root prior (MRP) algorithm [9], firstly proposed by Alenius, is an application of OSL algorithm. This algorithm is good at coping with those images that have locally monotonic structures. However, the images reconstructed by MRP are still noisy because median filter cannot remove Gaussian and Poisson noise effectively, which dominate in PET images [7].

The anisotropic diffusion (AD) filter [11] is a nonlinear partial differential equation (PDE) based on diffusion process. Overcoming the undesirable effects of linear smoothing filter, such as blurring or dislocating the useful edge information of the images, AD has been widely used in image smoothing, image reconstruction and image segmentation [12–15]. The basic idea of AD is adaptively choosing diffusive coefficients in diffusion process so that diffusion is maximal within smooth regions and minimal near the edges.

In order to remove noise and preserve edge information at the same time, image reconstruction based on AD has become the research focus [7, 8, 15, 16]. Yan proposed a PML algorithm [16] that combined MLEM with AD filter (called MELM-PDE) and could obtain acceptable reconstructed results. However, MLEM-PDE cannot remove the isolated noise and preserve edge information accurately due to the defects of P-M diffusion model [15].

In this paper, we proposed a new PML algorithm for PET image reconstruction based on AMD. The proposed algorithm can effectively remove noise while preserving edge information accurately. In Section 2, PET image reconstruction algorithms such as MLEM, OSL, and MRP are introduced. The AMD filter is presented in Section 3. In Section 4, our proposed algorithm is described. Simulation experiments are given in Section 5. Finally, Section 6 is the conclusion.

2. Image Reconstruction Algorithms for PET

In PET, the maximum likelihood (ML) algorithm seeks a solution that makes the measured data most likely to occur and maximizes the conditional probability , where is the measured data and is the emission image. It is described in the following:

In (2), is the probability of photons emitted by pixel , which can be detected by the detector , is the number of photons emitted by the pixel ,  and   is the number of photons captured by the detector .

In order to solve (1), Shepp and Vardi [1] have proposed the EM algorithm. The iterative formula can be described in the following:

where is the iteration.

Although MLEM algorithm is better than filtered back-projection (FBP) algorithm [17], its convergence rate is extremely slow, and as the iteration number increases, the reconstructed results suffer from noise artifacts. The usual method to solve this problem is to introduce a regularization term, and the objective function is

where has been explained above and is regularization term or penalty term.

The OSL algorithm uses the current image when calculating the value of the derivative of the energy function, and the iterative formula can be defined as [6]

where is the energy function and is the Bayes weight of the prior.

The MRP algorithm can be coped with monotonic structures in a neighborhood by comparing the pixel against the local median. Edge preservation is an intrinsic characteristic of median filter, and its update equation is [9]

where is the median of pixel within its neighborhood.

3. Anisotropic Median-Diffusion Filter

The AD filter (usually called P-M diffusion model), firstly proposed by Perona and Malik, is a nonlinear filter that purports to remove noise without blurring edges, and the basic equation is [11]

where is the time parameter, is the divergence operator, is the original image, is the gradient of the image at time , and is the diffusion coefficient, which is a function of local gradient. This function should be satisfied:

so that the diffusion is more in smooth regions and less near the edges. They put forward the following two diffusion coefficients:

where is a gradient threshold that judges if there is a local edge. Gradient threshold determines the quality of the filter [18], Perona and Malik suggested using Canny’s “noise estimator” [19] to determine , Torkamani-Azar and Tait used the mean of the absolute gradient as [20], and Black et al. determined from the median absolute deviation [21].

Although AD filter can suppress the noise and preserve the edge information to same extent, it cannot preserve the detail edges effectively and accurately [8]. Ling and Bovik proposed AMD filter [22] on the basis of AD filter. This filter incorporates a median filter into the diffusion step, and the discrete form is defined as

where controls the rate of diffusion, is iteration number, is the gray value of pixel , is the number of neighbor at pixel (usually four directions, north, south, east, and west, resp.), , is the window for the median operator (such as a square), and the diffusion coefficient is

The flux function is defined as .

The diffusion coefficient and the flux function are presented, respectively, in Figures 1 and 2. These figures together with (9) suggest that diffusion is maximal within smooth regions and stopped near the edges.

Figure 1: The plot of diffusion function with .

Figure 2: The plot of flux function with .

In AMD, diffusion process can smooth the regions with small gradient between current pixel and its neighborhood, while regions with large gradient (edge or noise spike) are left unchanged. Noise spike that generated the large gradients will be removed effectively by the median filter subsequently. However, those large gradients generated by the edges will not be affected by the median filter. Thus, low-level noise is smoothed by diffusion, and high-level noise is removed by median filter. In brief, the AD and median filter are always mutually complementary to gradually eliminate noise without blurring the edges [22].

4. Proposed Algorithm

Although MLEM-PDE algorithm is better than classical reconstructing algorithms, the images reconstructed by MLEM-PDE are still noisy because P-M diffusion model is not good enough for the removal of isolated noise and accurate preservation of edges. Therefore, it is natural that a novel diffusion model should be introduced into the MLEM algorithm. In this paper, we proposed a new PML algorithm (called MELM-AMD) by fusing AMD filter to the MLEM algorithm. The discrete form of the proposed algorithm can be defined as follows:

where is diffusion number and diffusion coefficient is defined by (10).

The proposed algorithm can select appropriate diffusion coefficients adaptively during the diffusion process. At each diffusion step, noise with small gradients will be smoothed by (12), while the large gradients generated by noise spike will be removed effectively by (13). If large gradients generated by edges, the algorithm will not affect them. So, the proposed algorithm can preserve the edges more validly than other methods.

5. Simulation Experiments

In order to examine the effectiveness of the proposed algorithm, we compared the reconstructed image produced by MLEM-AMD algorithm with those brought by different algorithms, such as MLEM, MRP, and MLEM-PDE in computer simulations. These algorithms tested by computer-generated modified Shepp-Logan phantom are shown in Figure 3, and the size of the image phantom is pixels. The distribution of projection data has the Poisson characteristics and is assumed to be generated by 128 angular views (averagely distributed in the range of 180 degrees), while each angle has 128 radial bins. The Poisson noise is added into projection data and the number of photons captured by radial bins is about .

Figure 3: Modified Shepp-Logan phantom.

The reconstructed images generated by different algorithms that are shown in Figures 4 and 5 are the zoomed images of Figure 4. From the figures, we can see that the proposed algorithm has well performance of noise removal and edges preservation especially the thin edges and detailed information. At the same time, we can observe that the MLEM-AMD algorithm overcomes the shortcoming of streak artifacts and the reconstructed image is more similar to the original phantom. In a word, the proposed algorithm is good at producing high quality image with low noise and preserving edge information.

(a)

(b)

(c)

(d)

Figure 4: The modified Shepp-Logan phantom reconstructed by different algorithms after 50 iterations: (a) MLEM-AMD with and ; (b) MLEM-PDE with and ; (c) MRP with ; (d) MLEM.

(a)

(b)

(c)

(d)

Figure 5: The zoomed images of Figure 4: (a) MLEM-AMD; (b) MLEM-PDE; (c) MRP; (d) MLEM.

Normalized root-mean-square error (NRMSE) and signal-to-noise ratios (SNR) are used to evaluate quantitatively the performance of the proposed algorithm, which are defined as

where and denote the value of pixel in reconstructed image and the original image, respectively. The smaller the value of NRMSE is, the better the performance of the algorithm embodies:

where and have been explained above and is the average grayscale of all pixels in reconstructed image . The better the algorithm is, the bigger the value of SNR is.

Figure 6 is the plots of NRMSE along with iterations for different algorithms. As it is shown, the NRMSE produced by the proposed algorithm is always smaller than those produced by other algorithms distinctly. The MLEM-AMD algorithm significantly improves the quality of reconstructed image in term of NRMSE. Figure 7 is the plots of SNR along with iterations for different algorithms, and it also proves that the MLEM-AMD algorithm can obtain the better reconstructed image than other algorithms.

Figure 6: The plots of NRMSE along with iterations for different algorithms.

Figure 7: The plots of SNR along with iterations for different algorithms.

Simulated data generated by the real thorax phantom (see Figure 8) is applied to further test the validity of proposed algorithm. The projection data were generated by the method mentioned above, and the number of photons captured by radial bins is about .

Figure 8: Real thorax phantom.

The images reconstructed by MLEM-AMD, MLEM-PDE, MRP, and MLEM that are shown in Figures 9 and 10 are the zoomed images of Figure 9. From the figures, we can see that proposed algorithm can also remove noise and preserve edge information effectively.

(a)

(b)

(c)

(d)

Figure 9: The real thorax phantom reconstructed by different algorithms after 50 iterations: (a) MLEM-AMD; (b) MLEM-PDE; (c) MRP; (d) MLEM.

(a)

(b)

(c)

(d)

Figure 10: The zoomed images of Figure 9: (a) MLEM-AMD; (b) MLEM-PDE; (c) MRP; (d) MLEM.

Figures 11 and 12 are the plot of NRMSE and SNR along with iteration number for different algorithms, respectively. From these quantitative results, it is also obvious that the application of proposed algorithm can obtain better reconstructed image.

Figure 11: The plots of NRMSE along with iterations for different algorithms.

Figure 12: The plots of SNR along with iterations for different algorithms.

6. Conclusion

In this paper, a novel PML algorithm called MLEM-AMD was proposed for improving the quality of reconstructed image, which fuses AMD filter into MLEM. Iterations of the proposed method can be divided into two steps: firstly, reconstruct image with MLEM algorithm; secondly, remove noise with the AMD filter. From simulation experiments in Section 5, we can see that the MLEM-AMD algorithm can remove noise and preserve edge information effectively in PET image reconstruction. In our experiment, MLEM-AMD shows better performance than the classic algorithms, for example, MLEM-PDE, MRP, and MLEM. The reconstructed images by the proposed algorithm show significant improvement in visual quality. In conclusion, the proposed algorithm absorbs the advantage of AMD filter and can reconstruct the images of high quality.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgments

The authors would like to thank anonymous reviewers for their helpful work of this paper. This work was supported by Fund Project of Hunan Provincial Education Department no. 12C0573 and Yiyang Municipal Science and Technology Bureau no. 2012JZ50.

References

L. A. Shepp and Y. Vardi, “Maximum likelihood reconstruction for emission tomography,” IEEE Transactions on Medical Imaging, vol. 1, no. 2, pp. 113–122, 1982.

A. Gaitanis, G. Kontaxakis, G. Spyrou, G. Panayiotakis, and G. Tzanakos, “PET image reconstruction: a stopping rule for the MLEM algorithm based on properties of the updating coefficients,” Computerized Medical Imaging and Graphics, vol. 34, no. 2, pp. 131–141, 2010.

C.-T. Chen, V. E. Johnson, W. H. Wong, X. Hu, and C. E. Metz, “Bayesian image reconstruction in positron emission tomography,” IEEE Transactions on Nuclear Science, vol. 37, no. 1, pp. 636–641, 1990.

S. Somayajula, C. Panagiotou, A. Rangarajan, Q. Li, S. R. Arridge, and R. M. Leahy, “PET image reconstruction using information theoretic anatomical priors,” IEEE Transactions on Medical Imaging, vol. 30, no. 3, pp. 537–549, 2011.

G. Wang and J. Qi, “Penalized likelihood PET image reconstruction using patch-based edge-preserving regularization,” IEEE Transactions on Medical Imaging, vol. 31, no. 12, pp. 2194–2204, 2012.

P. J. Green, “Bayesian reconstructions from emission tomography data using a modified EM algorithm,” IEEE Transactions on Medical Imaging, vol. 9, no. 1, pp. 84–93, 1990.

J. H. Yan and J. Yu, “Median-prior tomography reconstruction combined with nonlinear anisotropic diffusion filtering,” Journal of the Optical Society of America A, vol. 24, no. 4, pp. 1026–1033, 2007.

Z. G. Gui and J. W. He, “Regularized maximum likelihood algorithm for PET image reconstruction using a detail and edges preserving anisotropic diffusion,” Optik, vol. 123, no. 6, pp. 507–510, 2012.

S. Alenius and U. Ruotsalainen, “Bayesian image reconstruction for emission tomography based on median root prior,” European Journal of Nuclear Medicine, vol. 24, no. 3, pp. 258–265, 1997.

I.-T. Hsiao, A. Rangarajan, and G. Gindi, “A new convex edge-preserving median prior with applications to tomography,” IEEE Transactions on Medical Imaging, vol. 22, no. 5, pp. 580–585, 2003.

P. Perona and J. Malik, “Scale-space and edge detection using anisotropic diffusion,” IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 12, no. 7, pp. 629–639, 1990.

K. Krissian and S. Aja-Fernández, “Noise-driven anisotropic diffusion filtering of MRI,” IEEE Transactions on Image Processing, vol. 18, no. 10, pp. 2265–2274, 2009.

D. Kazantsev, S. R. Arridge, S. Pedemonte, et al., “An anatomically driven anisotropic diffusion filtering method for 3D SPECT reconstruction,” Physics in Medicine and Biology, vol. 57, no. 12, pp. 3793–3810, 2012.

Y. Wu, B. Tracey, P. Natarajan, and J. P. Noonan, “James-Stein type center pixel weights for non-local means image denoising,” IEEE Signal Processing Letters, vol. 20, no. 4, pp. 411–414, 2013.

Z. G. Gui, Y. Liu, and J. W. He, “PML algorithm for positron emission tomography combined with nonlocal fuzzy anisotropic diffusion filtering,” IEEE Transactions on Nuclear Science, vol. 59, no. 5, pp. 1984–1989, 2012.

J. H. Yan, Investigation of Positron Emission Tomography Image Reconstruction, Huazhong University of Science & Technology, Wuhan, China, 2007.

G. T. Herman, Image Reconstruction from Projections: the Fundamentals of Computerized Tomography, Academic Press, New York, NY, USA, 1980.

C. Tsiotsios and M. Petrou, “On the choice of the parameters for anisotropic diffusion in image processing,” Pattern Recognition, vol. 46, no. 5, pp. 1369–1381, 2013.

J. Canny, “A computational approach to edge detection,” IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 8, no. 6, pp. 679–698, 1986.

F. Torkamani-Azar and K. E. Tait, “Image recovery using the anisotropic diffusion equation,” IEEE Transactions on Image Processing, vol. 5, no. 11, pp. 1573–1578, 1996.

M. J. Black, G. Sapiro, D. H. Marimont, and D. Heeger, “Robust anisotropic diffusion,” IEEE Transactions on Image Processing, vol. 7, no. 3, pp. 421–432, 1998.

J. Ling and A. C. Bovik, “Smoothing low-SNR molecular images via anisotropic median-diffusion,” IEEE Transactions on Medical Imaging, vol. 21, no. 4, pp. 377–384, 2002.