Evaluation of the Feasibility and Quantitative Accuracy of a Generalized Scatter 2D PET Reconstruction Method

Abstract

Scatter degrades the contrast and quantitative accuracy of positron emission tomography (PET) images, and most methods for estimating and correcting scattered coincidences in PET subtract scattered events from the measured data. Compton scattering kinematics can be used to map out the locus of possible scattering locations. These curved lines (2D) or surfaces (3D), which connect the coincidence detectors, encompass the surface (2D) or volume (3D) where the decay occurs. In the limiting case where the scattering angle approaches zero, the scattered coincidence approaches the true coincidence. Therefore, both true and scattered coincidences can be considered similarly in a generalized scatter maximum-likelihood expectation-maximization reconstruction algorithm. The proposed method was tested using list-mode data obtained from a GATE simulation of a Jaszczak-type phantom. For scatter fractions from 10% to 60%, this approach reduces noise and improves the contrast recovery coefficients by 0.5–3.0% compared with reconstructions using true coincidences and by 3.0–24.5% with conventional reconstruction methods. The results demonstrate that this algorithm is capable of producing images entirely from scattered photons, eliminates the need for scatter corrections, increases image contrast, and reduces noise. This could be used to improve diagnostic quality and/or to reduce patient dose and radiopharmaceutical cost.

1. Introduction

Compton scattering degrades image contrast and compromises quantitative accuracy in positron emission tomography (PET) [1, 2]. Scattered coincidences are typically considered as noise which reduces PET image quality. This issue is more serious when operating in 3D mode without slice-defining septa and in large patients, where the scatter fraction can be as high as 40–60% [3–6]. Consequently a number of approaches for estimating and correcting scattered coincidences in PET data have been proposed [3, 7–19]. Most of these techniques estimate a scatter sinogram, which is used to subtract the scatter from the projection data [20] in precorrection methods [21] or as a constant additive term incorporated in a reconstruction algorithm [22–24]. Inaccuracy in the estimation of the scatter sinogram will introduce significant biases in the activity distribution [6]. The subtraction-based correction methods destroy the Poisson nature of the data, reduce the system’s sensitivity, and amplify image noise [3, 25].

With list-mode acquisitions in modern PET and improved detector technology, the use of the energy of individual photons becomes feasible and some authors have proposed new approaches that include the energy information in the estimation of the scatter distribution and image reconstruction process [6, 26–28]. These approaches attempt to improve the accuracy of the rejection of scattered coincidences from measured data but are limited by the energy resolution of the detectors [29].

Since the energy of detected photons carries some probabilistic information about the spatial distribution of the annihilation and by taking advantage of Compton scattering, scattered coincidences are also a potential source of latent information and can be included in PET image reconstruction. Incorporating scattered coincidences directly into the reconstruction algorithm eliminates the need for scatter correction and could improve both image quality and system sensitivity since a lower energy threshold can be used. In this work, we propose a reconstruction method similar to that of Conti et al. [30] that directly includes scattered photons in the image reconstruction algorithm instead of correcting for them as do conventional emission imaging methods.

An important difference between our work and that of Conti et al. [30] is that the method proposed by Conti et al. makes use of both time-of-flight (TOF) information and the energy of individual photons to reconstruct scattered and unscattered PET coincidences [31]. This paper focuses on a non-TOF PET algorithm which does not make use of the time difference between the two detected annihilation photons. Another difference is that we use Compton kinematics to predict the locus of possible scattering positions to confine the annihilation position instead of using a pixel-driven approach, which can reduce the computational workload. This work has been presented, in part, at the World Congress on Medical Physics and Biomedical Engineering in 2012 [32].

In this paper, we evaluate the feasibility of this approach, the effects of the scatter fraction, and the choice of a new coincidence selection threshold, which uses both the scattered photon energy and detector position for each coincidence, on the image quality obtained using a generalized PET reconstruction algorithm.

2. Materials and Methods

2.1. Reconstruction Theory

In clinical PET systems, more than 99.7% of the scattered events undergo Compton interactions in water at 511 keV [33]. Coherent scattering is neglected because its contribution to the total cross-section is small for the energies of interest in nuclear medicine [25]. In approximately 80% of detected scattered events, only one of the two annihilation photons undergoes a single Compton interaction [5, 8, 34]. This assumption has been validated with both Monte Carlo (MC) simulation and experiment measurements in the single-scatter simulation (SSS) technique [35, 36].

