PET Image Reconstruction with Multiple Kernels and Multiple Kernel Space Regularizers

Let $\bm{p}=[p_{1},\cdots,p_{i},\cdots,p_{N}]^{T}$ be the PET projection data vector with its elements being Poisson random variables having the log-likelihood function [18]

where $N$ is the total number of lines of response and $\bar{p}_{i}$ is the expectation of $p_{i}$. The expectation of $\bm{p}$, $\bar{\bm{p}}=[\bar{p}_{1},\cdots,\bar{p}_{i},\cdots,\bar{p}_{N}]^{T}$, is related to the unknown emission image $\bm{x}=[x_{1},\cdots,x_{i},\cdots,x_{M}]^{T}$ by the forward projection

where $x_{i}$ is the $i$th voxel intensity of the image $\bm{x}$, $M$ is the total number of voxels, $\bm{H}\in{\cal R}^{N\times M}$ is the PET system matrix, and $\bm{r}=[r_{1},\cdots,r_{i},\cdots,r_{N}]^{T}$ is the expectation of random and scattered events.

An estimate of the unknown image $\bm{x}$ can be obtained by the maximum a posteriori (MAP) estimation of $\bm{x}$ [4, 1]

Here $P(\bm{p})$ is a normalizing constant independent of $\bm{x}$; $P(\bm{x})$ represents the prior knowledge of unknown image $\bm{x}$ and is assumed to have the Gibbs distribution [29]

where $W$ is a normalizing constant and $V(x)$ the energy function. In this paper, we consider $V(\bm{x})$ as a weighted sum of multiple prior functions $v_{i}(\bm{x})$

with the weighting parameters $\beta_{i}$.

Solving $\hat{\bm{x}}$ from (3) is equivalent to maximizing the log-posterior probability with

where $L(\bm{p}|\bm{x})$ is as defined in (1), and the contributions of $P(\bm{p})$ and $W$ have been excluded since they are independent of $\bm{x}$. When all $\beta_{i}=0$, (6) reduces to the ML estimation [18]. The multiple prior functions $\beta_{i}v_{i}(\bm{x})$ in $V(\bm{x})$ introduce multiple regularizers in the ML estimation.

The solution of (6) can be readily found by the regularized MLEM algorithm [10, 4, 13]

where $\partial v_{i}(\bm{x}^{n})/\partial\bm{x}^{n}:=[\partial v_{i}({\bm{x}})/\partial{\bm{x}}]_{{\bm{x}}=\bm{x}^{n}}$, $\bm{x}^{n}$ is the estimate of $\bm{x}$ at the $n$th iteration, $\bm{1}_{N}\in R^{N}$ is an all 1 vector, $\circ$ and $\oslash$ are Hadamard product and division, defined as $A\circ B=[a(i,j)b(i,j)]$ and $A\oslash B=[a(i,j)/b(i,j)]$ for the matrices $A=[a(i,j)]$ and $B=[b(i,j)]$ of the same dimension.

The algorithm (7) includes the MLEM algorithm as a special case when all $\beta_{i}=0$. The summation term of $(\beta_{i}\partial v_{i}(\bm{x}^{n})/\partial\bm{x}^{n})$’s in (7) regularizes the MLEM algorithm, which greatly accelerates convergence and reduces the variance at convergence [1, 4] as discussed before.

Using the kernel theory [30], the unknown emission image $\bm{x}$ can be presented in the kernel space as [18]

where $\bm{K}\in{\cal R}^{M\times M}$ is a kernel (basis) matrix and ${\bm{a}}\in{\cal R}^{M}$ is the coefficient vector of $\bm{x}$ under $\bm{K}$.

