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

Authors: Qian He; Lihong Huang
Journal: Mathematical Problems in Engineering
 (2014)

---

## 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.


---

## Body

## 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 Section2, 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 probabilityp
(
g
∣
f
), where g is the measured data and f is the emission image. It is described in the following:
(1)
f
^
=
arg
f
≥
0
max
⁡
L
(
f
)
,
L
(
f
)
=
log
⁡
(
p
(
g
∣
f
)
)
=
∑
i
=
1
N
(
-
∑
j
=
1
M
H
i
,
j
f
j
+
g
i
log
⁡
⁡
(
∑
j
=
1
M
H
i
,
j
f
j
)
)
.
In (2), H
i
j is the probability of photons emitted by pixel j, which can be detected by the detector i, f
j is the number of photons emitted by the pixel j,  and  g
i  is the number of photons captured by the detector i.In order to solve (1), Shepp and Vardi [1] have proposed the EM algorithm. The iterative formula can be described in the following:
(2)
f
j
k
+
1
=
f
j
k
∑
i
=
1
N
H
i
j
(
g
i
/
∑
l
=
1
M
H
i
l
f
l
k
)
∑
i
=
1
N
H
i
j
,
where k 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
(3)
f
^
=
arg
f
≥
0
max
⁡
[
L
(
f
)
+
P
(
f
)
]
,
where L
(
f
) has been explained above and P
(
f
) is regularization term or penalty term.The OSL algorithm uses the current imagef
k when calculating the value of the derivative of the energy function, and the iterative formula can be defined as [6]
(4)
f
j
k
+
1
=
f
j
k
∑
i
=
1
N
H
i
j
(
g
i
/
∑
l
=
1
M
H
i
l
f
l
k
)
∑
i
=
1
N
H
i
j
+
β
(
∂
/
∂
f
j
)
U
(
f
)
|
f
k
,
where U
(
·
) 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]
(5)
f
j
k
+
1
=
f
j
k
∑
i
=
1
N
H
i
j
(
g
i
/
∑
l
=
1
M
H
i
l
f
l
k
)
∑
i
=
1
N
H
i
j
+
β
(
(
f
j
k
-
M
(
f
j
k
)
)
/
M
(
f
j
k
)
)
,
where M
(
f
j
k
) is the median of pixel j 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]
(6)
∂
f
(
x
,
y
,
t
)
∂
t
=
div
⁡
[
g
(
|
∇
f
(
x
,
y
,
t
)
|
)
∇
f
(
x
,
y
,
t
)
]
,
where t is the time parameter, div
⁡ is the divergence operator, f
(
x
,
y
,
0
) is the original image, ∇
f
(
x
,
y
,
t
) is the gradient of the image at time t, and g
(
·
) is the diffusion coefficient, which is a function of local gradient. This function should be satisfied:
(7)
lim
⁡
x
→
0
g
(
x
)
=
1
,
lim
⁡
x
→
∞
g
(
x
)
=
0
,
so that the diffusion is more in smooth regions and less near the edges. They put forward the following two diffusion coefficients:
(8)
g
1
(
x
)
=
exp
⁡
[
-
(
x
K
)
2
]
,
g
2
(
x
)
=
1
1
+
(
x
/
K
)
2
,
where K is a gradient threshold that judges if there is a local edge. Gradient threshold K determines the quality of the filter [18], Perona and Malik suggested using Canny’s “noise estimator” [19] to determine K, Torkamani-Azar and Tait used the mean of the absolute gradient as K [20], and Black et al. determined K 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
(9)
f
j
k
+
1
=
f
j
k
+
w
|
N
j
|
∑
j
′
∈
N
j
(
g
(
|
∇
f
j
,
j
′
k
|
)
∇
f
j
,
j
′
k
)
,
f
j
k
+
1
=
Median
(
f
j
k
,
W
)
,
where w
∈
[
0,1
] controls the rate of diffusion, k is iteration number, f
j is the gray value of pixel j, |
N
j
| is the number of neighbor at pixel j (usually four directions, north, south, east, and west, resp.), ∇
f
j
,
j
′
=
f
j
′
-
f
j, W is the window for the median operator (such as a 3
×
3 square), and the diffusion coefficient g
(
·
) is
(10)
g
(
x
)
=
{
25
16
K
[
1
-
(
x
5
K
)
2
]
2
,
|
x
|
≤
5
K
,
0
,
otherwise.
The flux function ϕ
(
x
) is defined as ϕ
(
x
)
=
g
(
x
)
·
x.The diffusion coefficient and the flux function are presented, respectively, in Figures1 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 withK
=
0.2.Figure 2
The plot of flux function withK
=
0.2.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:(11)
f
j
k
+
1
,
h
=
f
j
k
,
h
∑
i
=
1
N
H
i
j
(
g
i
/
(
∑
l
=
1
M
H
i
l
f
l
k
,
h
)
)
∑
i
=
1
N
H
i
j
,
(12)
f
j
k
+
1
,
h
+
1
=
f
j
k
+
1
,
h
+
w
|
N
j
|
∑
j
′
∈
N
j
(
g
(
|
∇
f
j
,
j
′
k
+
1
,
h
|
)
∇
f
j
,
j
′
k
+
1
,
h
)
,
(13)
f
j
k
+
1
,
h
+
1
=
Median
(
f
j
k
+
1
,
h
+
1
,
W
)
,
where h is diffusion number and diffusion coefficient g
(
·
) 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 Figure3, and the size of the image phantom is 128
×
128 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 6
×
1
0
5.Figure 3
Modified Shepp-Logan phantom.The reconstructed images generated by different algorithms that are shown in Figures4 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.The modified Shepp-Logan phantom reconstructed by different algorithms after 50 iterations: (a) MLEM-AMD withh
=
40 and K
=
1.5; (b) MLEM-PDE with h
=
40 and K
=
40; (c) MRP with β
=
0.1; (d) MLEM.
(a)
(b)
(c)
(d)The zoomed images of Figure4: (a) MLEM-AMD; (b) MLEM-PDE; (c) MRP; (d) MLEM.
(a)
(b)
(c)
(d)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(14)
(
1
)
NRMSE
=
∑
i
=
1
N
∑
j
=
1
M
[
f
(
i
,
j
)
-
f
^
(
i
,
j
)
]
2
∑
i
=
1
N
∑
j
=
1
M
[
f
(
i
,
j
)
]
2
,
where f
^
(
i
,
j
) and f
(
i
,
j
) denote the value of pixel (
i
,
j
) in reconstructed image and the original image, respectively. The smaller the value of NRMSE is, the better the performance of the algorithm embodies:
(15)
(
2
)
SNR
=
10
lo
g
10
∑
i
=
1
N
∑
j
=
1
M
(
f
^
(
i
,
j
)
-
f
-
)
2
∑
i
=
1
N
∑
j
=
1
M
(
f
^
(
i
,
j
)
-
f
(
i
,
j
)
)
2
,
where f
^
(
i
,
j
) and f
(
i
,
j
) have been explained above and f
- is the average grayscale of all pixels in reconstructed image f
^. The better the algorithm is, the bigger the value of SNR is.Figure6 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 Figure8) 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
×
1
0
5.Figure 8
Real thorax phantom.The images reconstructed by MLEM-AMD, MLEM-PDE, MRP, and MLEM that are shown in Figures9 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.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)The zoomed images of Figure9: (a) MLEM-AMD; (b) MLEM-PDE; (c) MRP; (d) MLEM.
(a)
(b)
(c)
(d)Figures11 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 Section5, 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.



---

*Source: 491239-2014-01-20.xml*