Double scatter simulation for improved reconstruction of positron emission tomography data

Abstract

Methods for simulating, and correcting for, doubly scattered annihilation gamma-ray photons in both time-of-flight (TOF) and non-TOF positron emission tomography scan data are disclosed.

Claims

1. A method for applying scatter correction to a scan data acquired in a time-of-flight positron emission tomography (TOF PET) scanner, the method comprising: (A) selecting a pair of detector positions in the TOF PET scanner's detector ring; (B) numerically estimating the double scatter coincidence event rate of annihilation photons in the selected pair of detector positions, wherein the annihilation photons in a double scatter coincidence event are assumed to be scattered once, each photon by a different scatter point in the TOF PET scanner's image volume, during a TOF PET acquisition scan period; (C) selecting another pair of detector positions; (D) repeating step (B) wherein the selected pair of detector positions is the another pair of detector positions; (E) repeating steps (C) and (D) until all detector pairs of interest have been selected; and (F) scatter correcting the acquired scan data obtained from the pair of detector positions during the TOF PET acquisition scan period to reconstruct a TOF PET image based on the scan data.

2. The method of claim 1, wherein the selected pair of detector positions are D.sub.1 and D.sub.2, wherein the step (B) comprises: a) selecting the two scatter points, S.sub.1 and S.sub.2, in the TOF PET scanner's image volume; b) sampling values of the emission rate density between the two scatter points S.sub.1 and S.sub.2, determining the time of flight offset bin, n, for each contribution, and storing these values in a vector R.sub.n.sup.12; c) calculating an exponential function of the negative of the line integral of the linear attenuation coefficients at energy E.sub.0 between the two scatter points S.sub.1 and S.sub.2; d) calculating the exponential function of the negative of the line integral of the linear attenuation coefficients at energy E.sub.1 between S.sub.1 and D.sub.1; e) calculating the exponential function of the negative of the line integral of the linear attenuation coefficients at energy E.sub.2 between S.sub.2 and D.sub.2; f) calculating the reciprocals of the squares of the distances from S.sub.1 to S.sub.2, from S.sub.1 to D.sub.1, and from S.sub.2 to D.sub.2; g) calculating the linear attenuation coefficients at energy E.sub.0, at the points S.sub.1 and S.sub.2; h) calculating the ratio of the differential Compton cross section for a photon at energy E.sub.0 to scatter to an energy E.sub.1 to its total Compton cross section; i) calculating the ratio of the differential Compton cross section for a photon at energy E.sub.0 to scatter to an energy E.sub.2 to its total Compton cross section; j) calculating the geometrical cross section of detector D.sub.1 for the photon traveling from S.sub.1; k) calculating the geometrical cross section of detector D.sub.2 for the photon traveling from S.sub.2; l) calculating the detection efficiency of detector D.sub.1 for the photon traveling from S.sub.1 with energy E.sub.1; m) calculating the detection efficiency of detector D.sub.2 for the photon traveling from S.sub.2 with energy E.sub.2; n) calculating the image sample volumes associated with S.sub.1 and with S.sub.2; o) multiplying the quantities computed in b) through n) together to get one contribution to the double scatter count rate for the detector pair D.sub.1, D.sub.2 for each sampled time of flight offset bin in; p) selecting another pair of distinct scatter points S.sub.1 and S.sub.2, repeating steps a) through o) and adding the results to the previous results in R.sub.n.sup.12; and q) repeating step p) until the entire image volume has been adequately sampled by both S.sub.1 and S.sub.2.

3. The method of claim 2, wherein the computed line integrals, Compton cross sections, geometrical cross sections, detector efficiencies and scatter point positions are cached for reuse.

4. The method of claim 2, wherein each scatter point is chosen randomly within a cell of a regular spatial grid, one point per cell.

5. The method of claim 4, wherein a subset of all possible scatter point pairs is used for the calculation.

6. A method for applying scatter correction to a scan data acquired in a time-of-flight positron emission tomography (TOF PET) scanner, the method comprising: (A) selecting a pair of detector positions in the TOF PET scanner's detector ring; (B) numerically estimating the double scatter coincidence event rate of annihilation photons in the selected pair of detector positions, wherein one of the pair of annihilation photons in a double scatter coincidence event is assumed not to be scattered and the other of the two photons is scattered twice by a scatter point in the TOF PET scanner's image volume, during a TOF PET acquisition scan period; (C) selecting another pair of detector positions; (D) repeating step (B) wherein the pair of detector positions is the another pair of detector positions; (E) repeating steps (C) and (D) until all detector pairs of interest have been selected; and (F) scatter correcting the acquired scan data obtained from the pair of detector positions during the TOF PET acquisition scan period to reconstruct a TOF PET image based on the scan data.

7. The method of claim 6, wherein the selected pair of detector positions are D.sub.1 and D.sub.2, wherein the step (B) comprises: a) selecting the two scatter points S.sub.1 and S.sub.2 in the TOF PET scanner's image volume; b) sampling values of the emission rate density between D.sub.1 and S.sub.1, determining the time of flight offset bin, n, for each sample, and storing these values in a vector quantity R.sub.n.sup.11; c) calculating an exponential function of the negative of the line integral of the linear attenuation coefficients at energy E.sub.0 between D.sub.1 and S.sub.1; d) calculating the exponential function of the negative of the line integral of the linear attenuation coefficients at energy E.sub.1 between S.sub.1 and S.sub.2; e) calculating the exponential function of the negative of the line integral of the linear attenuation coefficients at energy E.sub.12 between S.sub.2 and D.sub.2; f) calculating the reciprocals of the squares of the distances from S.sub.1 to S.sub.2, from S.sub.1 to D.sub.1, and from S.sub.2 to D.sub.2; g) calculating the linear attenuation coefficients for energy E.sub.0 at S.sub.1 and for E.sub.1 at S.sub.2; h) calculating the ratio of the differential Compton cross section for a photon at energy E.sub.0 to scatter to an energy E.sub.1 to its total Compton cross section; i) calculating the ratio of the differential Compton cross section for a photon at energy E.sub.1 to scatter to an energy E.sub.12 to its total Compton cross section; j) calculating the geometrical cross section of detector D.sub.1 for a photon pair produced between D.sub.1 and S.sub.1; k) calculating the geometrical cross section of detector D.sub.2 for the photon traveling from S.sub.2; l) calculating the detection efficiency of detector D.sub.1 for the photon traveling from the emission point with energy E.sub.0; m) calculating the detection efficiency of detector D.sub.2 for the photon traveling from S.sub.2 with energy E.sub.12; n) calculating the image sample volumes associated with S.sub.1 and with S.sub.2; o) multiplying the quantities computed in b) through n) together to get one contribution to the double scatter count rate for the detector pair D.sub.1, D.sub.2 for each sampled time of flight offset bin in R.sub.n.sup.11; p) selecting another pair of distinct scatter points S.sub.1 and S.sub.2, repeating steps a) through o) and adding the results to the previous results in R.sub.n.sup.11; and q) repeating step p) until the entire image volume has been adequately sampled by both S.sub.1 and S.sub.2.

8. The method of claim 7, wherein the computed line integrals, Compton cross sections, geometrical cross sections, detector efficiencies, and scatter point positions are cached for reuse.

9. The method of claim 7, wherein each scatter point is chosen randomly within a cell of a regular spatial grid, one point per cell.

10. The method of claim 9, wherein a subset of all possible scatter point pairs is used for the calculation.

11. The method of claim 7, wherein the steps a) through Q) are performed with the following substitutions of the variables: D.sub.1 .Math.D.sub.2, S.sub.1 .Math.S.sub.2, E.sub.1.fwdarw.E.sub.2, and R.sub.n.sup.11.fwdarw.R.sub.n.sup.22.