The $ij$th element of $\bm{K}$ is given by $K(i,j)=k(\bm{\tilde{x}}_{i},\bm{\tilde{x}}_{j})=\left\langle\phi(\bm{\tilde{x}}_{i}),\phi(\bm{\tilde{x}}_{j})\right\rangle,i,j\in\{1,2,\cdots,M\}$, where $k(\bm{\tilde{x}}_{i},\bm{\tilde{x}}_{j})$ is a kernel (function), $\bm{\tilde{x}}_{i}$ and $\bm{\tilde{x}}_{j}$ are respectively the data vectors for the $i$th and $j$th elements of $\bm{x}$ in the training dataset $\bm{\tilde{X}}$ of $\bm{x}$ obtained from other sources, eg MR anatomical images of the same (class of) subject/s, $\phi(\cdot)$ is a nonlinear function mapping $\bm{x}$ to a higher dimensional feature space but needs not to be explicitly known, and $\left\langle\cdot,\cdot\right\rangle$ denotes the inner product.

Typical kernels include the radial Gaussian kernel and the polynomial kernel

where $\sigma$, $\gamma$ and $d$ are kernel parameters, and $\|\cdot\|$ is the 2-norm of a vector. For computation simplicity and performance improvement, the kernel matrix $\bm{K}$ in practice is often constructed using

where $\cal{J}$ is the neighborhood of the $i$th element of emission image $\bm{x}$, with the neighborhood size $J$. Other methods such as $\bm{\tilde{x}}_{j}\in kNN$ (the $k$ nearest neighbor) of $\bm{\tilde{x}}_{i}$, measured by $\|\bm{\tilde{x}}_{i}-\bm{\tilde{x}}_{j}\|$, can also be used [18].

Because the kernel matrix $\bm{K}$ constructed in (11) involves only a single kernel $k(\bm{\tilde{x}}_{i},\bm{\tilde{x}}_{j})$, we will call it the single-kernel matrix hereafter.

Using the kernel space image model (8), we can convert the forward projection (2) and the log-posterior probability maximization (6) to their respective kernel space equivalents

and obtain a single-kernel regularized EM (KREM) algorithm

where $V_{i}(\bm{a}^{n})=v_{i}(\bm{Ka}^{n})$ with $\bm{K}$ absorbed into its parameters. We can then use (14) to estimate the coefficient $\bm{a}$ and use $\bm{x}^{n}=\bm{Ka}^{n}$ to obtain an estimate of the image $\bm{x}$.

The KREM algorithm (14) includes the KEM algorithm of [18, 20] as a special case when all $\beta_{i}=0$. As shown in Sections 3 and 4, the KEM algorithm tends to have larger reconstruction error and is sensitive to the number of iterations and the neighborhood size $J$ used in (11) to construct the single-kernel matrix $\bm{K}$. These problems may significantly degrade KEM’s performance in PET image reconstruction, especially at the low counts. Whereas, the KREM algorithm (14) with regularizations can effectively overcome these problems.

The single-kernel matrix $\bm{K}$ in (14) is the core of KREM algorithm. Its effect can be tuned by the selection of kernel $k(\bm{x}_{i},\bm{x}_{j})$, kernel parameters, and neighborhood size $J$ in (11). However, the efficacy of such tuning is limited as a single kernel is often insufficient for complex problems, such as heterogeneous data [31]. To solve this problem, a multi-kernel approach has been developed in machine learning, which combines multiple base kernels, such as that of (9) or (10), to construct more complex kernel $k_{M}(\bm{x}_{i},\bm{x}_{j})$ and corresponding multi-kernel matrix $\bm{K}_{M}$ for significantly improved performance [32, 31].

Inspired by the success of multi-kernel approach in machine learning, we introduce this approach in PET image reconstruction. We construct the multi-kernel matrix $\bm{K}_{M}$ by

and replace the single-kernel matrix $\bm{K}$ with $\bm{K}_{M}$ in KREM algorithm (14) to derive the multi-kernel KREM (MKREM) algorithm for improved performance.

