1 Introduction

Positron emission tomography (PET), is a medical imaging procedure that is helpful in the diagnosis, staging, and monitoring of several diseases. The images obtained from PET scans reflect the spatial and temporal distribution of a radiotracer within the body. Radiotracers are biochemical compounds designed to be metabolized by certain tissues of interest, such as cancer cells; moreover, radiotracers are positron-emitting particles. PET scanners are constructed to detect photon coincidences or events that are products of the annihilation resulting from positron decay.

Thus, the data obtained from a PET scan com-prise the counts of detected event coincidences, which can be presented as either sinograms or lists. The reconstruction of these data provides an image adequate for diagnostics or research. However, data produced by PET are heavily corrupted by noise resulting from several error sources, such as non-collinearity, scatter events, and death time in sensors [⁴]. Even though methods exist to alleviate and to correct most of these errors, the final quality still suffers degradation, especially at low-count rates.

Two main approaches are used to reconstruct PET data: direct methods and iterative methods. Direct methods, which were commonly used in early PET-scan procedures, provide fast recon-struction times compared with iterative methods; however, the quality of the image suffers. As the processing power of computers increased, iterative methods became more accessible and provided better reconstruction quality. One of the first iterative methods proposed was the maximum likelihood expectation maximization (EM) [¹⁸], followed by the use of ordered subsets (OSEM) [⁷], an improved version commonly used on commercially available scanners.

Recently, with the introduction of hybrid scanners such as PET/computer tomography (PET/CT) units, which are able to perform both exams on the same study, anatomical data are easily available. The literature describes several studies that take advantage of these data to increase the quality of the reconstruction. In one study [¹³], a non-local method was introduced to selectively incorporate anatomical information, making the method reliable even when information is given by fused data from another scanner.

In another study [¹¹], the authors proposed the incorporation of anatomical information through an anatomical prior with characteristics determined by the voxel intensity differences of the anatomical image.

In yet another study [⁸], a sparse image representation jointly determined by the prior anatomical image and the data from the scanner were used in reconstruction to preserve image details and smoothness. The reviewed methods incorporated anatomical information suggesting new complex priors, which can slow down the reconstruction process and cause instabilities.

In this paper, we explore the variational framework for introducing anatomical data to the reconstruction process, proposing a reconstruction method with a hybrid potential function that is applied to the gradient of the image and incorporates anatomical information as constant weights during the algorithm iterations, without adding complexity to the reconstruction process. To the best of our knowledge, such an approach has not been reported in PET reconstruction to date.

The rest of the content is organized as follows: In Section 2, the PET reconstruction is revisited. The proposed regularization scheme is discussed in Section 3, results on simulated and real images are presented in Section 4, and finally conclusions and future works are discussed in Section 5.

2 Problem Formulation

PET reconstruction aims to obtain an image given the radiotracer distribution within the body. To accomplish this, the scanner has an array of discrete detector elements that count the events generated by the annihilation of positrons.

These counts are histogramed, by angle ϕ and distance r, from the center of the scanner to obtain the object projections Pϕ(r); each projection is then stacked to form sinograms (see figure 1).

Fig. 1 Parallel projection obtained from a ring detector of a PET scanner system 

The sinogram is the radon transform represen-tation of the object scanned; thus, with enough projections of the noise-free data, the object of interest can be reconstructed accurately by finding its inverse radon transform. However, due to the physical limitations of the scanner sensors and regulations on the radiotracer dose administered to patients, the acquired data from the scanner contain heavy noise and a limited number of counts, thus requiring the application of alternative reconstruction methods. The reconstruction of PET images can be cast as a system of linear equations:

where x is a vector of size N representing the image of the object that generated the sinogram y, represented as a vector of size M, A is the system matrix, and η is noise inherent to the system. Because of the ill posed nature of (1), it is necessary to regularize the problem to enforce stability [²²].

3 Adaptive Regularization

The total variation (TV) model [¹⁷], has become one of the most commonly used methods to solve ill-posed problems by regularization.

The TV scheme achieves effective noise removal and preservation of edges; however, it suffers from undesired staircasing effects in the homogeneous regions of the image [²⁰]. The functional of the TV model for the denoising problem is given by:

where the minimum x is an estimate of the unknown image to be recovered and y is the observed data, the matrix operator H represents degradation or transformation of the data, ‖⋅‖2 is the Euclidean norm, and ‖⋅‖TV is the TV seminorm. The TV model has been adapted to perform reconstruction in [²², ²³]. In this paper, we propose the inclusion of a Tikhonov term in addition to the TV term as:

This model has been explored in several works [¹⁰, ¹⁴] for denoising, debluring, and inpainting problems. Here, the model was adapted to perform reconstruction and deal with Poisson noise. For this end, the model was derived from a Bayesian framework by first formulating the likelihood under Poisson noise for the sinogram data y and accounting for the scanner with the system matrix A:

We then proposed the prior probability for the image as a Gibbs distribution of the form:

Where λ is a parameter of the model. Once the likelihood (4) and the prior (5) are defined, the posterior probability is given by:

Note that in (6), the probability, P(y), is constant because y is the observed image. The estimated image, x^, is found using a maximum a posteriori (MAP) on (6):

using the −log function to eliminate the exponen-tials and cast the problem as a minimization, we obtain:

the last term does not have effect in the minimization, and the estimator x^ can be found as follows:

Finally substituting (5) and (4) in (9) we obtain:

Since the problem is convex, a solution can be determined using gradient descent as follows:

Here τ is a regularization parameter and N is the number of pixels in the reconstructed image. We proposed a find the parameters λ in equation (3) by using anatomical data estimated from another medical image modality, such as CT or MRI [², ¹]. To achieve this, we performed a blurred process on the anatomical image [³]; an edge map was then calculated using a labeled image. Thus, the edge map, λ was estimated from anatomical image, I, using:

where G is a Gaussian kernel with zero mean, and for the standard deviation (STD), we suggest to use a estimator of the STD of the image given by the median absolute deviation (MAD) [⁹, ⁵] defined as σ^=1.4826⋅MAD, the resulting λ is then normalized to have values in the [0…1] set. Finally this anatomical information is incorporated as weights of the regularization terms.

4 Results

This section delineates the results obtained after applied the proposed method and provides com-parisons with well-known reconstruction methods. In all experiments, the proposed method ran with 15 iterations.

In a first experiment, aimed to evaluated the resolution of the reconstructed images, we used a phantom consisting of a circular case with four rows of cylinders of 2, 3, 4, and 5 mm in diameter. Each hole was filled with activity corresponding to a background with an ​18F-fluorodeoxyglucose concentration of 1:8. The phantom described is shown in Figure 2a. The reconstruction system matrix size was of 16768×8100.

Fig. 2 Points phantom. a) original phantom, b) edge map c) EM reconstruction, and d) Proposed method 

A simulated PET scan was performed via the use of the Simset (Simulation System for Emission Tomography) software [⁶], and the resulting data was reconstructed using the EM method. Figures 2(d), and (e), shows the reconstructed images using the EM and the method proposed, respectively.

Figure 3, shows profiles across the cylinders. The proposed method gives the better image reconstruction for the different circle sizes per row, as can be seen in the curves. We also obtained quantitative comparative results in terms of Peak signal-to-noise ratio (PSNR) and structural similarity (SSIM) [²¹], between the reconstructed images and ground truth (GT) of Figure 2a, the results are shown in Table 1. Low values of PSNR and SSIM are expected, since the GT data go through a process of scanner simulation before reconstruction, which added heavy noise and blur.

Fig. 3 Profiles of the phantom of Figure 2. a) 5 mm diameter, b) 4 mm diameter c) 3 mm diameter, and d) 2 mm diameter 

In a second experiment, we aimed to test the proposed method with anatomical structures using a phantom (patient 1) from the PET Simulation Of Realistic Tridimensional Emitting Objects (PETsorteo) database [¹⁶].

Image reconstructions using the different met-hods are shown in Figure 4; this includes filtered back projection (FBP) [¹⁹, ¹⁵], EM, and the proposed method. No corrections for scattering or post filtering were applied to the images.

Fig. 4 PET Sorteo phantom patient 1. a) Edge Map b) FBP reconstruction, c) EM reconstruction, and d) Proposed method 

We evaluate the the percentage of standard deviation (%STD) [¹²] on three regions with different activity; the results (see figure 5) illustrate that the proposed method attained a lower variation and sharper images.

Fig. 5 Results of %STD evaluated on the enclosed squares. a) FBP reconstruction, b) EM reconstruction, and c) Proposed method 

5 Conclusion and Future Work

This article presents a new reconstruction method for PET data. We regularized the problem of reconstruction by using a hybrid combination of TV and Tikhonov potentials and proposed to balance both terms by the use of anatomical information incorporated through constant weights. The phantom data results show that the proposed method improves the reconstruction performance more significantly than do EM or FBP.