12. The method of claim 11, wherein the computed line integrals, Compton cross sections, geographical cross sections, detector efficiencies and scatter point positions are cached for reuse.

13. The method of claim 11, wherein each scatter point is chosen randomly within a cell of a regular spatial grid, one point per cell.

14. The method of claim 13, wherein a subset of all possible scatter point pairs is used for the calculation.

15. A system for processing and reconstructing time-of-flight positron emission tomography (TOF PET) sinogram data, comprising: a processor capable of executing instructions; and a non-transitory, machine readable storage medium encoded with program instructions for controlling a TOF PET scanner, such that when the processor executes the program instructions, the processor performs a method for applying scatter correction to a scan data acquired in the TOF PET scanner, the method comprising: (A) selecting a pair of detector positions in the TOF PET scanner's detector ring; (B) numerically estimating the double scatter coincidence event rate of annihilation photons in the selected pair of detector positions, wherein the annihilation photons in a double scatter coincidence event are assumed to be scattered once, each photon by a different scatter point in the TOF PET scanner's image volume, during a TOF PET acquisition scan period; and (C) selecting another pair of detector positions; (D) repeating step (B) wherein the selected pair of detector positions is the another pair of detector positions; (E) repeating steps (C) and (D) until all detector pairs of interest have been selected; and (F) scatter correcting the acquired scan data obtained from the pair of detector positions during the TOF PET acquisition scan period to reconstruct a TOF PET image based on the scan data.

16. A non-transitory, machine readable storage medium encoded with program instructions for controlling a time-of-flight positron emission tomography (TOF PET) scanner, such that when a processor executes the program instructions, the processor performs a method for applying scatter correction to a scan data acquired in the TOF PET scanner, the method comprising: (A) selecting a pair of detector positions in the TOF PET scanner's detector ring; (B) numerically estimating the double scatter coincidence event rate of annihilation photons in the selected pair of detector positions, wherein the annihilation photons in a double scatter coincidence event are assumed to be scattered once, each photon by a different scatter point in the TOF PET scanner's image volume, during a TOF PET acquisition scan period; and (C) selecting another pair of detector positions; (D) repeating step (B) wherein the selected pair of detector positions is the another pair of detector positions; (E) repeating steps (C) and (D) until all detector pairs of interest have been selected; and (F) scatter correcting the acquired scan data obtained from the pair of detector positions during the TOF PET acquisition scan period to reconstruct a TOF PET image based on the scan data.

17. A system for processing and reconstructing time-of-flight positron emission tomography (TOF PET) sinogram data, comprising: a processor capable of executing instructions; and a non-transitory, machine readable storage medium encoded with program instructions for controlling a TOF PET scanner, such that when the processor executes the program instructions, the processor performs a method for applying scatter correction to a scan data acquired in the TOF PET scanner, the method comprising: (A) selecting a pair of detector positions in the TOF PET scanner's detector ring; (B) numerically estimating the double scatter coincidence event rate of annihilation photons in the selected pair of detector positions, wherein one of the pair of annihilation photons in a double scatter coincidence event is assumed not to be scattered and the other of the two photons is scattered twice by a scatter point in the TOF PET scanner's image volume, during a TOF PET acquisition scan period; (C) selecting another pair of detector positions; (D) repeating step (B) wherein the pair of detector positions is the another pair of detector positions; (E) repeating steps (C) and (D) until all detector pairs of interest have been selected; and (F) scatter correcting the acquired scan data obtained from the pair of detector positions during the TOF PET acquisition scan period to reconstruct a TOF PET image based on the scan data.

18. A non-transitory, machine readable storage medium encoded with program instructions for controlling a time-of-flight positron emission tomography (TOF PET) scanner, such that when a processor executes the program instructions, the processor performs a method for applying scatter correction to a scan data acquired in the TOF PET scanner, the method comprising: (A) selecting a pair of detector positions in the TOF PET scanner's detector ring; (B) numerically estimating the double scatter coincidence event rate of annihilation photons in the selected pair of detector positions, wherein one of the pair of annihilation photons in a double scatter coincidence event is assumed not to be scattered and the other of the two photons is scattered twice by a scatter point in the TOF PET scanner's image volume, during a TOF PET acquisition scan period; (C) selecting another pair of detector positions; (D) repeating step (B) wherein the pair of detector positions is the another pair of detector positions; (E) repeating steps (C) and (D) until all detector pairs of interest have been selected; and (F) scatter correcting the acquired scan data obtained from the pair of detector positions during the TOF PET acquisition scan period to reconstruct a TOF PET image based on the scan data.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

(1) The invention will now be more fully described by way of example with reference to the accompanying drawings in which:

(2) FIG. 1 is a diagram illustrating the relationship between PET projection data and a sinogram;

(3) FIG. 2 is a diagram illustrating the concept of time of flight in PET imaging;

(4) FIG. 3 is a diagram illustrating an example of a trajectory for a single scattering coincidence event.

(5) FIG. 4 is a diagram illustrating a possible trajectory for double scattering coincidence event with the three possible double scatter contributions to LOR (D.sub.1, D.sub.2) depending on the location of the emission source on the trajectory.

(6) FIG. 5 is a diagram illustrating a possible first scatter point S.sub.1 and three possible second scatter points, S.sub.2, S′.sub.2 and S″.sub.2, defining 3 possible double scatter pair samples.

(7) FIG. 6 is a diagram illustrating a supercell grid for N=16 and M=4.

(8) FIG. 7 is a diagram illustrating the 8 double scatter pairs associated with one supercell.

(9) FIG. 8 is a plot of average sinogram radial profiles comparing validation data against Monte Carlo simulations for a 20 cm cylinder phantom.

(10) FIG. 9 is a plot of average sinogram radial profiles comparing validation data against Monte Carlo simulations for a 30 cm cylinder phantom.

(11) FIGS. 10A-10D are coronal images from four different reconstructions of clinical PET data using the radiotracer .sup.18F-DCFPyL. ‘SSS’=single scatter simulation, ‘DSS’=double scatter simulation, ‘Rel’=relative scaling, ‘Abs’=absolute scaling. ‘Relative total activity’ refers to the total activity in the imaged portion of the body compared to the ‘SSS+DSS Abs’ case.

(12) FIGS. 11A-11C are coronal images from three different reconstructions of clinical PET data using the radiotracer .sup.18F-DCFPyL (same patient as in FIGS. 10A-10D). ‘SF’=scatter fraction, ‘it’=number of iterations used in the reconstruction, ‘sec’=seconds required for the reconstruction, ‘dbl_mult’=multiplicity factor used for selecting double scatter points.

(13) FIGS. 12A-12D are coronal images from four different reconstructions of a second clinical PET study using the radiotracer .sup.18F-DCFPyL. The notations have the same meaning as in FIGS. 10A-10D.

(14) FIGS. 13A-13C are coronal images from three different reconstructions of a clinical PET study using the radiotracer .sup.18F-FDG. The notations have the same meaning as in FIGS. 10A-10D. ‘M’=the multiplicity factor used for selecting double scatter points.

(15) FIGS. 14A-14D are coronal images from two different reconstructions for each of two clinical PET studies using the radiotracer .sup.18F-FDG. The notations have the same meaning as in FIGS. 10A-10D.

(16) FIGS. 15A-15D are coronal images from two different reconstructions for each of two clinical PET studies using the radiotracer .sup.68Ga-PSMA. The notations have the same meaning as in FIGS. 10A-10D.

(17) FIGS. 16A-16B are coronal images from two different reconstructions of a clinical PET study using the radiotracer .sup.68Ga-PSMA. The notations have the same meaning as in FIGS. 10A-10D.