In (15), $\bm{K}_{l}=[K_{l}(i,j)]=[k_{l}(\bm{x}_{i},\bm{x}_{j})]$, $i,j=1,2,\cdots,M$, are constructed using (11) with different or same $k(\bm{x}_{i},\bm{x}_{j})$ having different or same parameters and neighborhood size $J$ for different $l\in\{1,2,\cdots,G\}$; and $\bm{K}_{M}=[K_{m}(i,j)]$, $i,j=1,2,\cdots,M$, is the product of $\bm{K}_{l}$’s thus constructed.

The multi-kernel matrix $\bm{K}_{M}$ constructed using (15) includes single-kernel matrix $\bm{K}$ as a special case when $G=1$. Hence, the MKREM algorithm derived from replacing $\bm{K}$ by $\bm{K}_{M}$ in (14) includes KREM algorithm as a special case. We therefore call (14) the MKREM algorithm and KREM algorithm alternately when there is no confusion.

It can be readily proven that for all $i,j\in\{1,2,\cdots,M\}$ the $ij$th element of $\bm{K}_{M}$, $K_{M}(i,j)=k_{M}(\bm{\tilde{x}}_{i},\bm{\tilde{x}}_{j})$, is a new kernel resulting from multiplications and summations of the base kernels $k_{l}(\bm{\tilde{x}}_{i},\bm{\tilde{x}}_{j})$, $l=1,2,\cdots,G$. Therefore, $\bm{K}_{M}$ constructed from (15) is indeed a multi-kernel matrix [32, 31]. The proof is omitted due to space limit.

The regularizers $V_{i}(\bm{a}^{n})$ in (14) are symbolic. In this and the next subsections, we will present two specific regularizers to be used in this work. To facilitate presentation, we introduce the kernel space image dictionary for PET emission image $\bm{x}$.

Let $\bm{K}_{al}$, $\bm{K}_{bl}\in R^{M\times M}$, $l=1,2,\cdots,G$, be two groups of single-kernel matrices constructed using (11) with the kernel $k_{l}(\cdot,\cdot)$ and neighborhood size $J=J_{al}$ for $\bm{K}_{al}$ and $J=J_{bl}$ for $\bm{K}_{bl}$. Let $\bm{K}_{Ma}$, $\bm{K}_{Mb}\in R^{M\times M}$ be two multi-kernel matrices constructed using (15) with $\bm{K}_{al}$ for $\bm{K}_{Ma}$ and $\bm{K}_{bl}$ for $\bm{K}_{Mb}$. Here $\bm{K}_{Ma}$ is used as $\bm{K}_{M}$ in MKREM algorithm (14) and $\bm{K}_{Mb}$ is used in the dictionary learning described below. When $G=1$, $\bm{K}_{Ma}=\bm{K}_{a1}:=\bm{K}_{a}$ and $\bm{K}_{Mb}=\bm{K}_{b1}:=\bm{K}_{b}$ reduce to the single-kernel matrices. As shown by the experiment results in Section 3, $J_{b}<J_{a}$ and hence $\bm{K}_{Mb}\neq\bm{K}_{Ma}$ are generally needed for the best reconstruction performance.

Since $\bm{K}_{Ma}$, $\bm{K}_{Mb}\in R^{M\times M}$ and $\bm{K}_{Mb}\neq\bm{K}_{Ma}$, $\bm{K}_{Mb}$ admits the kernel matrix factorization

Replacing $\bm{K}$ in (8) by $\bm{K}_{Ma}$ and $\bm{K}_{Mb}$ respectively gives the representations of emission image $\bm{x}$ under the multi-kernel matrices $\bm{K}_{Ma}$ and $\bm{K}_{Mb}$

where $\bm{a},\bm{b}\in R^{M}$ are respectively the coefficient vectors of $\bm{x}$ under $\bm{K}_{Ma}$ and $\bm{K}_{Mb}$ and are related by $\bm{a}=\bm{\tilde{K}}\bm{b}$.

