A Novel Kernel-Based Regularization Technique for PET Image Reconstruction

Abdelwahhab Boudjelal 1,2,*, Zoubeida Messali 3 and Abderrahim Elmoataz 2 1 Electronics Department, University of Mohammed Boudiaf, 28000 M’sila, Algeria 2 Image Team, GREYC Laboratory, University of Caen Normandy, 14050 Caen CEDEX, France; abderrahim.elmoataz@unicaen.fr 3 Electronics Department, University of Mohamed El Bachir El Ibrahimi, 34030 Bordj Bou Arréridj, Algeria; messalizoubeida@yahoo.fr * Correspondence: abdelwahhab.boudjelal@unicaen.fr


Introduction
Efforts are underway to improve the image reconstruction quality achieved by positron emission tomography (PET) [1].The purpose of PET is to provide 3D imaging of the administered tracer distribution, which is believed to be associated with the distribution of activity of certain molecular mechanisms underlying pathology or disease models [2].More detailed PET scans require the injection of a small amount of biologically-relevant material, such as glucose or oxygen, that has marked radionuclides [3].All of the isotopes used are radioactive with a rapid decomposition time by positron emission.The most common isotope employed in PET is fluorine-18 [4].This tracer is very useful because of its long half-life (110 min) and because it disintegrates by emitting positrons having a low positron energy that contributes to high-resolution imaging acquisition.The molecule accumulates in cells that are most metabolically active, or if the molecule is receptor-specific, it collects in cells where the receptors are present.This evaluation of the cellular and physiological function by nuclear medicine can be effectively used to locate and determine the extent of a disease in the body.Nowadays, a major set of PET applications is in oncology [5], cardiac imaging [6], and brain imaging [7].With the developments of specific radiotracers, several applications of PET in neurology have been found during the last few years, such as dopamine neurotransmitter imaging in Parkinson's disease and tau imaging [8] in Alzheimer's disease.PET works by a radioactive tracer that has been introduced into the body, decaying to emit a positron β+, which can only travel a few millimeters (the positron range [9]), before colliding with an electron.The antimatter-matter collision annihilates both particles and results in a pair of 511-keV gamma photons that travel in nearly opposite directions.A pair of detectors is connected by a line of response (LOR) and is tuned to pick up photons within a coincidence-timing window.
The annihilation events can be stored in list-mode format or in a 3D matrix according to the position of their corresponding line of response.The process of acquiring the projection data is called emission scan.The histogram data are said to lie in a sinogram or projection domain, and the object to be scanned is said to lie in an image domain.The emission tomography uses the projection data to reconstruct the spatial distribution of the radioisotope by taking into account the geometric factors, noise properties, and physical effects.To estimate the radioactivity distribution accurately, we must consider the effects of the attenuation that result from the photoelectric absorption effects as a consequence of the interaction of the emitted gamma photons through the subject matter and the gantry.
The PET image is reconstructed from the acquired projection data, but as the count of the projection data declines, there is an increase in noise in the reconstructed PET image.The most simple method used to reconstruct the PET image is the analytical method called filtered back-projection (FBP) [10].Yet, because of the Poisson-distributed noise that is characteristic of projection data, FBP-reconstructed PET images are often of low quality.Alongside regular reconstruction methods, noise smoothing can be performed in the pre-reconstruction phase of the sinogram domain [11][12][13][14][15], in the post-reconstruction phase of the image [16][17][18], or during the iterative process of statistically-based reconstruction [19].An alternative method that incorporates prior knowledge into de-noising images is non-local means (NLM) [20,21].In general, applying denoising methods correctly and accurately post-reconstruction is more difficult than within the reconstruction.
This paper uses the maximum-likelihood expectation-maximization (MLEM) algorithm to reconstruct PET images [22,23].Using MLEM for PET image reconstruction that has a noisy dataset results in a noisy image.Therefore, to improve the quality of the image, we propose a novel kernel-based image regularization technique.Of the various possible options of the kernel function, an exponentially-modified Gaussian kernel was used [24][25][26]; this is comparable to conventional MLEM reconstruction in terms of simplicity.Noise is modeled in the projection domain, where independent Poisson random variables effectively model the PET data.The results indicated that the proposed kernel-based MLEM (κ-MLEM) method provides effective reconstruction and resolution results that exceed those generated by conventional MLEM and advanced pre-reconstruction methods.
In Section 2 of this paper, the novel MLEM kernel-based regularization image reconstruction method is outlined.Section 3 presents the PET setup and the evaluation criteria for the quality of the reconstruction.Results are presented in Section 4, together with a comparison of the κ-MLEM and other methods.Section 5 is the conclusion.

