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

Nowadays, positron emission tomography (PET) iswidely used in engineering. In this paper, a novel penalizedmaximum likelihood (PML) algorithm is presented for improving the quality of PET images. The proposed algorithm fuses an anisotropic mediandiffusion (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.


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][4][5][6][7][8][9][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][13][14][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.

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 backprojection (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.

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, div is the divergence operator, (, , 0) 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  ∈ [0, 1] 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 3 × 3 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.
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].

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.

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 computergenerated modified Shepp-Logan phantom are shown in Figure 3, and the size of the image phantom is 128 × 128 pixels.The distribution 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 6 × 10 5 .
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.
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 f(, ) 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 f(, ) and (, ) have been explained above and  is the average grayscale of all pixels in reconstructed image f.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 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.
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 5.2 × 10 5 .
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.
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.

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.

Figure 2 :
Figure 2: The plot of flux function with  = 0.2.

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

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

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