Evaluation of Image Quality Improvements When Adding Patient Outline Constraints into a Generalized Scatter PET Reconstruction Algorithm

Scattered coincidences degrade image contrast and compromise quantitative accuracy in positron emission tomography (PET). A number of approaches to estimating and correcting scattered coincidences have been proposed, but most of them are based on estimating and subtracting a scatter sinogram from the measured data. We have previously shown that both true and scattered coincidences can be treated similarly by using Compton scattering kinematics to define a locus of scattering which may in turn be used to reconstruct the activity distribution using a generalized scatter maximum-likelihood expectation maximization (GSMLEM) algorithm. The annihilation position can be further confined by taking advantage of the patient outline (or a geometrical shape that encompasses the patient outline). The proposed method was tested on a phantom generated using GATE. The results have shown that for scatter fractions of 10–60% this algorithm improves the contrast recovery coefficients (CRC) by 4 to 28.6% for a source and 5.1 to 40% for a cold sourcewhile the relative standard deviation (RSD)was reduced. Including scattered photons directly into the reconstruction eliminates the need for (often empirical) scatter corrections, and further improvements in the contrast and noise properties of the reconstructed images can be made by including the patient outline in the reconstruction algorithm as a constraint.


Introduction
Scattered photons are a significant source of image quality degradation and lead to quantitative errors in positron emission tomography (PET) [1].The scatter fraction can be as high as 40-60% when a tomograph operates in 3D mode without slice-defining septa and in large patients, making this of greater consequence for cardiac imaging [2][3][4].In conventional PET reconstruction methods, scattered coincidences are assumed to be noise and consequently a number of ways for estimating and correcting scattered coincidences in measured data have been proposed.However, most of these techniques are based on the estimation and subtraction of scatter from the projection data instead of exploring the possibilities of using scattered coincidences in the reconstruction [2,5].During the scatter correction process, artifacts may be introduced in the source distribution due to inaccurate estimation of the scatter sonogram [1] while at the same time, the system's sensitivity will be reduced and image noise will be amplified [2,6].
The energy resolution of PET detectors has improved in [7,8] recently, making it conceivable to use the energy and detected location of the coincident photons to develop a spatial distribution of the annihilation positions.Scattered coincidences are thus a latent source of concealed information which can be used to improve PET image quality.
In our previous study, we have shown that true coincidences can be considered to be a subset of scattered coincidences and that a GS-MLEM algorithm can use both true and scattered events to reconstruct the source distribution [9].This method takes advantage of the kinematics of Compton scattering by connecting the coincidence detectors with two arcs (2D) (rather than the conventional straight line assumption which is only true for unscattered photons) which describe the locus of scattering.It can be shown that the point of annihilation is encompassed in 2D by the two circular arcs (TCA) or in 3D by a surface, described by the rotation of these arcs around an axis connecting the detectors; see Figure 1.Including scattered coincidences directly in the image reconstruction eliminates the need for scatter correction.It also improves the image quality, particularly when the data is sparse and may increases the system's sensitivity since a lower energy threshold can be used and more data will be measured [10].This work has been presented, in part, at the World Congress on Medical Physics and Biomedical Engineering in 2012 [11] and the 2012 IEEE Nuclear Science Symposium and Medical Imaging Conference [12].
We hypothesize that the annihilation position can be further confined by making use of the patient outline (or a geometrical shape that encompasses the patient outline) as a further spatial constraint.In this work, we evaluated the contrast and noise properties of the constrained GS-MLEM algorithm with a patient/phantom outline constraint as well as the dependency of the proposed method on the accuracy of the patient/phantom outline constraints employed.

Outline Constraint Reconstruction Theory.
In PET, the three main sources of photon scatter are (1) the patient, (2) detectors, and (3) the gantry and surrounding environment.Patient scatter generally dominates in human imaging, and in this work we assume that the patient is the only scattering source and ignore the relatively small contribution due to scattering in the gantry and detectors [13,14].
In such a scenario, the possible annihilation positions can be further constrained by connecting the intersection points between the TCA and the patient outline, C (the furthest from the unscattered photon detector A) and D (closer to the unscattered photon detector A), with the unscattered photon detector A. The position of annihilation is confined to the area encompassed by the TCA, the patient/phantom outline, and the line AC (area CDE as shown in Figure 2).If both circular arcs of a TCA for a coincidence intersect with the patient/phantom outlines, the areas used to confine the annihilation positions can be calculated for each circular arc separately in a similar way.
The patient outline can be estimated by a variety of means.One way is to use an external X-ray source or optical system (either laser based or photogrammetric).In PET/CT or PET/MRI systems the anatomical image provided by the CT or MRI could be used.Alternatively it may be possible to estimate the constraints using an approximation of the patient outline based on initial iterations or by using a basic geometric shape (say a circle or ellipse) as shown in Figure 3.The size of the circle could be chosen based on a simple measurement (or estimate) of the patient size or based on the early iteration images.