(18) FIG. 17 is a schematic illustration of a transaxial slice of the simulated phantom.

(19) FIG. 18A shows a comparison of Monte Carlo simulation double scatter TOF sinograms to the corresponding DSS TOF sinograms for 8 TOF bins with significant data.

(20) FIG. 18B show a comparison of Monte Carlo simulation total scatter (all orders) TOF sinograms to the corresponding SSS+DSS TOF sinograms.

(21) FIG. 19 is a high-level block diagram of a PET scanner system.

DETAILED DESCRIPTION

(22) This description of the exemplary embodiments is intended to be read in connection with the accompanying drawings, which are to be considered part of the entire written description. The disclosed embodiments are merely exemplary of the invention and the invention may be embodied in various and alternative forms.

(23) Disclosed herein is a new double scatter simulation (DSS) algorithm that, when combined with the previously known single scatter simulation (SSS) technique, accurately approximates the total Compton scattering contribution to PET data for both TOF and non-TOF PET scanning. These double scatter contributions are efficiently computed by considering a subset of pairs of the single scatter points. By fully accounting for the physics, an absolute scaling is achieved, reducing or eliminating the need for a sometimes problematic scaling to the emission data.

(24) The DSS algorithm disclosed herein is predicated on the prior existence of an initial PET emission image (positron emissions/volume/time) of the person or object being scanned in the PET camera, and the corresponding image of the linear attenuation coefficients for the annihilation radiation within that person or object. This initial emission image may be uncorrected for scattered radiation, and thus inaccurate. The contribution of scattered radiation to this image, estimated by means of the SSS+DSS simulations, can be used to more accurately model the emission data, and thereby allow a more accurate reconstruction of the emission image. This process of scatter estimation followed by a corrected reconstruction can be iterated until the scattered radiation is fully accounted for.

(25) FIG. 4 shows a schematic illustration of a possible trajectory for the gamma photons from a double scattering coincidence event. D.sub.1 and D.sub.2 are detectors on a detector ring. S.sub.1 and S.sub.2 are possible scatter points. There are three distinct contributions, depending on the location of the emission point (the gamma annihilation point) relative to the two scatter points. For emission points on the rays (S.sub.1, D.sub.1) and (S.sub.2, D.sub.2) one photon will be unscattered and the other will be scattered twice. For emission points on the ray (S.sub.1, S.sub.2) both photons will be scattered once. According to an embodiment of the present disclosure, the total contribution from double scattering to the LOR (D.sub.1, D.sub.2) in a TOF offset bin n, from emission points along the ray R.sup.11 between S.sub.1 and D.sub.1, is estimated as in the following equation (1):

(26) $\begin{array}{cc}\begin{array}{c}{R}_n^{11}=\underset{V_S}{\int }\underset{V_S}{\int }{{\left[{\int }_{D_1}^{S_1}{h}_n\left[\Delta t\left(s\right)\right]n_e\left(s\right)\mathrm{ds}\right]}_1\left[\frac{{\sigma }_{D_1}{\mathrm{dV}}_{S_1}^{2/3}}{2\pi R_{S_1{D}_1}^2}\right]}_2\\ {{\left[e^{{\int }_{D_1}^{S_1}\mu \left(E_{0,s}\right)\mathrm{ds}}\right]}_3\left[\mu \left(E_0,S_1\right){\mathrm{dV}}_{S_1}^{1/3}\right]}_4\\ {{\left[\frac{1}{{\sigma }_c\left(E_0\right)}\frac{d{\sigma }_c\left(E_0,E_1\right)}{d\Omega }\right]}_5\left[\frac{{\mathrm{dV}}_{S_2}^{2/3}}{R_{S_1{S}_2}^2}\right]}_6\\ {{\left[e^{-{\int }_{S_1}^{S_2}\mu \left(E_1,s\right)\mathrm{ds}}\right]}_7\left[\mu \left(E_1,S_2\right){\mathrm{dV}}_{S_2}^{1/3}\right]}_8\\ {{\left[\frac{1}{{\sigma }_c\left(E_1\right)}\frac{d{\sigma }_c\left(E_1,E_{12}\right)}{d\Omega }\right]}_9\left[\frac{{\sigma }_{D_2}}{R_{S_2{D}_2}^2}\right]}_{10}\\ {{\left[e^{-{\int }_{S_2}^{D_2}\mu \left(E_{12},s\right)\mathrm{ds}}\right]}_{11}\left[{ϵ}_{D_1}\left(E_0\right){ϵ}_{D_2}\left(E_{12}\right)\right]}_{12}\end{array}& \left(1\right)\end{array}$

(27) Only Compton scattering from free electrons, with total cross-section σ.sub.c and differential scattering cross-section dσ.sub.c/dΩ, is modeled. E.sub.0 is the initial annihilation photon energy (511 keV). If θ.sub.i is the scattering angle at S.sub.i, the energy of a singly scattered photon is E.sub.i=E.sub.0/(2−cos θ.sub.i) and the energy of a doubly scattered photon is E.sub.12=E.sub.0/(3−cos θ.sub.1−cos θ.sub.2). E.sub.12 does not depend on the order of the two scatter events. There is a double volume integral over the two scatter points, and dV.sub.S is the incremental scattering volume. n.sub.e(s) is the emission rate density (annihilations/volume/time) at the spatial position s and μ(E, s) is the linear attenuation coefficient for energy E at s. R.sub.SD and R.sub.SS are the detector-scatter and scatter-scatter ray lengths, respectively. ϵ.sub.D(E) is the detection efficiency at energy E, and σ.sub.D is the geometrical cross section of the detector. Δt(s) is the coincidence detection time offset for two photons emitted from the points, which follow the double scatter trajectory. It is equal to their path length difference divided by the speed of light. h.sub.n[ ] is an index function whose value is 1 if Δt(s) falls in the n.sup.th discrete TOF offset bin, and 0 otherwise. The time resolution of the detectors is not modeled at this stage, but is applied after all contributions to the LOR (D.sub.1, D.sub.2) have been computed. The number and widths of the discrete TOF offset bins may be chosen so that they fully cover the total coincidence window τ. For the non-TOF case, there is only a single time bin to consider, which is the total coincidence window τ, and the function h.sub.n[ ] is always 1, so the time offsets Δt(s) do not need to be computed.

(28) Referring to the numbered brackets in equation (1), each numbered bracket being referred to as a term in the following discussion, Term 1 is the integral of emission density along (D.sub.1, S.sub.1). Term 2 is the geometrical cross-section to intersect D.sub.1 and S.sub.1. Term 3 is the attenuation along (D.sub.1, S.sub.1). Term 4 is the probability of interacting at the location S.sub.1. Term 5 is the probability of Compton scatter toward S.sub.2. Term 6 is the solid angle for scattering toward the location S.sub.2. Term 7 is the attenuation along (S.sub.1, S.sub.2). Term 8 is the probability of interacting at the location S.sub.2. Term 9 is the probability of Compton scatter toward D.sub.2. Term 10 is the solid angle for scattering toward D.sub.2. Term 11 is the attenuation along (S.sub.1, D.sub.2) or (S.sub.2, D.sub.2). Term 12 is the detector efficiencies.

(29) Referring to FIG. 4, in order to properly estimate the total double scatter contributions, we must also account for the other two distinct contribution sources mentioned above. The estimation of the contribution from the emission points on the ray R.sup.22 between S.sub.2 and D.sub.2 is provided by the following equation (2):