The $\bm{a}$ and $\bm{b}$ can be viewed as kernel space image signals and hence both admit dictionary representation [33]. Defining a dictionary $\bm{D}_{b}\in R^{M\times S}$ with $S$ $(S\geq M)$ dictionary atoms and a coefficient vector $\bm{c}\in R^{S}$ for $\bm{b}$, we can represent $\bm{b}$ and $\bm{a}$ respectively as

where $\bm{D}_{a}:=\bm{\tilde{K}}\bm{D}_{b}$ is a dictionary for $\bm{a}$. It can be seen that $\bm{\tilde{K}}$ is a mapping of $\bm{b}$ and $\bm{D}_{b}$ in the kernel space of $\bm{K}_{Mb}$ to $\bm{a}$ and $\bm{D}_{a}$ in the kernel space of $\bm{K}_{Ma}$.

Let $\bm{\tilde{X}}=[\bm{\tilde{x}}_{1},\cdots,\bm{\tilde{x}}_{i},\cdots,\bm{\tilde{x}}_{T}]\in R^{M\times T}$ be the matrix consisting of $T$ training data vectors $\bm{\tilde{x}}_{i}$ for $\bm{x}$, obtained from other sources, eg MR anatomical images of the same (class of) subject/s. We use the following sparse dictionary learning to obtain the dictionary $\bm{D}_{b}$ [33]

where $\bm{B}:=\bm{K}_{Mb}^{-1}\bm{\tilde{X}}$ is the data matrix for dictionary learning, $\bm{C}:=[\bm{c}_{1},\cdots,\bm{c}_{i},\cdots,\bm{c}_{T}]$ is the matrix consisting of $T$ estimates of the coefficient vector $\bm{c}$, $\|\cdot\|_{F}$ is the Frobenius norm, and $s$ is a sparsity constraint constant. Theoretically, $\bm{K}_{Mb}$ is always positive semidefinite [30]. In the rare case it is singular, $(\bm{K}_{Mb}+\epsilon I)^{-1}$, with $\epsilon>0$ a small constant, can be used as $\bm{K}_{Mb}^{-1}$ without affecting image reconstruction result.

In case the number of training data vectors $T$ is limited, eg $T=1$, or the row dimension of $\bm{B}$ is too high, we can use the patch-based dictionary learning to increase training data or reduce computation complexity [33]

where $\mathbb{P}$ is the operator extracting patches from $\bm{B}$ and $\bm{D}_{b}\bm{C}$.

A number of well developed methods can be used to solve (17) ((18)) for $\bm{D}_{b}$ and $\bm{C}$ [33]. All these methods are based on the alternating optimization: 1) fix $\bm{D}_{b}$ and solve (17) ((18)) for $\bm{C}$, 2) fix $\bm{C}$ and solve (17) ((18)) for $\bm{D}_{b}$, 3) repeat 1) and 2) until the maximum number of iterations is reached. In this work, we use the orthogonal matching pursuit (OMP) in 1) to compute $\bm{C}$ and the K-SVD in 2) to compute $\bm{D}_{b}$ [33].

With above learned dictionary $\bm{D}_{b}$, we can obtain

and introduce the first regularizer and its derivative

where $\bm{a}^{n}$ is from the $n$th MKREM iteration in (14), and $\bm{c}^{n}$ is from solving the optimization below for $\bm{c}^{n}$ using $\bm{a}^{n}$ and OMP before the $(n+1)$th MKREM iteration.

As mentioned above, the neighborhood size $J_{b}<J_{a}$ is generally needed for the multi-kernel matrix $\bm{K}_{Mb}$ used in dictionary learning. This is to allow the dictionary to capture the local texture details embedded in the training image data which is usually the images with anatomical details. To avoid over-localization of the image reconstruction, we treat the image data as the signal on graph and use the graph Laplacian quadratic [34] of the estimated kernel image $\bm{a}^{n}$, given in (23), as the second regularizer to smooth the estimated PET image $\bm{x}^{n}=\bm{K}_{Ma}\bm{a}^{n}$ on the image graph at each iteration of MKREM algorithm (14).