Constrained GS-MLEM Algorithm.
Before building the patient outline constraint into the GS-MLEM algorithm, the expected number of coincidences in which an unscattered photon is observed at A while the other photon is observed at B following a Compton scattering through an angle  is determined from In the above expression the total number of photons which could reach point  is obtained by integrating the source   from line segment  to .The possibility that the photons at  undergo a Compton scattering event is proportional to the electron density at this point weighted by the differential Klein-Nishina electronic cross section.The outer integral is calculated over all possible scattered points which have been constrained by the patient outline, where  is the acquisition time and  is the linear attenuation coefficient.
To reduce the large computational workload, we ignore the attenuation and the electron density and assume that the expected number of detected coincidences is linearly proportional to the product of the differential Klein-Nishina electronic cross section and the total activity within area CDE as given by where  is a constant.Putting this data model into the system model and following the same maximization process as in [4], we can derive the constrained GS-MLEM algorithm in list mode [15,16] as follows: where  is the total number of pixels in the image,  is the total number of coincidences,  is the iteration number, and  , is the system matrix accounting for the probability that the th coincidence detected in the list mode entry comes from pixel .To incorporate the patient outline as a constraint in the algorithm, the sum over   in the denominator of the second factor in the right hand side is only calculated within the confined areas.

Phantom Evaluation.
To make the comparison consistent with our previous study [9], the same PET scanner and phantom (see Figure 4) are simulated by GATE in this work [17,18].The simulated system had a perfect energy resolution with the energy window of the system set from 170 to 511 keV.Noise and dead time were not modeled in the simulation process.In this initial evaluation of the proposed method, only the coincidences scattered within the phantom were used in the reconstruction.The image quality was evaluated using the contrast recovery coefficient (CRC) and relative standard deviation (RSD) which reflects the contrast and noise properties of the reconstructed images.The local contrast recovery coefficient (CRC) for hot disk was defined by where  is the average value within the hot disk,  is the mean value of the background, and  is the experimentally set hotto-background ratio.Similarly, CRC local is defined for a cold area by where  represents the mean value in the cold region.
An evaluation point on the CRC curves was determined by the shortest distance from the CRC curve to the point defined by the ideal CRC = 1 and RSD = 0.

Evaluation of Images Reconstructed Using the Proposed
Method with Different Scattering Fraction Data.The image quality of the image reconstructed with the phantom outline constraint was compared to the results using GS-MLEM without the phantom outline constraint and to the conventional LOR-MLEM method given in [9].Coincidences with scattering fractions ranging from 0 to 60% were randomly simulated by GATE.To be consistent with our previous results, 3 × 10 5 true coincidences and scattered coincidences with the required scatter fraction were generated.Here, the actual phantom outline (a circle with a radius 40 mm) was employed in the reconstruction process as the outline constraint.