(30) $\begin{array}{cc}\begin{array}{c}{R}_n^{22}=\underset{V_S}{\int }\underset{V_{\mathrm{.Math.}}}{\int }{{\left[{\int }_{D_2}^{S_2}{h}_n\left[\Delta t\left(s\right)\right]n_e\left(s\right)\mathrm{ds}\right]}_1\left[\frac{{\sigma }_{D_2}{\mathrm{dV}}_{S_2}^{2/3}}{2\pi R_{S_2{D}_2}^2}\right]}_2\\ {{\left[e^{-{\int }_{D_2}^{S_2}\mu \left(E_{0,s}\right)\mathrm{ds}}\right]}_3\left[\mu \left(E_0,S_2\right){\mathrm{dV}}_{S_2}^{1/3}\right]}_4\\ {{\left[\frac{1}{{\sigma }_c\left(E_0\right)}\frac{d{\sigma }_c\left(E_0,E_2\right)}{d\Omega }\right]}_5\left[\frac{{\mathrm{dV}}_{S_1}^{2/3}}{R_{S_1{S}_2}^2}\right]}_6\\ {{\left[e^{-{\int }_{S_2}^{S_1}\mu \left(E_2,s\right)\mathrm{ds}}\right]}_7\left[\mu \left(E_2,S_1\right){\mathrm{dV}}_{S_1}^{1/3}\right]}_8\\ {{\left[\frac{1}{{\sigma }_c\left(E_2\right)}\frac{d{\sigma }_c\left(E_2,E_{12}\right)}{d\Omega }\right]}_9\left[\frac{{\sigma }_{D_1}}{R_{S_1{D}_1}^2}\right]}_{10}\\ {{\left[e^{-{\int }_{S_1}^{D_1}\mu \left(E_{12},s\right)\mathrm{ds}}\right]}_{11}\left[{ϵ}_{D_2}\left(E_0\right){ϵ}_{D_1}\left(E_{12}\right)\right]}_{12}\end{array}& \left(2\right)\end{array}$

(31) The estimation of the contribution from the emission points on the ray R.sup.12 between S.sub.1 and S.sub.2, where the photons will each be scattered once, is provided by the following equation (3):

(32) $\begin{array}{cc}\begin{array}{c}{R}_n^{12}=\underset{V_S}{\int }\underset{V_{\mathrm{.Math.}}}{\int }{{\left[{\int }_{D_1}^{S_2}{h}_n\left[\Delta t\left(s\right)\right]n_e\left(s\right)\mathrm{ds}\right]}_1\left[\frac{{\mathrm{dV}}_{S_1}^{2/3}{\mathrm{dV}}_{S_2}^{2/3}}{2\pi R_{S_1{S}_2}^2}\right]}_2\\ {{\left[e^{-{\int }_{S_1}^{S_2}\mu \left(E_{0,s}\right)\mathrm{ds}}\right]}_3\left[\mu \left(E_0,S_1\right){\mathrm{dV}}_{S_1}^{1/3}\right]}_4\\ {{\left[\frac{1}{{\sigma }_c\left(E_0\right)}\frac{d{\sigma }_c\left(E_0,E_1\right)}{d\Omega }\right]}_5\left[\frac{{\sigma }_{D_1}}{R_{S_1{D}_1}^2}\right]}_6\\ {{\left[e^{-{\int }_{S_1}^{D_1}\mu \left(E_1,s\right)\mathrm{ds}}\right]}_7\left[\mu \left(E_0,S_2\right){\mathrm{dV}}_{S_2}^{1/3}\right]}_8\\ {{\left[\frac{1}{{\sigma }_c\left(E_0\right)}\frac{d{\sigma }_c\left(E_0,E_2\right)}{d\Omega }\right]}_9\left[\frac{{\sigma }_{D_2}}{R_{S_2{D}_2}^2}\right]}_{10}\\ {{\left[e^{-{\int }_{S_2}^{D_2}\mu \left(E_2,s\right)\mathrm{ds}}\right]}_{11}\left[{ϵ}_{D_1}\left(E_1\right){ϵ}_{D_2}\left(E_2\right)\right]}_{12}\end{array}& \left(3\right)\end{array}$

(33) The total double scatter contribution to the LOR (D.sub.1, D.sub.2) in a TOF offset bin n is obtained by summing together the partial estimates above: R.sub.n=R.sub.n.sup.11+R.sub.n.sup.22+R.sub.n.sup.12, which has units of (number of counts)/cm.sup.2/sec. R.sub.n is the ideal TOF response. To account for the actual time resolution of the detectors, a discrete time resolution function, p.sub.mn, is applied. p.sub.mn is the probability that an event actually occurring in time bin n would be measured by the detectors in time bin m. The measured TOF response in the LOR is therefore modeled as R.sub.m=Σ.sub.n p.sub.mnR.sub.n. The ideal and measured time bins may have different structures. For example, in the ideal case there may be many small bins, while the measured data may have fewer, wider time bins.

(34) The integrals in equations (1)-(3) are computed by sampling scatter point positions, but the 1/R.sup.2.sub.SS factors in these equations suggest the contributions could become arbitrarily large when the two scatter points get close together, making it difficult to sample them adequately. This is dealt with by using an approximate analytical integration over a small volume ΔV.sub.S around S.sub.1 with a specified radius R.sub.S when R.sub.SS<R.sub.S. The possibility that R.sub.SS=0 is allowed, in which case the scatter angle is chosen randomly. When the emission points lie between S.sub.1 and S.sub.2, there are three terms in equation (3) that involve the distance between S.sub.1 and S.sub.2, as shown in the first line in the set of equations (4) below. One needs to estimate the average value of these terms over all possible scatter point pairs lying in the sphere around S.sub.1 with radius R.sub.S. Instead of estimating this by discrete sampling, it is computed analytically using certain approximations. The second line follows from the first under the approximation that the emitter density n.sub.e and the linear attenuation coefficient μ do not vary significantly between S.sub.1 and S.sub.2. The third line is equivalent to the second, and the fourth line is equivalent to the third. The fifth line follows from the fourth under the assumption that the photon attenuation over a distance R.sub.S is small, and the desired average contribution is found to be 3n.sub.e/R.sub.s. This value is used in place of sampled values of these terms whenever R.sub.SS<R.sub.S. In the TOF case, the TOF bins for the contributions are based on the positions of the sampled emission points as usual.

(35) $\left[{\int }_{S_1}^{S_2}{n}_e\left(s\right)\mathrm{ds}\right]\left[\frac{1}{R_{S_1{S}_2}^2}\right]\left[e^{-{\int }_{S_1}^{S_2}\mu \left(E_0,s\right)\mathrm{ds}}\right]\approx \frac{1}{\Delta V_s}\stackrel{R_s}{\underset{0}{\int }}\stackrel{\pi }{\underset{0}{\int }}\stackrel{2_{\pi }}{\underset{0}{\int }}\left(n_{e}R\right)\left(\frac{1}{R^2}\right)\left(e^{-\mu R}\right)R^2\sin \theta \mathrm{dR}d\theta d\varphi =\frac{4\pi n_e}{\Delta V_s}{\int }_0^{R_s}{\mathrm{Re}}^{-\mu R}\mathrm{dR}=\frac{4\pi n_e}{\Delta V_s}\frac{1}{{\mu \left(E_0\right)}^2}\left[1-e^{-\mu \left(E_0\right)R_s}\left(1+\mu \left(E_0\right)R_s\right)\right]\approx \frac{4\pi n_e{R}_s^2}{\Delta V_s}=\frac{3n_e}{R_s}\mathrm{for}\mu R_s⪡1$

(36) When the emission points lie between D.sub.1 and S.sub.1 or D.sub.2 and S.sub.2, there are only two terms in equations (1) or (2) that involve the distance between S.sub.1 and S.sub.2, as shown in the first line in the set of equations (5) below. Following the same process as above, the desired average contribution is found to be 3/R.sub.s.sup.2. This value is used in place of sampled values of these terms whenever R.sub.SS<R.sub.S.

