We developed an empirical PET model taking into account system blurring and a blind iterative reconstruction scheme that estimates both the actual image and the point spread function of the system. Reconstruction images of high quality can be acquired by using the proposed reconstruction technique for both synthetic and experimental data. In the synthetic data study, the algorithm reduces image blurring and preserves the edges without introducing extra artifacts. The localized measurement shows that the performance of the reconstruction image improved by up to 100%. In experimental data studies, the contrast and quality of reconstruction is substantially improved. The proposed method shows promise in tumor localization and quantification.

1. Introduction

In the past decade, positron emission tomography (PET) has rapidly become a popular diagnostic imaging modality for tumor detection and cancer staging. PET generates three-dimensional (3D) tomographic images of the distribution of positron-emitting radiotracers within an object, which provides quantitative functional information in vivo. However, the physics of the photon emission and detection process limits the resolution of PET images. The sensitivity and resolution of the PET system is a complex function of many factors, such as positron range, scattering, medium attenuation, photon noncollinearity, detector size, and intercrystal crosstalk. Iterative methods are often used in PET reconstruction because of the ability to handle the imaging system more precisely and to incorporate a prior information in the process.

Figure 1 illustrates a PET imaging system [1]. A positron is emitted following a nuclear transmutation. After traveling a short distance (positron range), the positron annihilates with an electron and produces two photons of 0.511 MeV almost 180∘ to each other. The photons are then detected by the detector ring. Finally, by recognizing the coincidence photons, a line of response (LOR) can be established to determine the origin of the positron emission. It should be noted that LOR is not a real line due to the finite size of the detector elements. As shown in Figure 1, photons may be scattered to travel along a totally different direction upon reaching the detector ring, and the origin of the photon pair cannot easily be described by the LOR. Therefore, reconstructions obtained from iterative algorithms that determine the system matrix solely based on LOR or extended LOR are often blurry and contaminated by artifacts.

We recently developed a blind deblurring reconstruction technique to reduce blurring in pinhole SPECT imaging [2]. It is natural to apply a similar technique to PET imaging, because SPECT and PET share a similar underlying physic. Our approach is to model the discrepancy of LOR, or the measurement, and the true emission by an independent random variable E, as demonstrated in Figure 1, and to derive the maximum likelihood (ML) solution of the new system incorporating measurement error. In this report, we describe the development of the reconstruction algorithm.

2. Blind Deblurring Reconstruction

Equation (1) illustrates the widely used Expectation-Maximization (EM) algorithm for PET reconstruction [3]:

Μ‚πœ†π‘–π‘›+1=Μ‚πœ†π‘›π‘–βˆ‘π‘—π‘π‘–π‘—ξ“π‘—π‘›π‘—π‘π‘–π‘—βˆ‘π‘–β€²Μ‚πœ†π‘›π‘–β€²π‘π‘–π‘—,(1) where pij, element (𝑖,𝑗) of system matrix 𝑃, denotes the probability of detecting an emission from voxel 𝑖,𝑖=1,…,𝑆, at detector pair 𝑗,𝑗=1,…,𝑇. The emitted photon pairs Y, radiotracer concentration Ξ›, and detected photon pairs 𝑁 are as follows:

𝑛𝑁=π‘ƒπ‘Œ,𝑗=𝑖𝑝𝑖𝑗𝑦𝑖,π‘¦π‘–ξ€·πœ†βˆΌPoisson𝑖.(2) The system matrix P can be factorized and correction factors can be applied as follows [4]:

𝑃=𝑃det.sens𝑃det.blur𝑃attn𝑃geom𝑃positron,(3) where 𝑃geom is the geometric projection matrix with each element (𝑖,𝑗) representing the probability of a photon pair produced in voxel i reaching detector pair j, ignoring attenuation, and assuming perfect photon-pair collinearity. The values in the matrix are defined by the geometries of each LOR. 𝑃attn and 𝑃det.sens are the attenuation correction and detector sensitivity normalization factors, respectively.𝑃positron, the positron range factor, is usually omitted in 18F studies in which the positron range is submillimeter, and 𝑃det.blur is the detector blurring factor used to model photon-pair noncollinearity, intercrystal scatter, and penetration, which cause the same photon detected by several adjusting detector elements. In general, the correction can be viewed as weighting and broadening the LOR on top of 𝑃geom.

In this work, instead of adding correction factors to the system matrix P, we model the uncertainty caused by positron range, photon emission angle, scatter, and detector response as a blurring factor E that introduces measurement error, which is a random variable independent from emission Y, and modify the system model accordingly. The measurement error introduces blur in the reconstruction image, and a blurred PET reconstruction can be viewed as the convolution of a low-pass point spread function (PSF) with the actual image, where both the PSF and the actual image are unknown in practice. Inspired by the blind deconvolution algorithm introduced by Holmes [5] and along the same line of research we have applied to SPECT imaging [2], we formulated a blind deblurring reconstruction algorithm. The algorithm is a modified EM algorithm that includes a convolution kernel to model the blurring factor 𝐸. The algorithm consists of two iterative updates, instead of one, to reconstruct both the object and the PSF.