with the derivative

where $\bm{Q}_{a}:=\bm{K}_{Ma}^{T}\bm{Q}\bm{K}_{Ma}$ and $\bm{Q}$ is the Laplacian matrix of the graph of training image data described below.

Inspired by [35], we treat a training image $\bm{\tilde{x}}=[\tilde{x}_{1},\cdots,\tilde{x}_{i},\cdots,\tilde{x}_{M}]^{T}$ as a graph signal and construct its Laplacian matrix as follows: We first use a window of size $m$ centered at the $i$th element $\tilde{x}_{i}$ of $\bm{\tilde{x}}$ to extract $m$ elements and use them to form the vector $\bm{y}_{i}$, for $i=1,2,\cdots,M$. Then we construct a distance matrix $\bm{E}$ with elements $E(i,j)=\|\bm{y}_{i}-\bm{y}_{j}\|$, $i,j=1,2,\cdots,M$, and convert $\bm{E}$ to an affinity matrix $\bm{W}$ by an adaptive Gaussian function $W(i,j)=exp(-E(i,j)/b_{i})^{2}$, where $b_{i}$ is the ($knn+1$) smallest $E(i,j)$ in the $i$th row of $\bm{E}$. Next, we symmetrize $\bm{W}$ to $\bm{\tilde{W}}=\bm{W}+\bm{W}^{T}$ and construct a Markov transition matrix $\bm{Z}$ with $Z(i,j)=\tilde{W}(i,j)/\sum_{j}\tilde{W}(i,j)$. If there are $T>1$ training images available, we apply above procedure to each image $\bm{\tilde{x}}_{i}$ to get $\bm{Z}_{i}$ and take $\bm{Z}=\frac{1}{T}\sum_{i=1}^{T}\bm{Z}_{i}$.

Using above constructed $\bm{Z}$, we construct the Laplacian matrix $\bm{Q}=\bm{I}-\bm{Z}^{t},$ where $\bm{I}$ is identity matrix and $\bm{Z}^{t}$ is the optimal Markov transition matrix. The optimal power $t$ of $\bm{Z}^{t}$ is determined by $\frac{||\bm{Y}\bm{Z}^{t}-\bm{Y}\bm{Z}^{t-1}||_{F}^{2}}{||\bm{Y}\bm{Z}^{t-1}||_{F}^{2}}\leq\epsilon,$ with $\epsilon>0$ a small constant, $\bm{Y}:=[\bm{y}_{1},\cdots,\bm{y}_{i},\cdots,\bm{y}_{M}]^{T}$ and $\bm{y}_{i}$ as described above.

Substituting $\partial V_{1}(\bm{a}^{n})/\partial\bm{a}^{n}$ and $\partial V_{2}(\bm{a}^{n})/\partial\bm{a}^{n}$ derived in (21) and (24) into (14), we obtain the full mathematical expression of MKREM algorithm proposed for PET image reconstruction

where $\bm{K}_{Ma}$ has replaced $\bm{K}_{M}$ in (14) as explained in Subsection 2.4. When $\bm{K}_{Ma}$ is constructed with $G=1$ in (15), it reduces to a single-kernel matrix $\bm{K}_{a}$ and (25) becomes the full mathematical expression of KREM algorithm (14). The pseudo-code of the proposed method is given in Algorithm 1.

The same method as that of [18] was used to generate 3D PET activity images. The 3D anatomical model from BrainWeb database consisting of 10 tissue classes was used to simulate the corresponding PET images [36],[37]. The regional time activities at 50 minutes were assigned to different brain regions to obtain the reference PET image. The matching T1-weighted 3D MR image generated from BrainWeb database was used as training image data. A slice of the anatomical model and the reference PET image are shown in Fig. 2 (a) and (b), respectively. The bright spot in the periventricular white matter is a simulated lesion created by BrainWeb.