Figure 1 illustrates a scatter coincidence in a patient. The two annihilation photons are generated at source X and the unscattered photon is detected at A while the second photon undergoes a scatter at S and is detected at B with energy 𝐸^′. The scattering angle 𝜃 is given by the Compton equation 𝑚cos𝜃=1+^𝑒𝑐²𝐸−𝑚_𝑒𝑐²𝐸^′, (1) where the scattering angle 𝜃 is defined as the angle between the direction of the incident and the scattered photon, 𝑚_𝑒𝑐² is the electron rest mass energy equal to 511 keV, 𝐸 is the incident photon energy, and 𝐸^′ is the scattered photon energy in keV. The energy of the scattered photon decreases as the scattering angle increases with the energy of the scattered photon ranging from ~170 to 511 keV, corresponding to backscatter and forward scatter, respectively. By taking advantages of the kinematics of Compton scattering, a 2D surface described by two circular arcs (TCA) connecting the coincidence detectors (rather than the conventional straight line assumption which is true only for unscattered photons) which scribes the possible scattering locus and encompasses the annihilation position can be identified as shown in Figure 1. The locus defined above ignores the positron range, and the accuracy is closely related to the energy resolution of the detectors. Our preliminary investigation assumes that the detectors have perfect energy resolution and data was simulated with the Monte Carlo code.

Figure 1 

A diagram of a Compton scattering event in a patient. The two antiparallel photons are generated at the annihilation position X shown by the green dot. The unscattered photon is observed by detector A while the other photon suffers a Compton scattering event at S, deviates from its initial path by an angle 𝜃, and is detected by detector B. This figure also illustrates the two circular arcs (TCA), shown as blue dotted curves, which describe the locus of possible scattering locations for the scattered photon and which enclose the annihilation position.

The size and shape of the area encompassed by the TCA are a function of the scattering angle and the detected positions of the coincidence photons. Figure 2 illustrates that, in the limit, where the scattering angle approaches zero and the energy of the scattered photon approaches 511 keV, the shape approaches the line of response (LOR) for unscattered photons. Thus, the true coincidences may be considered as a subset of the scattered coincidences with zero scattering angle. We therefore propose a generalized scatter (GS) approach which generalizes the conventional meaning of scatter to include both true and scattered events and will show in Section 2.2 that they can be used to reconstruct an image in a similar way.

Figure 2 

The TCA shapes versus the scattering angles for the same two detectors. Each TCA consists of two symmetrically circular arcs. The inner to the outer arcs correspond to scattering angles of 0, 30, 60, 90, 120, and 150 degrees, respectively. When the scattering angle is smaller than 90 degrees, the TCA is made up of two minor arcs; when the scattering angle is 90 degrees, the TCA is a circle; when the scattering angle is larger than 90 degrees, the TCA is made up of two major arcs. As the area encompassed by the TCA increases with scattering angle, the ability to accurately determine the annihilation position decreases.

The work in this paper is carried out in two dimensions, although the same approach can in general be carried out in three dimensions where it will ultimately be of greater value.

2.2. The Generalized Scatter Algorithm

𝑃_ab(𝜃) represents the mean number of coincidences, where an unscattered photon is observed at A while the other annihilation photon undergoes a Compton scattering through an angle 𝜃 and is observed at B with an energy 𝐸^′ (see Figure 1).  𝑃_ab(𝜃) can be calculated using the following equation: 𝑃_ab(𝜃)=𝜏^B_A^S_A𝑓_𝑥2𝜋𝑑𝑥∗𝜌_𝑒(S)∗𝑑𝜎^KN_𝐶∗1𝑑𝜃2𝜋∗𝑒^{−(∫SA𝜇𝑑𝑙+∫BS𝜇𝑑𝑙)}𝑑𝑠. (2) The inner integral in the above expression calculates the total photon flux that can reach one of the possible scattering points (say S) due to the source 𝑓_𝑥 lying on the line segment AS. If the effect of the noncollinearity of the annihilation photons is neglected, the photons are emitted isotropically and no distinction can be made between the two annihilation photons. Thus, the probability of the annihilation photons traveling along AS is 1/2𝜋 instead of 1/4𝜋. Once these photons arrive at one of the possible scattering positions S, the possibility of photons suffering a Compton interaction with a scattering angle 𝜃 is linearly proportional to the differential Klein-Nishina electronic cross-section  𝑑𝜎^KN_𝐶/𝑑𝜃 and the electron density 𝜌_𝑒 at this point. The second 1/2𝜋 takes into account that only a portion of the photons scattered into the solid angle around 𝜃 are detected at B. The outer integral sums the scattered coincidences over all the possible scattering positions along the TCA.  Tau (𝜏) is the acquisition time, and μ is the linear attenuation coefficient.