In a PET imaging system, suppose π‘ π‘˜,π‘˜=1β‹―Ξ›,π‘ π‘˜βˆˆ{1⋯𝑆}, where Ξ› is the total number of emitted photons, denotes the index of location from which the π‘˜th pair of photons are emitted, and π‘‘π‘˜,π‘‘π‘˜βˆˆ{1⋯𝑇} denotes the location where the kth pair of photons are detected, as shown in Figure 1. We call these emission locations β€œtrue emission points.” A finite number of these points form an inhomogeneous Poisson random-point process having the intensity function πœ†π‘–. With the presence of measurement error, the positional measurement of each emission point is corrupted by a random translation. Let π‘’π‘˜,π‘’π‘˜βˆˆ{1⋯𝑆} denote this error vector, then the measured data for detector pixel j is related to π‘ π‘˜ and π‘’π‘˜ by

π‘ƒξ€·π‘‘π‘˜=π‘—βˆ£π‘ π‘˜+π‘’π‘˜ξ€Έ=𝑖=𝑝𝑖𝑗,𝑛𝑗=Ξ›ξ“π‘˜=1π‘ξ€·π‘ π‘˜+π‘’π‘˜ξ€Έπ‘—.(4) Here, π‘’π‘˜ is statistically independent of all π‘ π‘˜β€™s, and they are assumed to be all statistically independent of each other for all photons emitted and identically distributed with a probability density 𝑔𝑖, which is also the PSF of the PET system. In addition, the set of error vectors 𝐄={𝑒1⋯𝑒Λ} also constitutes an inhomogeneous Poisson random point process with intensity function 𝛾𝑖=Λ𝑔𝑖 [6].

Now let 𝑦𝑖 be the actual number of photon pairs emitted from voxel i, and let 𝑏𝑖 be the corresponding error vectors within voxel 𝑖; they then follow the Poisson distribution with mean πœ†π‘– and 𝛾𝑖, respectively, from the results above. We then use the EM algorithm to find the maximum likelihood solution of the system. In our application, 𝑛𝑗 and Ξ› are known measured data, whereas 𝑦𝑖 and 𝑏𝑖 are unknown data. Here we note the set of true emission vectors 𝐘 and the set of error vectors 𝐁 as

ξ€½π‘¦π˜=1,𝑦2,𝑦3ξ€Ύ,𝑏,…𝐁=1,𝑏2,𝑏3ξ€Ύ.,…(5) Then the log likelihood of I can be expressed as

𝐿(Iβˆ£π€)=βˆ’π‘–πœ†π‘–+ξ“π‘–ξ€·πœ†ln𝑖𝑦𝑖,(6) where 𝝀 is the vector notation for all πœ†π‘–. Similarly, the log likelihood of B can be written as

𝐿(𝐘∣𝐠,Ξ›)=βˆ’π‘–π›Ύπ‘–+𝑖𝛾ln𝑖𝑏𝑖,(7) where g is the vector notation for all 𝑔𝑖. The log likelihood of the complete data then can be equivalently expressed in two ways, assuming Y or B being known, that is,

𝐿(𝐈,π˜βˆ£π€,𝐠,Ξ›)=𝐿(πˆβˆ£π€,𝐠,Ξ›)+𝑙(π˜βˆ£π€,𝐠,Ξ›)=βˆ’π‘–πœ†π‘–+ξ“π‘–ξ€·πœ†ln𝑖𝑦𝑖+Ξ›ξ“π‘˜=1𝑔𝑠lnπ‘˜ξ“ξ€Έξ€Έ=βˆ’π‘–π›Ύπ‘–+𝑖𝛾ln𝑖𝑏𝑖+Ξ›ξ“π‘˜=1ξƒ©πœ†ξ€·π‘’lnπ‘˜ξ€ΈΞ›ξƒͺ.(8) By maximizing (8) using a derivation similar to those of Holmes [5] and Li [2], the following iteration can be shown to converge to the maximum likelihood estimate of πœ†π‘– and 𝑔𝑖:

Μ‚πœ†π‘–π‘›+1=Μ‚πœ†π‘›π‘–ξ“π‘—π‘›π‘—βˆ‘π‘˜π‘”π‘–βˆ’π‘–β€²β€²π‘π‘–β€²β€²π‘—βˆ‘π‘–'ξ‚€Μ‚πœ†π‘›π‘–'βˆ—π‘”π‘–'𝑝𝑖'𝑗,(9)̂𝑔𝑖𝑛+1=Μ‚π‘”π‘›π‘–Ξ›ξ“π‘—π‘›π‘—βˆ‘π‘–β€²β€²πœ†π‘–βˆ’π‘–β€²β€²π‘π‘–β€²β€²π‘—βˆ‘π‘–β€²ξ‚€πœ†π‘–'βˆ—Μ‚π‘”π‘›π‘–β€²ξ‚π‘π‘–β€²β€²π‘—,(10) where βˆ— denotes convolution. The initial πœ†0𝑖 is an image of all 1’s, and 𝑔0𝑖 is the same image normalized to 1. Equations (9) and (10) are then evaluated to acquire a new pair of estimates of πœ† and 𝑔. The PSF of the system is assumed to be real, nonnegative, band limited, and limited in extend. Letting 𝐹𝑧 be the frequency components of the PSF that are known to be zero and πΉπ‘Ÿ be the space components of the PSF that are known to be zero, the band-limited and limited-extend constraints are incorporated by executing the following steps in each iteration.

(1)The Fourier transform of ̂𝑔𝑛+1 is taken, and any frequency components that lie within 𝐹𝑧 are set to zero.(2)The inverse Fourier transform in step 1 above is taken, and any negative or complex values or values within πΉπ‘Ÿ in the spatial domain are set to zero.

The first step of the process ensures the band-limited constraint, and the second step ensures the reality, nonnegativity, and limited-extend of the PSF. Realness and nonnegativity are implicitly applied to πœ†. Equations (9) and (10) and steps (1) and (2) are then iterated until convergence occurs.