(37) $\begin{array}{cc}\left[\frac{1}{R_{S_1{S}_2}^2}\right]\left[e^{-{\int }_{S_1}^{S_2}\mu \left(E_1,s\right)\mathrm{ds}}\right]\approx \frac{1}{\Delta V_s}\stackrel{R_s}{\underset{0}{\int }}\stackrel{\pi }{\underset{0}{\int }}\stackrel{2_{\pi }}{\underset{0}{\int }}\left(\frac{1}{R^2}\right)\left(e^{-\mu R}\right)R^2\sin \theta \mathrm{dR}d\theta d\varphi =\frac{4\pi }{\Delta V_s}{\int }_0^{R_s}{e}^{-\mu R}\mathrm{dR}=\frac{4\pi }{\Delta V_s}\frac{1}{\mu \left(E_1\right)}\left[1-e^{-\mu \left(E_1\right)R_s}\right)\approx \frac{4\pi R_s}{\Delta V_s}=\frac{3}{R_s^2}\mathrm{for}\mu R_s⪡1& \left(5\right)\end{array}$

(38) The integrals in equations (1)-(3) are computed by discrete sampling techniques. In particular, the volume integrals over V.sub.S1 and V.sub.S2 are discretized so that each volume element dV.sub.S is associated with one scatter point. For efficiency, the same set of scatter points used for single scatter estimation may be used for both double scatter points. In principle, for each dV.sub.S1, all dV.sub.S2 should be sampled. Therefore, for N scatter points, one should in principle compute contributions from all N.sup.2 ordered pairs of scatter points (allowing for the possibility of 2 scatters within one dV.sub.S cell). In practice, however, complete sampling may not be practical or necessary. The pairs only need to adequately sample the different materials present, and scattering angles. Referring to FIG. 5, it can be seen that for a first scattering point S.sub.1 some of the possible second scattering points, for example S.sub.2 and S′.sub.2, will give very similar contributions, while others, for example S.sub.2 and S″.sub.2, will give very different contributions. It is therefore advantageous for the efficiency of the simulation to choose a subset of the possible scatter point pairs that reduces the redundancy of the sampling without introducing significant bias.

(39) To this end, a reduced sampling algorithm is used: Suppose there are N scatter points. Each scatter point has a randomly chosen position within a cell of a defined regular three-dimensional grid covering the imaging volume of interest. There is one scatter point per cell. Let M (1≤M≤N) be a pair multiplicity such that each scatter point is to be associated with M other points to form M pairs per original point. The total number of scatter point pairs to be defined is thus NM. Define a super-grid having approximately √{square root over (NM)} super-cells. Each super-cell contains approximately N/√{square root over (NM)}=√{square root over (N/M)} scatter points per grid cell. Each point belongs to one and only one cell. These √{square root over (NM)} cells are chosen to be contiguous and grouped so as to have minimal spatial extent. There are approximately M×√{square root over (N/M)}=√{square root over (NM)} pairs originating in each of the √{square root over (NM)} super-cells. For each super-cell, the second points per pair are chosen randomly subject to the constraint that all super-cells must be sampled, giving exhaustive sampling at the super-cell level. The scatter contributions computed using this reduced sampling are scaled up to account for the under-sampling. When M.fwdarw.N, the super-cell is equal to a single scatter cell, and the sampling is exhaustive over all possible N.sup.2 ordered pairs.

(40) FIG. 6 shows an example of the scatter grid (dotted and solid lines) and super-grid (solid lines only), for the case N=16, and M=4. There are 16 scatter point cells and 8 super-cells, each containing 2 points. In FIG. 15, one of the super-cells is shaded, and 8 possible scatter point pairs associated with the two points in the super-cell are indicated. Out of the 256 possible scatter pairs in this example, only 64 will be sampled, but these are distributed with good spatial uniformity.

(41) The DSS technique described herein can be efficiently implemented in conjunction with the SSS technique since they can share much of the same infrastructure, such as the scatter point distribution, detector sampling, detector to scatter point distances and angles, emission and attenuation ray sums, detector efficiencies, the time-of-flight resolution model, photon cross-section tables, and other components. To these common components, the DSS technique adds the scatter pair definition, ray sums between the scatter points, cross-sections for the second scatter, analytical approximations for contributions from closely spaced pairs, an efficient technique for sampling scatter pairs, and other components. Three new parameters are introduced to control the behavior of the DSS: a minimum energy for tracking a doubly scattered photon (e.g., 350 keV), the pair multiplicity M (e.g., 100), and the scatter point separation below which the analytical approximation is used (e.g., 2.5 cm).

(42) DSS (like SSS) uses an “absolute scaling” technique to make its estimates of scatter quantitatively compatible with measured PET data. “Absolute scaling” means that a theoretically complete physical model for estimating scatter events/sec in a detector pair is used, as described in equations (1)-(3). But using this to directly estimate measured data would require a very accurate model of the detector efficiencies and count losses. However, the corresponding expression for unscattered events (below) involves similar efficiencies (last two terms):

(43) $T=\left[{\int }_{D_1}^{D_2}{n}_e\left(s\right)\mathrm{ds}\right]\left[e^{-{\int }_{D_1}^{D_2}\mu \left(E_0,s\right)\mathrm{ds}}\right]\left[\frac{{\sigma }_{D_1{D}_2}^2}{2\pi R_{D_1{D}_2}^2}\right]{ϵ}_{D_1{D}_2}^2\left(E_0\right)$

(44) The ratio of these efficiencies can be estimated much more accurately than the individual efficiencies, so these true coincidence efficiencies are estimated using the same physics models used for the scattered event efficiency calculations, and include in the SSS and DSS estimates:

(45) $e^{-{\int }_{D_1}^{D_2}\mu \left(E_0,s\right)\mathrm{ds}}{\int }_{D_1}^{D_2}{n}_c\left(s\right)\mathrm{ds}\sim {\left(\left[\frac{{\sigma }_{D_1{D}_2}^2}{2\pi R_{D_1{D}_2}^2}\right]{ϵ}_{D_1{D}_2}^2\left(E_0\right)\right)}^{-1}\left(S_1+S_2\right)=\mathrm{SSS}+\mathrm{DSS}$

(46) In other words, an estimated trues normalization is applied. This means that the sinogram computed by SSS+DSS is directly quantitatively (“absolutely”) comparable to a normalized trues sinogram that represents the attenuated forward projection of the emission-rate density, with units of events/cm.sup.2/sec.

(47) An initial validation of the DSS algorithm was performed by comparing it to Monte Carlo (MC) simulation of 20 and 30 cm diameter cylindrical water-filled phantoms, 26 cm long, uniformly activated with the positron-emitting isotope .sup.18F. These MC simulations were performed with Penelope-2011 using a cylindrical geometry approximating the Siemens mMR PET/MR scanner (detector ring diameter 65.6 cm, axial FOV 26 cm). The calculations assumed an energy resolution of 13.5% at 511 keV and a photopeak energy window of 430-610 keV. These first results are for the non-TOF case. The DSS (and SSS) simulations are absolutely scaled—there is no cross-scaling to the MC results.

(48) FIG. 8 shows profiles for the MC, SSS and DSS sinogram estimates for the 20 cm phantom. These are radial profiles averaged over azimuthal and polar angles, and over z. The data labeled ‘MC Trues’ are unscattered events, and the curve labeled ‘Model Trues’ is a simple attenuated forward projection of the phantom. The data labeled ‘MC Tot Scat’ includes all scatter orders. The curve labeled ‘SSS+DSS’ is the result of the SSS and DSS simulations combined. The data labeled ‘MC Single’ includes only single scatter events. The curve labeled ‘SSS’ is the result of the SSS simulation. The data labeled ‘MC Double’ includes only double scatter events. The curve labeled ‘DSS’ is the result of the DSS simulation. FIG. 9 shows similar profiles for the 30 cm phantom simulation.