Data Acquisitions and Sinograms
PET scans detect vast numbers of pairs of annihilation photons.The total number of pairs of photons measured by a particular detector pair is proportional to the integrated radioactivity that occurs along the line joining the detector pair.LORs pass through a single point in the patient body and trace a sinusoid curve in the raw data histogram; hence the term sinogram for the raw data format.A simple example of an object in the spatial domain and its corresponding sinogram are illustrated in Figure 1.

PET Data Corrections
Collecting PET data is not a perfect process.The photons emitted within the patient are attenuated by photoelectric interactions; the efficiency of detector elements is variable, and together with recording true coincidence events, random and scattered coincidences are accrued.To create a clinically-valuable image and reliable quantitative information, these artifacts and errors need to be corrected.
Attenuation correction (AC) is the most important data adjustment [27].The effect of encountering dense or a greater quantity of material on the path between the annihilation site and the detector leads to photon absorption or scattering.The result is an attenuation of photons in those environments compared to photons that do not encounter such material.Images that are reconstructed from projection data of photons in less dense areas of the body (e.g., the lungs) appear dark, as they have emitted more photons than the surrounding tissue (e.g., the mediastinum), which is dense.
Random coincidences are also corrected [28].These need to be estimated and minimized from the collected measurements in each LOR to find the sum of scattered and true coincidences, which is necessary for quantitative measurement.

Maximum-Likelihood Expectation-Maximization Algorithm
Following the acquisition and correction of PET data in the sinograms, an estimate of the in vivo tracer distribution needs to be determined.This is a mathematically complex stage, and is covered in detail in [29][30][31].It is assumed that the object being imaged has a radionuclide distribution represented by the matrix x j , where the index j denotes the location of an individual volume element ("voxel") in the object.
The purpose of the tomographic reconstruction algorithm is to calculate the radionuclide image x j from the set of counts b i recorded by the external detector array circling the patient, where the index i identifies the position of each detector element used in making our measurement from the patient.In this mathematical derivation, the radionuclide image x j is represented as a matrix in which each matrix element x j represents the radionuclide concentration at a point identified by the index j.Similarly, the recorded radionuclide projection data b i are represented as a matrix, in which each matrix element b i identifies both the angle and the location of the detector, where the data are recorded.
We assume that we know the probability P ij of a photon emitted in voxel j being detected by detector i.The value of the detection probability P ij is directed by physical phenomena, such as the geometric efficiency of the collimator, the attenuation provided by the material surrounding the radioactive voxel j, and the detection efficiency of the imaging system.
The behavior of the imaging system can be represented by using a matrix equation that couples the radionuclide concentration x j in the object to the complete set of detector measurements b i using the known value of the probability elements P ij : As the values of b i can be measured with the detector array and since the value of the probability elements P ij can be calculated from the physics of the imaging system, the matrix equation Equation ( 1) can be mathematically inverted to solve for the radionuclide distribution x j .In reference to matrix Equation ( 1), the radionuclide concentration x j is unknown, but a solution can be assumed.Consequently, x j is the best estimate of the true radionuclide distribution.
Since the probability matrix P ij is known, it is possible to estimate the detector counts (b i ) that would be obtained for the estimated radionuclide distribution x j using the forward projection operation: Reconstruction is primarily founded on a mathematical model that relates parameters to the observed and unobserved data.The mathematical model is based on the notion that emissions occur according to a spatial Poisson point process inside the region of interest in the source.For each detector position i, the appropriate model is a Poisson-distributed random variable b i [22].The probability matrix P ij is known from the detector array geometry, the attenuation provided by the material system and the detection efficiency of the imaging system.Suppose y i is the expected values of Poisson-distributed random variables b i where: The expectation of the detected data is combined with the model parameters: Since the variables λ ij -the number of emissions in voxel j detected in detector i-are independent Poisson-distributed random variables with the mean, the likelihood function is then: Although the sum in ( 6) is complicated, it shows that the log-likelihood is concave as a function of x: The E-step is computed as: Since we know that λ ij are Poisson-distributed random variables with mean x ij , after some manipulations: In the M-step, we maximize (10) Combining the E and the M steps: We start the algorithm with an initial guess x [0] j , and then, in each iteration, if x [k] j denotes the current estimate of x j , the new estimate is defined as: The MLEM algorithm is summarized as: 1.
Start with an initial estimate x [0] If x [k] j denotes the estimate of x j at the k th iteration, define a new estimate x [k+1] j by using Equation (11), 3.
If the required accuracy for the numerical convergence has been achieved, then stop.