To reduce the computational load, the electron density was assumed to be uniformly distributed, and for this geometry, the attenuation was taken to be negligibly small. Thus, 𝑃_ab(𝜃) is approximately proportional to the integral of the source intensity within the area confined by the TCA and 𝑃_ab(𝜃) can be written as 𝑃_ab(𝜃)=𝐶∗𝑑𝜎^KN_𝐶𝑑𝜃^B_A𝑓_𝑥𝑑𝑥𝑑𝑦, (3) where 𝐶 is a constant normalization factor. Considering the sparse number of coincidences contributing to each ⟨𝑃_ab(𝜃)⟩, we derived the maximum-likelihood expectation maximization (ML-EM) [37, 38] in a list mode in which each coincidence will be stored and calculated successively. Thus, 𝑓^→_𝑖^𝑛+1=𝑓^𝑛_𝑖1∑^𝑁_𝑗′=1𝑎_𝑗′^𝑁_,𝑖𝑗=1𝑎_𝑗,𝑖1∑^𝑝_𝑖′=1𝑎_𝑗,𝑖′𝑓^𝑛_𝑖′,𝑖=1,…,𝑝, (4) where 𝑝 is the total number of pixels in the image, 𝑁 is the total number of detected coincidences, and 𝑎_𝑗,𝑖 is the element of system matrix representing the probability, that the two annihilation photons emitted from pixel 𝑖 will be detected as the 𝑗th coincidence (whether scattered or not) in a list mode entry weighted by the Compton differential cross section. In practice, 𝑎(𝑖,𝑗) is proportional to the Compton cross-section if the 𝑖th pixel falls within the 𝑗th TCA and is zero if 𝑖 is outside the TCA. The main difference between this GS-method-based ML-EM algorithm (GS-MLEM) and the conventional ML-EM algorithm (LOR-MLEM) is that the summation for each coincidence is over the area confined by the TCA, instead of being along the LOR. This algorithm also can account for random coincidences 𝑅 by adding them to the projector as follows: ∑^𝑝_𝑖′=1𝑎_𝑗,𝑖′𝑓^𝑛_𝑖′+𝑅_𝑗. The GS-MLEM was implemented in MATLAB as shown in Figure 3. The ratio of the measured number of coincident events (which is 1 with a list-mode algorithm) to the sum of the pixel values within each TCA for the current estimate is backprojected into a “ratio” matrix. This backprojection is repeated for all coincident events. The “ratio” matrix is used to modify the current estimate. This process is repeated in an iterative fashion, updating the current image until a good convergence is obtained.

Figure 3 

The GS-MLEM algorithm flow chart. The additions and changes to the conventional MLEM (LOR-MLEM) algorithm are highlighted by the boxes with a black thick outline.

2.3. The GATE Platform and Phantom Measurements

The proposed algorithm was evaluated by simulating a 2D small animal PET system using GATE [39]. GATE is an object-oriented Monte Carlo simulation platform based on Geant4 libraries (a generic Monte Carlo code) [40].  The simulated PET scanner had a 24 cm diameter detector ring made up of 42 detector blocks. Each detector block is arranged in a 9×5 lutetium oxyorthosilicate (LSO) crystal array. Each crystal element had a surface area of 2×2 mm² and was 18 mm thick. The energy resolution of the simulated scanner in this preliminary investigation was 0.1% full width at half maximum (FWHM) at 511 keV, and a 170–511 keV energy window, in which all single scattered coincidences can be selected, was used.  Noise and dead time were not included in the simulation process in order to focus on the role of that scattered coincidences played in the image quality.