The reference PET image was forward projected to generate noiseless sinograms using the PET system matrix generated by Fessler toolbox. The count number of noiseless sinograms was then set to 370k and 64k, respectively, to obtain the noiseless high count sinograms and the noiseless low count sinograms. The thinning Poisson process was applied to the noiseless high count and low count sinograms to generate the noisy high count and noisy low count sinograms for image reconstruction tests. Attenuation maps, mean scatters and random events were included in all reconstruction methods to obtain quantitative images. Ten noisy realizations were simulated and each was reconstructed independently for comparison.

To simulate the real scanning process and test the robustness of different methods, the BrainWeb generated MR images with noise were used as the training data $\bm{\tilde{}}{X}$ to construct the kernel matrices, dictionaries and Laplacian matrix of the proposed MKREM and KREM methods.

We used $G=2$, $\bm{K}_{l}\equiv\bm{K}_{a},\forall l$, for $\bm{K}_{Ma}$ and $\bm{K}_{l}\equiv\bm{K}_{b},\forall l$, for $\bm{K}_{Mb}$ in (15) to construct the multi-kernel matrices $\bm{K}_{Ma}=\bm{K}_{a}\bm{K}_{a}$ for the MKREM algorithm (14) and $\bm{K}_{Mb}=\bm{K}_{b}\bm{K}_{b}$ for dictionary learning. The above single-kernel matrices $\bm{K}_{a}$ and $\bm{K}_{b}$ were constructed using the Gaussian kernel (9) with $\sigma=0.5$, $J_{a}=21$ and $J_{b}=3$ in (11). The $\bm{K}_{a}$ thus constructed was used as the kernel matrix $\bm{K}$ in the KREM algorithm (14).

For kernel space dictionary learning, the training data $\bm{B}=K_{Mb}\bm{\tilde{}}{X}$ was divided into 12960 overlapping 5 $\times$ 5 patches. To ensure the generalization and efficiency of learned dictionary, 400 patches were randomly chosen from the 12960 patches. The matrix of patched training data $\mathbb{P}(\bm{B})$ was $25\times 400$. The number of dictionary atoms was set to 50 to obtain over-complete dictionaries $\bm{D}_{b}$ and $\bm{D}_{a}$. The sparsity constraint constant $s=50$ and the total iteration number = 50 were used in dictionary learning.

The data window of size $m=9$ was used to construct the Markov transition matrix $\bm{Z}$, and $\epsilon=10^{-4}$ was used to find the optimal $t=3$ for the Laplacian matrix $\bm{Q}=\bm{I}-\bm{Z}^{t}$.

For the low count image reconstruction, we tuned the weighting parameters to $\beta_{1}=0.003$ and $\beta_{2}=0.13$, and for the high count image, $\beta_{1}=0.001$ and $\beta_{2}=0.04$.

The same $\bm{K}_{a}$ as described above was used as the kernel matrix for the KEM algorithm. The post-filter of MLEM was Gaussian filter with window size 5$\times$5.

For fair comparison, the total iteration numbers and parameters of the compared methods were tuned separately to achieve their respective minimum mean square error (MMSE) in image reconstruction.

The mean square error (MSE), $\text{MSE}(\bm{x},\bm{f})=\frac{\sum_{j}^{M}\left(\bm{x}_{j}-\bm{f}_{j}\right)^{2}}{\sum_{j}^{M}(\bm{f}_{j})^{2}}$, the ensemble mean square bias, $\text{Bias}=\frac{\sum_{j}^{M}\left(\bm{x}_{j}-\bm{f}_{j}\right)}{\sum_{j}^{M}\left(\bm{f}_{j}\right)}$, and the standard deviation, $\operatorname{Var}={\frac{1}{N}\frac{\left.\sum_{i}^{O}\sum_{j}^{M}\left(\bm{x}_{j}^{i}-\frac{1}{O}\sum_{i}^{O}\bm{x}_{j}^{i}\right)\right)^{2}}{\sum_{j}^{M}\bm{f}_{j}^{2}}}$, were used for the quantitative performance evaluation and comparison of the algorithms, where $\bm{x}_{j}$ is the $j$th element of reconstructed PET image, $\bm{f}_{j}$ is the $j$th element of the reference PET image (the true image) and $O$ is the total number of noisy realizations.