Regularization
The expectation maximization algorithm is ill-conditioned, which means that enforcing the maximization of the likelihood function may result in a solution with severe oscillations and noisy behavior.

The Non-Local Means Filter
The NLM algorithm is a nonlinear spatial filter that uses a weighted average of the intensity of adjacent pixels to adjust the intensity value of each pixel [20,21].While comparisons of neighboring pixels can be performed between pixels separated in the image by any distance, for the benefit of computation speed, comparisons are typically limited to a narrow search window around the pixel undergoing value adjustment.
Mathematically, this method can be described as: where where P is the patch size around each pixel that is being compared, K G is a Gaussian kernel used to reduce the weight in the similarity calculation given to patch pixels further from the center, and h adjusts the intensity similarity-to-weight relationship.

The Anisotropic filtering
The anisotropic filter (AF) attempts to avoid the blurring effect of the Gaussian by convolving the image x at j only in the direction orthogonal to Dx(j).The idea of such a filter goes back to [32].It is defined by: for j such that Dx(j) = 0 and where (j, l) ⊥ = (−l, j) and G h is the one-dimensional Gaussian function with variance h 2 .The image method noise of an anisotropic filter AF h [20] is , where h is a filtering parameter (which usually depends on the standard deviation of the noise), and curv(x)(j) denotes the curvature, the signed inverse of the radius of curvature of the level line passing by j.This method noise is zero wherever x behaves locally like a straight line and large in curved edges or texture.Consequently, the straight edges are well restored, while flat and textured regions are degraded.
An extended discussion of the anisotropic filter for attenuation correction in PET can be found in [11].

Proposed Kernel-Based Exponentially-Modified Gaussian Regularization
Inspired by the kernel-based image representation for PET images [25], we propose a new kernel-based image regularization technique to improve the PET image reconstruction.We define x j as a feature map for pixel j that is the output of the kernel, which is a mixture of two sets of information.The first set is the input image, which has a matrix of pixels, where a pixel consists of an integer value between zero and 80.The second set is the convolution kernel, which consists of a single matrix of floating point numbers.
Convolution can be considered a mixing of information.This is achieved by taking an image box from the input image that is the size of the kernel.This experiment uses a 256 × 256 image and a 7 × 7 kernel, so it would take 7 × 7 boxes.An element-wise multiplication with the image box and convolution kernel is performed, and the sum of this multiplication results in a single pixel of the feature map, as depicted in Figure 2, by using a 3 × 3 kernel.The image intensity at pixel j is represented by a linear combination of the kernels,