A simplified Deluxe Jaszczak phantom [41, 42]  having 3 hot disks and 1 cold disk with radii of 1.5 mm, 3 mm, 4.5 mm, and 6 mm within a cylindrical water phantom with a 40 mm radius and a hot-to-background ratio of 4 as shown in Figure 4 was used.

Figure 4 

A simplified Deluxe Jaszczak phantom. The three yellow circles represent sources with radii of 1.5 mm, 3 mm, and 4.5 mm, respectively. The black circle represents a cold area with radius 6 mm. The red area represents the 40 mm water phantom. The hot-to-background ratio is 4.

Images were reconstructed within a 99×99 pixel matrix, with a pixel size of 1×1 mm². The contrast recovery ratio of the reconstructed images was analyzed using the mean value of the disks relative to the local background while the noise was calculated using the relative standard deviation (RSD) of a 10 mm diameter circular region of interest (ROI) located in the center of the phantom image. Two concentric circular ROI’s with diameters of 6 mm and 8 mm were defined within the #3 and #4 disks, respectively. Annular ROI’s with a 20 mm outer diameter and an inner diameter of 9 mm and 12 mm respectively, were drawn in the background around disks #3 and #4. The local contrast recovery coefficient (CRC_local) was defined by CRC_local=(𝐻/𝐵−1),𝑅−1 (5) where 𝐻 is the average values within the ROI in the hot disk, 𝐵 is the average value in the annulus around the disk, and 𝑅 is the experimentally simulated hot-to-background ratio. For the cold source, CRC_local is defined by CRC_local𝐶=1−𝐵, (6) where 𝐶 was the average value of the ROI in the cold disk.

For the purpose of evaluation, a point on the CRC curve was determined which was the shortest distance from the CRC curve to the point defined by an ideal CRC of 1 and a RSD of 0.

2.3.1. The Feasibility Test of Images Reconstructed Using Only Scattered Coincidences

To illustrate that scattered coincidences can indeed be used to reconstruct images, 10⁶ scattered coincidences were used to reconstruct images using the GS-MLEM method. By comparison the conventional ML-EM (LOR-MLEM) method which takes scattered coincidences as true was used to reconstruct 10⁶ scattered coincidences as well as 10⁶ true coincidences generated using the GATE simulation.

2.3.2. Comparison of Images Reconstructed by GS-MLEM and LOR-MLEM with Different Scattering Fraction Data

The scatter fraction is a function of the geometry of scanner, the density, and the size of the patient as well as the energy window employed. To compare images reconstructed with different scatter fractions, not affected by the above factors, we simulated 3×10⁵ true coincidences and added scattered coincidences to obtain the required data with scatter fractions ranging from 0% to 60% using the same phantom.  The simulation data we used is more artificial than real; however, it serves the goal of this research to investigate the relative value of the proposed algorithm for different fractions of scattered coincidences to the image quality while at the same time removing unwanted complexities resulting from more realistic simulations.

2.3.3. Evaluating the Different Energy Thresholds of Scattering Coincidences on Image Quality

The scatter fraction affects contrast and noise of images. As more scattered coincidences are included in the image reconstruction, the noise will decrease but so will the contrast.  However, the strength of the proposed method lies in its ability to recover contrast while keeping noise low.  The tradeoff between contrast and noise (for a particular scatter fraction) can be adjusted by varying the scattered photon energy threshold. Photons scattered through large angles may still be of value close to the periphery of the phantom. In this paper, those coincident photons which were scattered in an area defined by the intersection of the area encompassed by the TCA and the image matrix which was less than some predetermined threshold were used.  This is similar to using a scattered energy threshold, but it enables both energy and detector position to be used to full advantage.  A total of 6×10⁵ coincidences with a 50% scattering fraction were used to test how the varying thresholds enabled the GS method to trade off between the contrast and noise. By setting different thresholds and calculating the intersection areas between TCA and image matrix, the data was generated and used for reconstruction.

3. Results

3.1. Feasibility Test of Reconstructing Images Only by Scattered Coincidences