Fig. 3(a) shows the true image and the images reconstructed by different algorithms for the slice of low count image. As seen from the figure, KREM and MKREM achieve lower MMSE and better visual effects than those of MLEM+F and KEM, and the KREM and MKREM images appear to be less noisy, higher contrast and clearer texture. The intensity (brightness) of the lesion area in KREM and MKREM images, especially that of MKREM images, is closer to the true image than in KEM and MLEM images. As shown in the enlarged cutouts, the edges and texture of the true image are better preserved by KREM and MKREM than by MLEM+F and KEM.

Fig. 3(b) shows the true image and the images reconstructed by different algorithms for the same slice of high count image. Due to the high count of this image, all algorithms achieve better results. But KREM and MKREM images are sharper and less noisy with lower MMSE and globally clearer image details than MLEM+F and KEM images. The lesion area in KREM and MKREM images is more pronounced than in MLEM+F and KEM images. The quantitative comparison is given in Fig. 4. As shown in the enlarged cutouts, KREM and MKREM can more effectively preserve the boundaries and details of the image than MLEM+F and KEM.

Fig. 4 plots the pixel line profiles across the lesion in the low count and high count images. The lesion area shows up as one peak indicated by the black arrow, which is significantly different from the white matter region in the high count image. Compared with other three methods, MKREM has the highest uptake in the lesion area in both low count and high count images.

The above results show that KREM and MKREM outperforms MLEM+F and KEM for both low count and high count images in terms of image visual effect and reconstruction MMSE.

Comparing the results in Fig.’s 3-4, it can be clearly seen that MKREM further outperforms KREM in image reconstruction MMSE, image visual effects and pixel uptakes for the low count and high count images.

Fig. 5 plots the average MSE (AMSE) versus the iteration number of MLEM+F, KEM, KREM and MKREM over ten noisy realizations. As seen from the plotted curves, all these algorithms approach their respective minimum AMSE’s rather quickly as the iteration number increases, but the minimum AMSE’s of KREM and MKREM are lower than those of MLEM+F and KEM, and that of MKREM is the lowest. When the iteration number further increases to beyond the minimum AMSE points, the AMSE’s of MLEM+F and KEM start to increase rapidly. As shown before in Fig. 1, such increase can last thousands iterations before the AMSE’s slowly converge to some very high values, resulting in extremely noisy and useless images. The same phenomenon has been observed and analyzed in [2] for non-kernelized MLEM in PET image reconstruction, and it is shown now for KEM in the same application. In stark contrast, the AMSE’s of KREM and MKREM are quite stable, with negligible increases from their minimum AMSE’s over thousands iterations. As shown in Fig. 5, these results hold for both low count and high count cases.

Fig. 6 shows the ensemble mean squared bias versus mean variance trade-off of different algorithms at the iterations 3 to 93 in increments of 10. As seen from the figure, the proposed KREM and MKREM always get much lower bias and variance simultaneously than those of MLEM+F and KEM for both low count and high count images.

Fig. 7 shows the AMMSE’s of KEM, KREM and MKREM versus the neighborhood size $J_{a}$ over ten noisy realizations. As seen from the figure, in both the low count and high count cases, KREM and MKREM achieve much lower AMMSE than that of KEM over a wide range of $J_{a}$ selections, and MKREM achieves the lowest. Hence, KREM and MKREM are less sensitive to the selection of $J_{a}$.