Kernel Input Output
where N is the total number of pixels in the image and χ is the kernel coefficient image.κ(•, • ) is a kernel function encoding image prior information.
x = Kχ In Equation ( 16), each column of K can be considered as a basis function for image representation.We proposed to use kernels including the exponentially-modified Gaussian [33]: where erfc is the complementary error function, defined as x e −t 2 dt, and σ, µ, and δ are the kernel parameters in the respective kernels.
In most instances, the kernel in the image is center-originated, indicating that the center point of a kernel is κ(0, 0).For example, if the kernel size is five, then the array index of five elements will be −2, −1, 0, 1, and 2. The origin is located at the middle of the kernel.
Because PET projection data b i are well modeled as independent Poisson random variables with the log-likelihood function [34], the expectation E[b i ] can be related to the unknown emission image x j through: with P being the detection probability matrix and r the expectation of random and scattered events.Substituting Equation ( 16) into (18), we obtain the following kernel-based forward projection model: Combining the kernel-based projection model and the standard MLEM algorithm [22] given in Equation (11), it can be directly applied to find the ML estimate of χ because of (20).The resulting kernelized EM update of χ at iteration (k + 1) is:

Simulated PET Data (Setup)
To assess the overall performance of the proposed kernel regularization algorithm, we simulated a PET scanner in a 2D mode, which has 420 crystals in each ring with each crystal having a cross-section of 6.3 × 6.3 mm 2 .A Hoffman brain phantom (Figure 3) [35] was used to simulate the radioactivity distribution in a single slice of PET tracer through grey matter, white matter, and three tumors.The digital design of the tumors in the phantom provides essential information about the reconstructed images.To create a projection data, the simulated phantom image was forward projected using the system matrix to generate the noise-free projection data.Once the simulated noise-free sinograms were produced, a 30% uniform background was added to simulate mean randoms and scatters.The fraction

Results and Discussion
In this study, the experimental kernel-based MLEM (κ-MLEM) algorithm was compared against conventional MLEM [22], the penalized likelihood reconstruction regularized by anisotropic diffusion filter (ADF) [17], the MLEM precondition reconstruction algorithm that was regularized by non-local means (pre-NLM) and post-reconstruction de-noising methods applied post-MLEM (post-NLM) [20].One hundred and fifty iterations were applied to the MLEM algorithm used in each of the regularization-based methods.Apart from the conventional MLEM algorithm, spatial smoothing was used to enhance the quality of the image.All images were sized at 256 × 256 pixels.The box and neighborhood size were 7 × 7.
We applied our technique to the Hoffman brain phantom, reconstructing the phantom erroneous projections.Figure 4 shows the reconstructed images.As is shown, the visual quality of the reconstructed image of the phantom using the κ-MLEM algorithm is comparable to the other methods.Indeed, compared to some, the experimental algorithm preserves edges better.There is less noisy projection data, and the reconstruction is more accurate.The method was able to detect the three tumors of different sizes, the smallest of which was not picked up by ADF and pre-NLM and only faint in the post-NLM and conventional MLEM images.Sensitive and reliable visualization of brain tumors is important for medical diagnosis and prognosis.
The effectiveness of noise removal for the test algorithm was comparable to that of the post-NLM method; however, the intensity in the region of interest was appreciably higher in the latter.Compared to the other methods, the κ-MLEM method generates a superior intensity profile, whilst preserving the edges.
The SNR, MSE, PSNR, and NCC quality measurements of the reconstructed images derived from the proposed algorithm and the other four reconstruction algorithms for the tumor regions of interest (ROIs) of the Hoffman brain phantom selected with a rectangle in (Figure 3) are shown in Figure 5.The minimum MSE achieved by all of the different methods for all iterations are shown Figure 5a.The peak signal-to-noise ratio (PSNR) for each iteration is presented in Figure 5b.The normalized cross-correlation (NCC) and the signal-to-noise ratio values are shown in Figure 5c,d, respectively.These data indicate that the measurements provided by κ-MLEM are of a higher quality than the other methods tested.Compared to the performance parameters of the other regularization methods, those of the proposed method were enhanced-especially at the lower number of iterations-and remained almost constant after 100 iterations; for the conventional ADF-MLEM and pre-NLM, however, the performance parameters began to decrease after 80 iterations.The data presented in these images indicate that the κ-MLEM-generated images have the highest PSNR, SNR, and NCC at all iterations and have a lower MSE; thus, the method outperforms the other methods, and this effect is heightened as the iterations increase.We can see the post-NLM algorithm is clearly demonstrating a better performance than all of the other alternative techniques and a performance quite close to the κ-MLEM method.Indeed, the quality measurements of conventional MLEM at 150 iterations are not dissimilar to those of κ-MLEM at 50 iterations (Figure 6).The lowest threshold number of iterations the technique needs to create an acceptable reconstructed image is 20, making it fast as well as efficient.κ-MLEM is also effective in removing star artifacts that are commonly generated by the conventional MLEM algorithm.The 1D line profile is the horizontal line that crosses the image in the two tumors, as shown in Figure 7 for the middle tumor and in Figure 8 for the small tumor.Noisy images in uniform regions are shown as spikes, as indicated by the conventional MLEM and post-NLM images; on the other hand, κ-MLEM line profiles are closer to being noise-free.The experimental method also nullifies aerial pixels, redistributing their values to pixels within the phantom.κ-MLEM reconstructed the ideal profile more effectively than the other methods, and the image produced was a close approximation of the original phantom image.To evaluate the performance of our proposed algorithm in the detection of tumors, we compare the size of tumors of the original phantom and tumors of reconstructed images from different algorithms.The reconstruction outcomes of the simulated phantom yielded the two rectangular areas in Figure 9, which highlight the two smaller tumors.Subsequently, we used a threshold equal to 75% to segment the tumors in PET images.The proposed algorithm successfully preserved the edges of tumors.For example, the middle tumor edges were effectively preserved when applying the proposed algorithm, whereas artifacts and deviations were more likely to occur with the NLM-MLEM and conventional MLEM algorithm and appear distorted with the ADF-MLEM algorithm, as shown in Figure 9. Furthermore, the intensity distribution within the middle tumor was made more homogeneous by the proposed algorithm as compared to the conventional MLEM algorithm.In addition to this, the proposed algorithm could enable the detection of the small tumor with the same size of the original tumor, unlike other algorithms, which did not show tumors at all.