(49) Table A reports the scatter fraction estimates resulting from the various simulations. The scatter fraction is defined as the ratio of the sum of the specified scattered events to the sum of all scattered and unscattered events from the MC simulation. Here S.sub.1, S.sub.2 and S.sub.tot are the MC single, double and total scatter fractions, respectively, and SSS and DSS refer to the single and double scatter fractions from the SSS and DSS algorithms.

(50) These results show that the SSS and DSS simulations agree very well with the MC simulations, both in shape and in the relative magnitude of double and single scatter. The scatter fractions listed in Table A confirm that a single plus double scatter simulation can account for all but a few percent of the total scatter, suggesting that absolutely scaled SSS+DSS simulations could provide a viable clinical alternative to current rescaling algorithms.

(51) TABLE-US-00001 TABLE A Scatter fractions 20 cm 30 cm MC S.sub.1 25.5% 30.9% SSS 25.1 ± .24% 30.7 ± .12% MC S.sub.2 4.05% 6.96% DSS 4.05 ± .20% 6.80 ± .17% MC S.sub.tot 30.0% 38.9% MC (S.sub.1 + S.sub.2) 29.6% 37.9% SSS + DSS 29.1% 37.4% MC (S.sub.1 + S.sub.2)/S.sub.tot 0.99 0.97 (SSS + DSS)/S.sub.tot 0.97 0.96

(52) A similar validation of DSS versus MC was performed for simulations that included time-of-flight (TOF) discrimination. In this case, the simulated phantom 100, as shown in cross-section in FIG. 17 consisted of a 30 cm diameter and 26 cm long water-filled cylinder, in which a 14 cm diameter by 26 cm long region, offset by 7.5 cm from the center, was activated with .sup.18F. The same detector ring diameter, and photopeak energy window, were used as previously, but in this case the energy resolution was set to 10%. A TOF resolution of 0.1 ns was assumed and 11 TOF bins each 0.2 ns wide were modeled.

(53) FIG. 18A compares the MC double scatter TOF sinograms to the corresponding DSS TOF sinograms, for the 8 TOF bins with significant data. In FIG. 18A, the row (A) are the MC double scatter sinograms and row (B) are the DSS TOF sinograms. FIG. 18B compares the MC total scatter (all orders) sinograms to the corresponding SSS+DSS TOF sinograms, for the 8 TOF bins with significant data. In FIG. 18B, the row (A) are the MC total scatter sinograms and row (B) are the SSS+DSS sinograms. These sinograms represent averages over polar angle and z position. Table B gives the ratio of the sums over these sinograms, as well as the ratio of their total over all TOF bins (non-TOF). The numbers in parentheses are the standard deviations of the ratios over 10 repeat SSS and DSS simulations. The good absolute agreement between the MC and SSS+DSS simulations confirms that the TOF SSS and DSS algorithms give an adequate representation of the physics of scattering in TOF PET, and further show that absolute SSS+DSS only slightly underestimates the total scatter.

(54) TABLE-US-00002 TABLE B TOP offset (ns) DSS/MC double SSS + DSS/MC total −0.4 0.873 (0.036) 0.853 (0.011) −0.2 1.008 (0.024) 0.971 (0.007) 0.0 1.014 (0.019) 0.985 (0.005) 0.2 1.045 (0.018) 1.003 (0.006) 0.4 1.038 (0.017) 0.999 (0.005) 0.6 1.027 (0.014) 0.982 (0.004) 0.8 1.018 (0.020) 0.953 (0.007) 1.0 0.729 (0.031) 0.689 (0.014) non-TOF 1.023 (0.015) 0.983 (0.005)

(55) Several tests of the DSS algorithm on clinical data were performed for the non-TOF case, to evaluate its effectiveness and robustness. FIGS. 10A-10D are coronal images from four different reconstructions of clinical PET data using the radiotracer .sup.18F-DCFPyL. They show that the absolutely scaled SSS+DSS technique results in images that are arguably superior visually and quantitatively to the other reconstruction options.

(56) FIGS. 11A-11C are zoomed image sections from three of the reconstructions shown in FIGS. 10A, 10C, and 10D, respectively, showing that the SSS+DSS Abs algorithm gives intermediate results for scatter fractions compared to the two algorithms currently available in clinical reconstruction software.

(57) FIGS. 12A-12D are coronal images from four different reconstructions of a second clinical PET study using the radiotracer .sup.18F-DCFPyL, showing that the SSS+DSS Abs images are comparable to those produced by the other reconstruction options.

(58) FIGS. 13A-13C are coronal images from three different reconstructions of a clinical PET study using the radiotracer .sup.18F-FDG. It is found that increasing the double scatter pair multiplicity factor M from 100 to 500 has very little effect on the image. This indicates the large values of M, which entail longer reconstruction times, will not be necessary.

(59) FIGS. 14A-14D are coronal images from two different reconstructions for each of two clinical PET studies using the radiotracer .sup.18F-FDG, showing that in these cases SSS+DSS Abs gives comparable results to the current standard reconstruction SSS Rel.

(60) FIGS. 15A-15D are coronal images from two different reconstructions for each of two clinical PET studies using the radiotracer .sup.68Ga-PSMA. They show that reconstructions using SSS+DSS Abs result in higher, and probably more accurate, estimates of the scatter fraction than the current clinical standard reconstruction for this type of study, SSS Abs, which is known to underestimate the total scatter.

(61) FIGS. 16A-16B are coronal images from two different reconstructions of a second clinical PET study using the radiotracer .sup.68Ga-PSMA. Again, they show that SSS+DSS Abs produces less underestimation of the total scatter than the standard SSS Abs technique, because it includes double as well as single scatter.

(62) Applying the DSS technique described above, methods for scatter correction in a TOF PET scanner will be described. It is to be understood that non-TOF DSS is included as a special case of TOF DSS.

(63) [Double Single Scatter Estimation]—

(64) In some embodiments, a method for applying scatter correction to a scan data acquired in a time-of-flight positron emission tomography (TOF PET) scanner assuming that both photons in the pairs of annihilation photons are scattered once is disclosed. The method comprises:

(65) (A) selecting a pair of detector positions in the TOF PET scanner's detector ring;

(66) (B) numerically estimating the double scatter coincidence event rate of annihilation photons in the selected pair of detector positions, wherein the annihilation photons in a double scatter coincidence event are assumed to be scattered once, each photon by a different scatter point in the TOF PET scanner's image volume, during a TOF PET acquisition scan period;

(67) (C) selecting another pair of detector positions;

(68) (D) repeating step (B) wherein the selected pair of detector positions is the another pair of detector positions;

(69) (E) repeating steps (C) and (D) until all detector pairs of interest have been selected; and

(70) (F) scatter correcting, or modeling scatter in, the acquired scan data obtained from the pair of detector positions during the TOF PET acquisition scan period to reconstruct a TOF PET image based on the scan data.

(71) Referring to equation (3), in some embodiments of the method, the selected pair of detector positions are D.sub.1 and D.sub.2, and the numerically estimating step comprises:

(72) a) selecting the two scatter points, S.sub.1 and S.sub.2, in the TOF PET scanner's image volume;

(73) b) sampling values of the emission rate density between the two scatter points S.sub.1 and S.sub.2, determining the time of flight offset bin, n, for each contribution, and storing these values in a vector R.sub.n.sup.12; [Term 1 in (3)]