The blind deblurring reconstruction algorithm estimates both the spatial radioactivity distribution and the system PSF from the set of blurred projection images. The iteration for reconstruction can be understood as replacing the forward projector in the original EM (denominator of (9) with the new projector using the convolved radioactivity map, and the iteration for solving the PSF can be understood as blind deblurring. This iteration differs from the general image-blind deconvolution in the sense that the kernel is partly known: 𝑝𝑖𝑗, the system matrix, is in fact part of the blurring kernel. The more precise the model 𝑝𝑖𝑗 is, the closer the remainder of blurring kernel, or g, is to a true delta function. In addition, instead of deconvolving an image where both the input and output are 2D images, the input of blind deblurring reconstruction is a series of projection images, and the output is a 3D-image array. This property gives us much more knowledge of the noise distribution within the object, because instead of a single-shot image, we now essentially have multiple samples for each point in the 3D array (although mixed with other points). Both simulation and experimental data were used to validate and evaluate the performance of the blind deblurring EM (BDEM) technique.

3. Methods

3.1. Monte Carlo Simulation

PET simulations were performed utilizing the NCAT phantom covering the chest region [7, 8]. The PET emission data were generated using the Monte Carlo method, with 1 million counts per slice using a geometry corresponding to the design of the GE discovery DSTE PET/CT scanner (GE Healthcare, Milwaukee, WI, USA). The images were reconstructed using the standard EM and the BDEM. The reconstruction image sets were then evaluated using visual inspections and line profiles. The EM reconstruction image was postsmoothed using a Gaussian kernel of 1 pixel full width at half medium, and no post processing was done for BDEM.

3.2. Convergence Study

One major concern with regard to the type of blind deblurring technique used is the effectiveness and convergence of the algorithm. Two hundred iterations of EM and BDEM reconstructions of the NCAT phantom were performed to evaluate the convergence property of the BDEM algorithm. In each iteration, the measurement log-likelihood L of the reconstructed image πœ† for EM was calculated as follows [6]:

𝐿(π€βˆ£π‘)=𝑗𝑛𝑗logπ‘–π‘π‘–π‘—πœ†π‘–βˆ’ξ“π‘–π‘π‘–π‘—πœ†π‘–βˆ’log𝑛𝑗!ξƒͺ(11) whereas for BDEM, the log-likelihood L of the reconstruction image was calculated as

=𝐿(π€βˆ£π‘)𝑗𝑛𝑗logπ‘–π‘π‘–π‘—ξ€·πœ†π‘–βˆ—π‘”π‘–ξ€Έβˆ’ξ“π‘–π‘π‘–π‘—ξ€·πœ†π‘–βˆ—π‘”π‘–ξ€Έβˆ’log𝑛𝑗!ξƒͺ.(12) Two different settings of limited-extend constraint were used for BDEM reconstruction, and the difference in log-likelihood was compared.

3.3. Comparison with Image Deconvolution

The blurred reconstruction image could always be modeled as a convolution of the true radioactivity with a blurring kernel β„Žπ‘–:

𝑓𝑖=πœ†π‘–βˆ—β„Žπ‘–,(13) where 𝑓𝑖 is the standard EM reconstruction image and * denotes 2D linear convolution. Assuming β„Žπ‘– can be estimated, one could first compute 𝑓𝑖 using standard EM iterations (with Μ‚πœ†π‘– being replaced by 𝑓𝑖) and then deconvolve 𝑓𝑖 with β„Žπ‘–. However, the Μ‚πœ†π‘– so obtained is not the maximum likelihood estimate of πœ†π‘– given the measurement error; also, the kernel β„Žπ‘– is generally a complex unknown function and is hard to measure. We used the Wiener filter and image blind deconvolution algorithm to denoise and deconvolve the EM reconstruction image, and compared the results with the BDEM reconstruction.

3.4. Physical Phantom Study

A physical phantom study with a water-tank phantom containing 4 spheres of diameters 1.5 cm, 2.0 cm, 3.0 cm, and 3.0 cm was conducted to assess the performance of a BDEM algorithm for a real scanner. Both the water tank and the spheres were filled with 18FDG, with activity concentrations of 0.20πœ‡Ci/cc and 1.01 0.20πœ‡Ci/cc, respectively. The tank was imaged with a GE Discovery DSTE PET/CT scanner (GE Healthcare) in 2D list-mode; the sinogram was extracted and reconstructed using the EM and BDEM techniques.

4. Results and Discussion

4.1. Monte Carlo Simulation

Figure 2 compares the image slices through a 2-cm-diameter tumor in the right lung produced by the EM and BDEM reconstruction techniques. It is clear that the image quality of the BDEM reconstruction is superior to that of the EM reconstruction. The contrast of the lesion of interest is greatly improved, edges are preserved, and artifacts are suppressed. The transverse-view line profile across the tumor also confirms the improvement. The peak contrast of the tumor in the BDEM reconstruction is double that of the EM reconstruction. Although no postsmoothing was performed for the BDEM reconstruction, the noise level in the background of the BDEM reconstruction is lower than that of the EM reconstruction, with postsmoothing by a 1-pixel full width at half medium Gaussian kernel. The mean and standard deviation of a small region of interest in the background region are 36.20 and 25.12 (resp.) for the Gaussian-smoothed EM and 34.83 and 21.28 (resp.) for BDEM, with no postprocessing.

4.2. Convergence Study

Figure 2 shows the log-likelihood as a function of iteration numbers as calculated in (11) and (12). The measurement log-likelihood of BDEM is a monotonic increasing function of iteration number, which confirms the convergence of the algorithm. It also varies with different PSF constraints and lies within the envelope defined by the EM algorithm, as demonstrated in the figure. Therefore, setting a proper constraint on PSF is important. If the initial PSF is a delta function and no other constraint is set, the BDEM would regress to the regular EM solution. Figure 2 also shows the log-likelihood of BDEM with a 10Γ—10- and a 20Γ—20-PSF spatial constraint being applied for a 128Γ—128 reconstruction. The small discrepancy in the convergence of two different constraints indicates that BDEM is not a strong function of PSF constraint: as long as the constraint is reasonable, BDEM would converge to very similar solutions. The figure also indicates that BDEM converges a little faster than conventional EM in terms of number of iterations. It should be noted, however, that the amount of computation for each BDEM iteration is about twice that of a computation of an EM iteration. We have not tested this premise, because the algorithm converges reasonable fast even without order-subset, but the principle of order-subset EM (OSEM) should be applied if necessary.

4.3. Comparison with Image Deconvolution

The results obtained by deblurring the EM reconstruction image with Wiener filter and the image-blind deblurring technique were displayed and compared with the results of the EM and BDEM reconstructions (Figure 4). These image deconvolution results are clearly inferior compared to BDEM reconstruction. The images are still noisy, and the visual improvement of deconvolution is limited. Our explanation is that the BDEM algorithm utilizes the statistical information of both emission and measurement noise in the projection data, which is loss in the reconstruction image. The noise in the projection data commonly appears as streaks or other local or global artifacts in the reconstruction image, which makes it much harder to either identify or clear just from the reconstruction image.

4.4. Physical Phantom Study

Figure 4 shows one slice of water phantom with spheres. As with the Monte Carlo study, the image contrast of the BDEM has been improved over the conventional EM, and the mass and edges are well preserved in the reconstruction. The image blurring is reduced, and the contrast is greatly improved (up to 50%), as observed from the line profile from the transverse slice. The shape of the lesions is not distorted, indicating that the algorithm preserves the general shape of objects.

5. Conclusion

In this work, we demonstrated highly desired reconstruction results with no complex assumption about the imaging system or the object. The blind deblurring reconstruction technique can significantly improve the quality and contrast of the reconstruction as demonstrated in both simulation and experimental scans. This algorithm does not only reconstruct the radiotracer map, but also determines the PSF of the system. The masses and edges are well preserved in the reconstruction image, which can be extremely useful when doctors need to localize, segment, or tally the activities in the possible tumor. Future studies would involve the application of this system in patient imaging and quantitative studies.