Figure 5(a) displays a typical image produced from true (unscattered) coincidences using a conventional ML-EM (LOR-MLEM) algorithm. Figure 5(b) shows that the same algorithm was incapable of producing recognizable images from scattered photons only since the LOR assigned for each scattered coincidence no longer passes through the annihilation point. This illustrates why scattered photons are traditionally considered contributing only noise to the resultant image and are, where possible, typically eliminated. Our initial experiments demonstrated that the proposed algorithm was capable of producing an image entirely from scattered photons and Figure 5(c), while not as good as the image in Figure 5(a), still adequately represents the activity distribution. There is still blurring evident in the scattered photon reconstruction and image converges at least 2-3 times slower. From the profiles under the corresponding images, we can see that the image produced using only scattered events has a decreased contrast and noisier background compared with that using the same number of true coincidences. Figure 5(c) was used to illustrate the feasibility of reconstructing an image from only scattered coincidences. However, in practice, an image would be reconstructed by incorporating both true and scattered events instead of only using scattered coincidences. We hypothesized that including scattered coincidences in the reconstruction provides more advantages than simply rejecting them and the results will be shown in Section 3.2.

(a)

(b)

(c)

Figure 5 

(a) The image created by reconstructing 106 true coincidences using an LOR-MLEM with 36 iterations. (b) The image reconstructed using 106 scattered coincidences and an LOR-MLEM approach which assumes that these photons are true coincidences.  (c) It uses the same scattered coincidences dataset as in (b) but reconstructs them using a GS-MLEM approach with 62 iterations. The second row plots the profiles through the above images in the horizontal and vertical direction passing through the center of the images, respectively.

3.2. Comparison of Images Reconstructed by GS-MLEM and LOR-MLEM with Different Scattering Fraction Data

Figure 6 shows the CRC curves for the #3 (largest hot disk) versus the relative standard deviation of the background obtained by varying the number of iterations. Figure 7 is similar except that it shows the results for #4 (the cold) source. The CRC curves for images reconstructed using the GS-MLEM algorithm were always above those created by the LOR-MLEM algorithm.  Unlike the CRC curves for the LOR-MLEM algorithm which decrease with increasing scatter fraction, the CRC curve for the GS-MLEM algorithm generally increased with increasing scatter fraction. For the hot disk, this trend reversed beyond the point where the CRC curve for the GS-MLEM intersected with that of the CRC for reconstructions using only true coincidences. For the cold disk, this reversal was not observed. This reduction in CRC was small and occurred only for iterations beyond the evaluation point.  For scatter fractions from 10% to 60%, the evaluation points on the CRC curve for the hot disk were 0.5–3.0% greater than those obtained using only true coincidence data only and were 3.0–24.5% greater than the corresponding curves reconstructed using the LOR-MLEM algorithm.  The noise was reduced by 0–1.7 % in comparison with the true coincidences data and was 2.0–12.0% less than that produced by LOR-MLEM method with the same scatter fraction. For scatter fractions between 10% and 60%, the evaluation points for the cold disk had a CRC 0–3.5% greater than the curves calculated using only true coincidences and were 4.0–24.0% greater than the LOR-MLEM method. The noise at the evaluation point was 0–1.6% less than that obtained using only true coincidences and was 0.5–8.3% less than that calculated using the LOR-MLEM method with the same scatter fraction. This trend is consistent with the position that the scattered coincidences contribute to noise in the conventional PET image reconstruction and that as the scatter fraction increases, the image quality deteriorates.

Figure 6 

The CRC curves for #3 disk versus the relative standard deviation of the background. The CRC curves shown with filled symbols represent the images reconstructed by GS-MLEM with different scattering fractions. The lines in the same style with empty symbols correspond to images with the same scatter fraction, reconstructed by LOR-MLEM. The CRC for true coincidences is shown as a solid line for comparison. The diamonds reflect the optimal points on the CRC-RSD curves.

Figure 7 

The CRC curves for the #4 (cold) disk versus the relative standard deviation of the background. The CRC curves shown with filled symbols represent the images reconstructed using the GS-MLEM with different scattering fractions. The lines in the same style with empty symbols correspond to images with the same scatter fraction, reconstructed by LOR-MLEM. The CRC for the true coincidences is shown as a solid line for comparison. The diamonds reflect the optimal points on the CRC-RSD curves.

3.3. Evaluating the Different Energy Thresholds of Scattering Coincidences on Image Quality