Evaluation on the Different Outline Constraint Approximations on Image
Quality.The improvement in the contrast and noise properties of the reconstructed images was related to the accuracy of the patient/phantom outline used.When the phantom outline constraints chosen are smaller than the actual phantom outline, the confined area may not encompass the annihilation position, which will reduce the image contrast and introduce artifacts.Thus, the minimum but still most accurate outline constraint would be the actual patient/phantom outline.A larger area can be used to constrain the annihilation position but will not be optimal.When the constraint outline approaches that of the detector positions, the outline constraint based GS-MLEM algorithm will approach the nonoutline constraint based algorithms.To evaluate the effect of different constraint sizes on the reconstructed image quality, we characterize the patient/phantom outline sizes as a function of the ratio of difference between the tested outline constraint and the actual phantom outline size.We therefore tested the effect of various phantom outline constraints using a circle with radii of 42 mm, 45 mm, 50 mm, and 60 mm, being 5%, 12.5%, 25%, and 50% larger than the actual phantom outline, respectively.A total of 6 × 10 5 coincidences with 50% scatter fraction were generated by GATE and reconstructed using the proposed method with the different outline constraints.The images were also reconstructed using the same dataset but with GS-MLEM without the phantom outline constraint and with the conventional LOR-MLEM algorithm, as well as with 3 × 10 5 true coincidences using conventional LOR-MLEM algorithm as a comparison.versus RSD for different methods/or data for the no. 3 (the largest hot) disk in Figure 5 and the no. 4 (cold) disk in Figure 6, respectively.The results show that the CRC curves for images reconstructed using GS-MLEM with phantom outline constraint method are always above those of the unconstrained GS-MLEM method and the conventional LOR-MLEM approach.The optimal CRC curve for the conventional LOR-MLEM method is at zero scatter fraction, and the curves decrease with increasing scattering fraction.In contrast to that, the CRC curve for the GS-MLEM, both with and without phantom outline constraint, generally increases with increasing scatter fraction.This trend, for the hot disk, changes beyond the point where the CRC curves for the constrained and unconstrained GS-MLEM intersect.However, this contrast reduction is not obvious and only occurs beyond the evaluation point.The same trend for the cold disk is not observed.The figures show that the GS-MLEM algorithm with phantom outline constraints improved the CRC properties of the reconstructed images by 0.6 to 3.8% compared with GS-MLEM without the phantom outline constraint for scatter fraction ranging from 10% to 60%.For the same scatter fractions, the results were 4% to 28.6% greater than the corresponding curves reconstructed using the LOR-MLEM algorithm.The noise of images reconstructed using the GS-MLEM algorithm with phantom outline constraint was 0.7% to 2.6% lower than the corresponding curves using an unconstrained GS-MLEM algorithm and was 2.4 to 14.7% less than that produced by the LOR-MLEM method.For the cold disk, with a scatter fraction of 10% to 60%, the evaluation point for the cold disk using GS-MLEM with phantom outline constraint had   (c) Figure 7: (a) shows the image using 6 × 10 5 coincidences with a 50% scattering fraction by the proposed method (GS-MLEM plus patient/phantom constraint method).(b) shows the image using the same data as (a) and by GS-MLEM without phantom outline constraint; (c) shows the image using the same true 3 × 10 5 coincidences plus 3 × 10 5 scattered coincidences that fall into the 350 to 511 keV energy window and was reconstructed using the conventional LOR-MLEM algorithm as a comparison.The second row shows profiles of the above images passing through the center of the images in the horizontal and vertical directions, respectively.a CRC of 1.1% to 11.6% greater than the curves calculated only using GS-MLEM and was 5.1% to 40% greater than the LOR-MLEM method for the corresponding scatter fraction.The noise at the evaluation point for images reconstructed using GS-MLEM with phantom outline constraint was 0.5% to 3.6% less than that obtained using an unconstrained GS-MLEM and was 1.5% to 11.8% less than that calculated using the LOR-MLEM method with the same scatter fraction.

Evaluation of Images
Figure 7(a) displays the image reconstructed from 6 × 10 5 coincidences with a 50% scatter fraction using the GS-MLEM method, constrained by the phantom outline.Figure 7(b) shows the image produced using the same number of coincidences and algorithm as in Figure 7(a) but without the outline constraint.As a comparison, the same 3 × 10 5 true coincidences were combined with scattered coincidences from a 3 × 10 5 dataset and reconstructed using a conventional LOR-MLEM algorithm with a 350 to 511 keV energy window.The results shown in Figure 7 demonstrate that constrained GS-MLEM (Figure 7(a)) has the best image quality with a more uniform background and sharper edges than those of the unconstrained GS-MLEM (Figure 7(b)) which in turn has a better image quality than that of the conventional LOR-MLEM method (Figure 7(c)).

Evaluation of the Dependency of the Proposed Algorithm on the Accuracy of the Outline Constraints.
Images with different phantom outline constraints were also reconstructed to evaluate the dependency of the proposed algorithm on the accuracy of the phantom outline constraints.The CRC curves of the no. 3 (the largest hot) disk for different phantom outline constraints as a function of the relative background standard deviation were obtained by varying the number of iterations as shown in Figure 8.The same result for the no. 4 (cold) disk can be seen in Figure 9.For the hot disk, the CRC curves for the GS-MLEM method with phantom outline constraints decrease as the phantom outline constraints increase up to the point where the CRC curves for the GS-MLEM method intersect with those of true coincidences.This intersection was not observed for the cold disk, but the same trend was observed.All the CRC curves with phantom outline constraints were above those that did not use phantom outline constraints.The CRC and RSD for the evaluation points on each CRC curve are plotted as a function of the relative increase of radii of the phantom outline constraints to the phantom's actual radius, for the hot and cold disks in Figures 10 and 11, respectively.As the relative radius increased from 5% to 50%, the CRC for the hot disk reduced by 2% while noise increased by 0.5%.For the cold disk, the CRC reduced by 4.5% and noise increased by 1.3%.For both the hot and cold disks, the CRC varied slowly when the relative radius was increased by less than 12.5% (corresponding to an outline constraint with radius 45 mm), when compared to the change for the relative radius increase from 12.5 to 50%.