(74) c) calculating an exponential function of the negative of the line integral of the linear attenuation coefficients at energy E.sub.0 between the two scatter points S.sub.1 and S.sub.2; [Term 3 in (3)]

(75) d) calculating the exponential function of the negative of the line integral of the linear attenuation coefficients at energy E.sub.1 between S.sub.1 and D.sub.1; [Term 7 in (3)]

(76) e) calculating the exponential function of the negative of the line integral of the linear attenuation coefficients at energy E.sub.2 between S.sub.2 and D.sub.z; [Term 11 in (3)]

(77) f) calculating the reciprocals of the squares of the distances from S.sub.1 to S.sub.2, from S.sub.1 to D.sub.1, and from S.sub.2 to D.sub.2; [Terms 2, 6 and 10 in (3)]

(78) g) calculating the linear attenuation coefficients at energy E.sub.0, at the points S.sub.1 and S.sub.2; [Terms 4 and 8 in (3)]

(79) h) calculating the ratio of the differential Compton cross section for a photon at energy E.sub.0 to scatter to an energy E.sub.1 to its total Compton cross section; [Term 5 in (3)]

(80) i) calculating the ratio of the differential Compton cross section for a photon at energy E.sub.0 to scatter to an energy E.sub.2 to its total Compton cross section; [Term 9 in (3)]

(81) j) calculating the geometrical cross section of detector D.sub.1 for the photon traveling from S.sub.1; [Term 6 in (3)]

(82) k) calculating the geometrical cross section of detector D.sub.2 for the photon traveling from S.sub.2; [Term 10 in (3)]

(83) l) calculating the detection efficiency of detector D.sub.1 for the photon traveling from S.sub.1 with energy E.sub.1; [Term 12 in (3)]

(84) m) calculating the detection efficiency of detector D.sub.2 for the photon traveling from S.sub.2 with energy E.sub.2; [Term 12 in (3)]

(85) n) calculating the image sample volumes associated with S.sub.1 and with S.sub.2; [Terms 2, 4 and 8 in (3)]

(86) o) multiplying the quantities computed in b) through n) together to get one contribution to the double scatter count rate for the detector pair D.sub.1, D.sub.2, for each sampled time of flight offset bin in R.sub.n.sup.12;

(87) p) selecting another pair of distinct scatter points S.sub.1 and S.sub.2, repeating steps a) through o) and adding the results to the previous results in R.sub.n.sup.12;

(88) q) repeating step p) until the entire image volume has been adequately sampled by both S.sub.1 and S.sub.2.

(89) In some embodiments of the method, the computed line integrals, Compton cross sections, geometrical cross sections, detector efficiencies, and scatter point positions can be cached for reuse.

(90) In some embodiments of the method, each scatter point is chosen randomly within a cell of a regular spatial grid, one point per cell. Here, ‘regular spatial grid’ refers to a three-dimensional rectilinear grid having uniform cell size. In some embodiments, a subset of all possible scatter point pairs is used for the calculation.

(91) [Single Double Scatter Estimation Case]—

(92) In some embodiments, a method for applying scatter correction to a scan data acquired in a time-of-flight positron emission tomography (TOF PET) scanner assuming that one of the two photons in the pairs of annihilation photons is scattered twice and the other of the two photons is not scattered is disclosed. The method comprises:

(93) (AA) selecting a pair of detector positions in the TOF PET scanner's detector ring;

(94) (BB) numerically estimating the double scatter coincidence event rate of annihilation photons in the selected pair of detector positions, wherein one of the pair of annihilation photons in a double scatter coincidence event is assumed not to be scattered and the other of the two photons is scattered twice, once at each of two scatter points in the TOF PET scanner's image volume, during a TOF PET acquisition scan period;

(95) (CC) selecting another pair of detector positions;

(96) (DD) repeating step (BB) wherein the pair of detector positions is the another pair of detector positions;

(97) (EE) repeating steps (CC) and (DD) until all detector pairs of interest have been selected; and

(98) (FF) scatter correcting, or modeling scatter in, the acquired scan data obtained from the pair of detector positions during the TOF PET acquisition scan period to reconstruct a TOF PET image based on the scan data.

(99) Referring to equation (1), in some embodiments of the method, the selected pair of detector positions are D.sub.1 and D.sub.2 and the step (BB) comprises:

(100) aa) selecting the two scatter points S.sub.1 and S.sub.2 in the TOF PET scanner's image volume;

(101) bb) sampling values of the emission rate density between D.sub.1 and S.sub.1, determining the time of flight offset bin, n, for each sample, and storing these values in a vector quantity R.sub.n.sup.11; [Term 1 in (1) or (2)]

(102) cc) calculating an exponential function of the negative of the line integral of the linear attenuation coefficients at energy E.sub.0 between D.sub.1 and S.sub.1; [Term 3 in (1) or (2)]

(103) dd) calculating the exponential function of the negative of the line integral of the linear attenuation coefficients at energy E.sub.1 between S.sub.1 and S.sub.2; [Term 7 in (1) or (2)]

(104) ee) calculating the exponential function of the negative of the line integral of the linear attenuation coefficients at energy E.sub.12 between S.sub.2 and D.sub.2; [Term 11 in (1) or (2)]

(105) ff) calculating the reciprocals of the squares of the distances from S.sub.1 to S.sub.2, from S.sub.1 to D.sub.1, and from S.sub.2 to D.sub.2; [Terms 2, 6 and 10 in (1) or (2)]

(106) gg) calculating the linear attenuation coefficients for energy E.sub.0 at S.sub.1 and for E.sub.1 at S.sub.2; [Terms 4 and 8 in (1) or (2)]

(107) hh) calculating the ratio of the differential Compton cross section for a photon at energy E.sub.0 to scatter to an energy E.sub.1 to its total Compton cross section; [Term 5 in (1) or (2)]

(108) ii) calculating the ratio of the differential Compton cross section for a photon at energy E.sub.1 to scatter to an energy E.sub.12 to its total Compton cross section; [Term 9 in (1) or (2)]

(109) jj) calculating the geometrical cross section of detector D.sub.1 for a photon pair produced between D.sub.1 and S.sub.1; [Term 6 in (1) or (2)]

(110) kk) calculating the geometrical cross section of detector D.sub.2 for the photon traveling from S.sub.2; [Term 10 in (1) or (2)]

(111) ll) calculating the detection efficiency of detector D.sub.1 for the photon traveling from the emission point with energy E.sub.0; [Term 12 in (1) or (2)]

(112) mm) calculating the detection efficiency of detector D.sub.2 for the photon traveling from S.sub.2 with energy E.sub.12; [Term 12 in (1) or (2)]

(113) nn) calculating the image sample volumes associated with S.sub.1 and with S.sub.2; [Terms 2, 4, 6 and 8 in (1) or (2)]

(114) oo) multiplying the quantities computed in bb) through nn) together to get one contribution to the double scatter count rate for the detector pair D.sub.1, D.sub.2, for each sampled time of flight offset bin in R.sub.n.sup.11;

(115) pp) selecting another pair of distinct scatter points S.sub.1 and S.sub.2, repeating steps aa) through oo) and add the results to the previous results in R.sub.n.sup.11; and

(116) qq) repeating step pp) until the entire image volume has been adequately sampled by both S.sub.1 and S.sub.2.

(117) In some embodiments of the method, the computed line integrals, Compton cross sections, geometrical cross sections, detector efficiencies, and scatter point positions can be cached for reuse. In some embodiments of the method, each scatter point is chosen randomly within a cell of a regular spatial grid, one point per cell. In some embodiments, a subset of all possible scatter point pairs is used for the calculation.