Conclusions
The new reconstruction algorithm presented in this paper aims to improve the noise in PET images.Based on conventional MLEM, the experimental algorithm uses a box-kernel of the exponentially-modified Gaussian.The kernelized image model can incorporate prior information in the forward projection model κ-MLEM.The achievements of κ-MLEM were quantitatively superior, and there was a significant reduction in noise at matched contrast when tested against conventional MLEM and NLM regularization-based methods in PET simulations.Promising results were also achieved when the algorithm was used to reconstruct a Hoffman phantom brain with tumors of varying sizes.We believe that our method can make an important contribution to the diagnosis and prognosis of brain tumors and has an application in clinical, education, and medical research.

Figure 1 .
Figure 1.A simple object and the sinogram.
y[k] is the new pixel value at coordinate k, C N is a normalization constant, SW is the search window centered at k, x[k − n] is the original pixel value at the (k − n) coordinate within the search window, and W[k − n] is the weight given the intensity value of x[k − n].The weights are calculated as:

Figure 2 .
Figure 2. Kernel feature extraction, calculating the box-wise multiplication with the image box and convolution kernel.

Figure 5 .Figure 6 .
Figure 5. Comparative analysis of reconstructing a Hoffman phantom using various reconstruction methods by varying iteration number.MSE: mean square error; NCC: normalized cross-correlation; PSNR: peak signal-to-noise ratio.

Figure 7 .
Figure 7. 1D profiles corresponding to reconstructed images for the middle tumor.

Figure 8 .
Figure 8. 1D profiles corresponding to reconstructed images for the small tumor.

Figure 9 .
Figure 9. Zoom-in of the tumors of the Hoffman phantom reconstruction achieved with various algorithms.