Discussion
In our previous study [9][10][11] and this work, we have shown that scattered coincidences are not just noise and can be used in the imaging reconstruction by taking advantage of the individual photon energy.Conti et al. have taken advantage of the time difference of scattered photons in a time of flight (TOF) PET and applied this knowledge to refine the physics model to describe the annihilation position [10].In this work, we included the patient/phantom outline as a constraint to the physics model to further confine the possible annihilation position.The results have shown that the contrast of the reconstructed image was improved while noise was reduced.Since this method no longer considers scatter coincidences as noise, more data is available, and a lower energy window can be applied, thus improving both the system sensitivity and image quality.
The results for the evaluation of different outline constraints on image quality have shown that uncertainties in defining the patient/phantom outline are not a significant obstacle when implementing the proposed method.This method is based on the assumption that the patient scattering dominates over detector, gantry, and surrounding environment scattering, which is true in human imaging.
In this initial test of the proposed method, we carried out this work in 2D, and only the coincidences scattered within phantom were selected to show the principle of the proposed method and at the same time to reduce the complexity of the mathematics.This approach could be implemented in 3D where it could be of even greater value.The TCA are calculated using the detector positions and the scattering angle which is closely related to the scattered photon energy using Compton equation.The accuracy with which the locus of the TCA assigns the possible scattering positions and encompasses the annihilation position will therefore depend on the energy resolution of the detectors.The area defined using a scattered photon energy with a large uncertainty may not encompass its annihilation position or may overestimate the confined area, resulting in artifacts or blurring of the reconstructed images.The effects of energy resolution on the proposed method will be investigated in the future work.

Conclusion
Previous work demonstrated that scattered coincidences can be directly included in the reconstruction process.The results of this study show that further improvements in the contrast and noise properties of the reconstructed images can be made by including the patient outline in the reconstruction algorithm as a constraint.While these results are promising, we are currently investigating ways of correcting for attenuation and improving the limited energy resolution of current detectors, both of which will need to be resolved before scatter reconstruction can be used for clinical applications.

Figure 1 :
Figure 1: A diagram illustrating a Compton scattering event.A positron annihilates at the green dot and generates two 511 keV photons.One is observed by detector A and the other one undergoes a Compton interaction and is finally detected by detector B. The possible scattering positions can be described by the two circular arcs (TCA) defined by the kinematics of Compton scattering.

Figure 2 :
Figure 2: In this case, one of the TCA interacts with the patient at point D (closer to A) and C (further from A).If the extent of the scatter volume is known (the patient outline), the possible annihilation area may be further confined to the area C-D-E encompassed by line CA, TCA, and the outline of the patient.

Figure 3 :
Figure 3: The patient outline is replaced by an ellipse which is slightly larger than the patient outline and intersects with TCA at point C (the furthest from the unscattered photon detector) and point D (closer to the unscattered photon detector).The possible annihilation positions are confined to the area C-D-E encompassed by line CA, TCA, and the ellipse calculated in the same way as the case of patient/phantom outline.

Figure 4 :
Figure 4: A cylindrical water phantom with three hot areas (yellow color) and one cold area (black color) was simulated.The radii for the 1 to 4 circles are 1.5 mm, 3 mm, 4.5 mm, and 6 mm, respectively.We set the activity ratio R between the hot and background to 4.

Figure 5 :
Figure5: The CRC curves of no. 3 (the largest hot) source were calculated with GS-MLEM with phantom outline constraints, GS-MLEM method and the conventional LOR-MLEM approach.

Figure 6 :
Figure6: The CRC curves of no. 4 (cold) source were calculated with GS-MLEM with phantom outline constraints, GS-MLEM method and the conventional LOR-MLEM approach.

Figure 8 :
Figure8: The CRC curves of no. 3 (largest hot) disk calculated using GS-MLEM with different phantom outline constraints.The CRC for 3 × 10 5 true coincidences and 6 × 10 5 coincidences with 50% scattered fraction using conventional LOR-MLEM were also reconstructed as a comparison.

Figure 9 :Figure 10 :
Figure9: The CRC curves of no. 4 (cold) disk calculated using GS-MLEM with different phantom outline constraints.The CRC for 3 × 10 5 true coincidences and 6 × 10 5 coincidences with 50% scattered fraction using conventional LOR-MLEM were also reconstructed as a comparison.

Figure 11 :
Figure 11: The CRC and noise properties of no. 4 (cold) disk as a function of relative increase of radii of the phantom outline constraints for the evaluation points of Figure 9.
Reconstructed Using the Proposed Method with Different Scattering Fraction Data.To illustrate the performance of the proposed method, we plot the CRC