(118) Referring to equation (2) instead of (1), in some embodiments of the method comprising the steps (AA), (BB), and (CC) described above, for the case that the annihilation photons are emitted between D.sub.2 and S.sub.2, the step (BB) can comprise the steps aa) through qq) described above but with the following substitutions of the variables: D.sub.1 .Math.D.sub.2, S.sub.1 .Math.S.sub.2, E.sub.1.fwdarw.E.sub.2, and R.sub.n.sup.11.fwdarw.R.sub.n.sup.22. In some embodiments of the method, the computed line integrals, Compton cross sections, geographical cross sections, detector efficiencies and scatter point positions can be cached for reuse. In some embodiments of the method, each scatter point is chosen randomly within a cell of a regular spatial grid, one point per cell. In some embodiments of the method, a subset of all possible scatter point pairs is used for the calculation.

(119) In some embodiments of the method, the total double scatter contribution for each detector pair of interest is computed as: R.sub.n=R.sub.n.sup.11+R.sub.n.sup.22+R.sub.n.sup.12. In some embodiments of the method, a discrete detector time resolution function, p.sub.mn, is applied according to R.sub.m=Σ.sub.np.sub.mnR.sub.n, to estimate the probability that an event actually occurring in time bin n would be measured by the detectors in time bin m.

(120) According to another aspect of the present disclosure, a method for applying scatter correction to a scan data acquired in a TOF PET scanner comprises (I) calculating the time of flight offset distribution of two detected photons in double scatter coincidence events; and (II) scatter correcting acquired scan data obtained during the TOF PET acquisition scan period to reconstruct a TOF PET image based on the scan data. In some embodiments, the step (I) can comprise:

(121) a) selecting two detector positions D.sub.1 and D.sub.2 in the TOF PET scanner's detector ring, defining them so that D.sub.1 is in the positive time offset direction;

(122) b) creating a time offset array initialized to zero;

(123) c) selecting two scatter points in the image volume, thereby defining a double scatter path;

(124) d) for each point sampled in the discrete line integrals of the emission rate density along the segments of the double scatter path, computing (1) the distance d.sub.1 along the double scatter path from the emission point to D.sub.1, (2) the distance d.sub.2 along the double scatter path from the emission point to D.sub.2, and (3) the contribution to the detector response from the emission rate density at this emission point;

(125) e) computing the time of flight offset as t=(d.sub.2−d.sub.1)/c, where c is the speed of light;

(126) f) adding the contribution to the double scatter count rate in (D.sub.1, D.sub.2) from this emission point into the bin of the time offset array corresponding to t;

(127) g) repeating steps c) through f) until all scatter point pairs have been adequately sampled;

(128) h) recording the values of the time offset array in the double scatter sinogram, in the appropriate time of flight bins for the detector pair (D.sub.1, D.sub.2); and

(129) i) repeating steps a) through h) for all detector pairs represented in the sinogram. In some embodiments of the method, the data in the time offset array are convolved with a time resolution function before being recorded.

(130) According to another aspect of the present disclosure, a method for applying scatter correction to a scan data acquired in a time-of-flight positron emission tomography (TOF PET) scanner is disclosed. The method comprises estimating contributions from annihilation photons from double scatter coincidence events arising from two scatter points S.sub.1 and S.sub.2 when the two scatter points lie close together; and scatter correcting acquired scan data obtained from a pair of detector positions during the TOF PET acquisition scan period to reconstruct a TOF PET image based on the scan data.

(131) In some embodiments of the method, the estimating step comprises: defining a minimum distance R.sub.s between the two scatter points S.sub.1 and S.sub.2; and when the distance between S.sub.1 and S.sub.2 is less than R.sub.s, using an analytical calculation of the average contribution from certain terms in the expressions for the double scatter event rates. These are terms 1, 2 and 3 in equation (3); and terms 6 and 7 in equations (1) and (2).

(132) In some embodiments of the method, for the case that both photons are scattered once, the three terms 1, 2 and 3 in equation (3) that involve the distance between the scatter points are replaced by the quantity 3n.sub.e/R.sub.s, where n.sub.e is the average emitter rate density at the scatter points.

(133) In some embodiments of the method, for the case that one photon is scattered twice and the other is not scattered, the two terms 6 and 7 in equations (1) and (2) that involve the distance between the scatter points are replace by the quantity 3/R.sub.s.sup.2.

(134) According to some embodiments, a method for applying scatter correction to a scan data acquired in a TOF PET scanner is disclosed. The method comprises: selecting a set of double scatter point pairs that contribute to double scattering of annihilation photons; and scatter correcting the acquired scan data to reconstruct a TOF PET image based on the scan data.

(135) In the various embodiments of the method, the step of selecting a set of double scatter point pairs that contribute to double scattering of annihilation photons comprises:

(136) a) defining a set of N single scatter points in the reconstructed PET images from which the double scatter contribution is to be computed;

(137) b) defining a pair multiplicity factor M which is between 1 and N;

(138) c) defining a grid in the images having approximately √{square root over (NM)} cells, with approximately √{square root over (N/M)} scatter points per cell, each point belonging to one and only one cell; and

(139) d) forming NM scatter point pairs by randomly choosing scatter points subject to the constraint that the scatter point pairs sample all cell pairs exhaustively. In some embodiments of the method, the N single scatter points are uniformly distributed over the entire volume of the images. In some embodiments of the method, the N single scatter points are randomly distributed over the entire volume of the images. In some embodiments of the method, the grid is a uniform rectilinear grid with a constant cell size, covering the entire volume of the images.

(140) FIG. 19 illustrates one embodiment of a nuclear imaging system 2. The nuclear imaging system 2 includes a scanner for a PET modality 12 provided in a first gantry 16a. A patient 17 lies on a movable patient bed 18 that can be movable with respect to the first gantry 16a. The PET modality 12 includes a plurality of detectors 50 configured to detect an annihilation photons. FIGS. 3 and 4 show the detectors 50 in their ring configuration.

(141) Scan data from the PET modality 12 is stored at one or more computer databases 40 and processed by one or more computer processors 60 of an accompanying computer system 30. The graphical depiction of the computer system 30 in FIG. 19 is provided by way of illustration only, and the computer system 30 may include one or more separate computing devices. The scan data can be provided by the PET modality 12, the second modality 14, and/or may be provided as a separate data set, such as, for example, from a memory coupled to the computer system 30. The computer system 30 can include one or more processing electronics for processing a signal received from the detectors 50.

(142) The methods and system described herein can be at least partially embodied in the form of computer-implemented processes and apparatus for practicing those processes. The disclosed methods may also be at least partially embodied in the form of tangible, non-transitory machine readable storage media encoded with computer program code. The media may include, for example, RAMs, ROMs, CD-ROMs, DVD-ROMs, BD-ROMs, hard disk drives, flash memories, or any other non-transitory machine-readable storage medium, wherein, when the computer program code is loaded into and executed by a computer, the computer becomes an apparatus for practicing the method. The methods may also be at least partially embodied in the form of a computer into which computer program code is loaded and/or executed, such that, the computer becomes a special purpose computer for practicing the methods. When implemented on a general-purpose processor, the computer program code segments configure the processor to create specific logic circuits. The methods may alternatively be at least partially embodied in a digital signal processor formed of application specific integrated circuits for performing the methods.

(143) In some embodiments, at least one non-transitory computer-readable storage medium is provided having computer-executable instructions embodied thereon, wherein, when executed by at least one processor, the computer-executable instructions cause the at least one processor to perform embodiments of the methods described herein.

(144) Although the subject matter has been described in terms of exemplary embodiments, it is not limited thereto. Rather, the appended claims should be construed broadly, to include other variants and embodiments, which may be made by those skilled in the art.