Figure 8 plots the CRC of #3 (the largest source) for different TCA area thresholds as a function of the relative background standard deviation obtained by varying the number of iterations. Figure 9 is similar, except that it shows the results for the #4 (cold) source. The results show that the value of the CRC for images reconstructed using the GS-MLEM method were always greater than those obtained with a conventional LOR-MLEM approach which assumed that all coincidences that fell into the typical 350–511 keV energy window were true. For the GS-MLEM method, the evaluation point (point #1 in Figure 8) for the source had a CRC that was 2.5% better than images reconstructed with true coincidence data only (which is an idealistic upper limit on existing algorithms) (point #2 in Figure 8) and was 13.0% greater than images reconstructed using the conventional LOR-MLEM algorithm (point #3 in Figure 8). The image noise was 2.0% less (point #1 in Figure 8) than the ideal reconstruction using true coincidences (point #2 in Figure 8) and was 7.0% less than the LOR-MLEM approach (point #3 in Figure 8). For the cold source, the evaluation point for the GS-MLEM approach (point #1 in Figure 9) had a CRC that was 5.0% greater than the CRC obtained with the ideal true coincidence calculation (point #2 in Figure 9) and 18.0% greater than the LOR-MLEM method (point #3 in Figure 9). For the cold source, the noise at the evaluation point using the GS-MLEM method (point #1 in Figure 9) was 2.0% less than the ideal reconstruction using trues (point #2 in Figure 9) and was 8.0% less than the LOR-MLEM approach (point #3 in Figure 9). Virtually no difference was found between the CRC curves for thresholds of 91.8% and 100%. The evaluation CRC/noise point for the source occurred for thresholds close to 100% of the matrix size while the evaluation point for the cold source was optimal when the threshold was set to 10.2% of the matrix size. The evaluation contrast recovery point for the source approached the true value as the threshold approaches zero, but the CRC for the cold source still showed an improvement for thresholds approaching 100%.  

Figure 8 

The CRC curves for the #3 (largest hot) disk calculated when only trues are present and with 50% scatter fraction calculated using the conventional LOR-MLEM approach and with the GS-MLEM approach using different thresholds defined by the ratio of the area of intersection between the TCA and image matrix to the total image matrix area. Points #1, #2, and #3 identify the evaluation points for reconstructions using GS-MLEM method, only trues, and the LOR-MLEM, respectively.

Figure 9 

The CRC curves for the #4 (cold) disk calculated when only trues are present and with 50% scatter fraction calculated using the conventional LOR-MLEM approach and with the GS-MLEM approach using different thresholds defined by the ratio of the intersection area between the TCA and image matrix to the total image matrix area. Points #1, #2, and #3 identify the evaluation points for reconstruction using GS-MLEM method, only trues, and the LOR-MLEM, respectively.

Figure 10(a) displays the image reconstructed from 3×10⁵ true coincidences using a conventional LOR-MLEM method. Figure 10(b) shows the image reconstructed using the same true 3×10⁵ coincidences plus 3×10⁵ scattered coincidences that fall into the 350–511 keV energy window and reconstructed using a conventional LOR-MLEM algorithm as a comparison. Figure 10(c) shows the image produced by 6×10⁵ coincidences with 50% scatter fraction and a threshold of 100% reconstructed using the GS-MLEM method.

(a)

(b)

(c)

Figure 10 

(a) The image using 3×105 true coincidences reconstructed using an LOR-MLEM reconstruction algorithm. (b) The image using the same true 3×105 coincidences plus 3×105 scattered coincidences that fall into the 350–511 keV energy window reconstructed using the conventional LOR-MLEM algorithm as a comparison. (c) The image using 6×105 coincidences with a 50% scattering fraction reconstructed using the GS-MLEM algorithm. The second row shows the profiles of the above images along horizontal and vertical direction passing through the center of the images, respectively.

4. Discussion

We have described a novel PET reconstruction algorithm which can incorporate both true and scattered coincidences into the reconstruction process.  Our initial results have shown that under the ideal conditions used in this study, including scattered coincidences in the reconstruction is more advantageous than simply rejecting them.

It has been shown that an image of the activity distribution can be reconstructed from scatter data only as shown in Figure 5(c). The image in Figure 10(c) is still idealistic and far from optimized, but these results show that the proposed method produces an image with a more homogeneous background, sharper edges, reduced noise, and improved contrast (for both hot and cold disks), relative to both conventional algorithm (which degrade in the presence of scatter) and the idealistic situation (only true events used). Algorithms which incorporate scatter correction will move closer to the idealistic response curve but will never exceed it. Thus, including scattered photons directly into the reconstruction could eliminate the need for (often empirical) scatter corrections required by conventional algorithms and increase image contrast and signal to noise ratio (SNR). This could be used to either improve the diagnostic quality of the images and/or to reduce patient dose and radiopharmaceutical cost.

The results in Section 3.3 show that an additional advantage of the GS-MLEM algorithm is that a threshold can be used to adjust contrast and noise in order to optimize the role that the scattered coincidences play in the image reconstruction process. At the optimal threshold, adding scattered coincidences into imaging reconstruction appears to result in minimal (if any) loss in contrast while significantly improving the SNR and contrast recovery. The results show that by appropriately including a greater number of scattered photons in the reconstruction, improvements in both contrast and noise can be achieved. In addition, this threshold also can remove partial coincidences scattered within detectors or in the gantry, which usually have relatively larger scattering angles and unusually short distances between scattering positions and detector positions.  The TCA for these situations usually covers the whole image space and can be removed from the reconstruction dataset by employing the proposed threshold. Our current work has assumed that intercrystal scattering is small when compared to scattering within the patient. While LSO crystals have been used in this simulation, in future we will use a dual layer detector to minimize inter-crystal scattering.

Similar to conventional energy-based scatter correction methods, the accuracy of the proposed method is limited by the energy resolution of the detector. In conventional scatter-correction-based methods, the ability to distinguish and reject scattered coincidences is limited by energy resolution.  In the proposed method, the energy resolution of detectors determines the confidence with which the TCA can be defined for each scattered coincidence, which in turn affects the locus of the possible scattering positions and hence the annihilation position. The locus is sensitive to the energy of the detected photons, and a small uncertainty in the energy of the scattered photon can result in a significant difference in the shape of the TCA and the position of the source.  For example, if the detector’s energy resolution is 4% at 511 keV, the maximum uncertainty in the calculated area between the TCA is approximately 21% for a scattered photon energy of 450 keV. Thus, developing a detector with high energy resolution and weighting the spatial probability distribution of the annihilation position by using both the patient outline [43] and the energy resolution of the detector will be important future work.

Guerin et al. indicated that neglecting the effect of multiple scattered coincidences was not a serious source of error. In practice, multiple scattered events cannot be distinguished from single scattered events on the basis of the energy of the scattered photons. Even though the relation between scattered energy and scattering angle cannot be connected by Compton scattering equation for multiple scattered coincidences, we propose to use a synthetic scattering angle for multiple scattered events based on the scattered photon energy. In addition, the uncertainty in the energy of the detected photon as a result of the limited energy resolution of the detector is more challenging for multiple scattered events when trying to determine the locus of possible scattering positions.  Corrections for randoms can be addressed by adding a factor to the forward projection in the algorithm as do most existing methods. The data model used to derive the reconstruction algorithm can be improved by including the electron density and tissue attenuation. The impact of these on the quality of the reconstructed images is beyond the scope of this paper and will be the focus of a separate study. Since the TCA is a function of scattering angle and detector positions, which will be calculated for each detected coincidences in the reconstruction algorithm, the computational time is currently 3-4 times longer when compared with a conventional MLEM algorithm. However, this algorithm still needs to be optimized for calculation efficiency.

The uniform attenuation phantom used in this work is simple in comparison to anthropomorphic phantoms and may not reflect the complexity of scatter in humans. However, the purpose of this work was to compare the results of the proposed algorithm as a function of the scatter fraction which would be difficult to analyze with an anthropomorphic phantom.

5. Conclusion

A new method which includes scattered photons directly into the reconstruction has been presented and evaluated in this paper. The results of the phantom study in this paper have shown that PET images can be reconstructed from scattered coincidences. Including scattered coincidences into the reconstruction eliminates the need for scatter corrections while increasing image contrast and reducing noise. A threshold, which depends on both energy and detector position, is used, adjust noise and contrast was also evaluated in this work. The optimal threshold was different for hot and cold sources, but the variation in CRC for the cold source was only weakly dependent on the threshold. Improvements in the CRC and noise for both hot and cold sources could be obtained by maximizing the use of all scattered photons.

Acknowledgments

This work was supported in part by CancerCare Manitoba Foundation, University of Manitoba, and the Natural Sciences and Engineering Research Council of Canada (NSERC).