Time efficient scatter correction for time-of-flight PET: the immediate scatter approximation

Utilization of time-of-flight (TOF) information allows us to improve image quality and convergence rate in iterative PET image reconstruction. In order to obtain quantitatively correct images accurate scatter correction (SC) is required that accounts for the non-uniform distribution of scatter events over the TOF bins. However, existing simplified TOF-SC algorithms frequently exhibit limited accuracy while the currently accepted reference method—the TOF extension of the single scatter simulation approach (TOF-SSS)—is computationally demanding and can substantially slow down the reconstruction. In this paper we propose and evaluate a new accelerated TOF-SC algorithm in order to improve this situation. The key idea of the algorithm is the use of an immediate scatter approximation (ISA) for scatter time distribution calculation which speeds up estimation of the required TOF scatter by a factor of up to five in comparison to TOF-SSS. The proposed approach was evaluated in dedicated phantom measurements providing challenging high activity contrast conditions as well as in representative clinical patient data sets. Our results show that ISA is a viable alternative to TOF-SSS. The reconstructed images are in excellent quantitative agreement with those obtained with TOF-SSS while overall reconstruction time can be reduced by a factor of two in whole-body studies, even when using a listmode reconstruction not optimized for speed.


Introduction
Utilization of time-of-flight (TOF) information provides clear benefits for iterative PET image reconstruction in terms of reduced image noise and improved convergence rate (Conti et al 2005, Vunckx et al 2010. Moreover, it has been shown that TOF image reconstruction is more robust in the presence of data inconsistencies such as erroneous attenuation map or detector normalization (Turkington and Wilson 2009, Conti 2011, ter Voert et al 2017, and provides better lesion detectability (Surti et al 2006, Kadrmas et al 2009, El Fakhri et al 2011, Mühlematter et al 2018. Finally, utilizing the TOF information allows us to reduce the cross-talk between reconstructed activity and attenuation maps in joint image reconstruction (Defrise et al 2012).
In order to fully utilize the benefits of TOF PET imaging and to obtain quantitatively correct images various corrections have to be applied prior to and during image reconstruction. Detector normalization as well as dead time correction are essentially the same as in non-TOF PET, and random coincidences are distributed uniformly over all TOF bins (Conti et al 2005). The TOF bin distribution of scattered events, however, is non-uniform. Therefore, conventional scatter correction (SC) methods, notably the single scatter simulation (SSS) approach (Watson 2000), cannot be used without modification. The most accurate way to estimate the scatter time distribution is to include TOF modeling directly into the scatter simulation process which leads to the TOF extension of the SSS algorithm (Werner et al 2006, Watson 2007. However, the increased complexity of the TOF-SSS algorithm causes an increase in computation time by an implementation dependent factor of about 3-7. This can result in significant image reconstruction slow-down for certain practically relevant choices of iteration scheme and reconstruction parameters. There are alternative TOF-SC approaches which allow to avoid this substantial computational overhead and to perform TOF scatter estimation only slightly slower than non-TOF algorithms. The key idea is to use a non-Time efficient scatter correction for time-of-flight PET: the immediate scatter approximation TOF scatter correction algorithm-in particular SSS-to estimate the spatial distribution of scattered events in each line of response (LOR) and to model scatter time distribution in a separate step that utilizes simplifying assumptions. The first two approaches of this kind were introduced in Conti et al (2005) and intended as a temporary solution since full TOF-SSS was still under development and not yet available: in the so-called simple scaling approach non-TOF scatter sinograms are produced first and then rescaled to fit the emission data tails individually for each TOF bin. In the radial distortion and scaling approach, the non-TOF scatter sinograms are modified prior to re-scaling using the assumption that the scatter is uniformly distributed in the scanner's field of view (FOV). As reported in Conti et al (2005), the first approach (simple scaling) significantly overestimates scatter in the center of large objects. The second approach (radial distortion and scaling) performs better in this regard but ignores dependency of the scatter TOF profile on the given activity distribution which is especially problematic under high contrast conditions (e.g. between bladder and surrounding tissue).
A further accelerated TOF-SC method was presented in Jin et al (2013) as part of the MOLAR image reconstruction for the Siemens Biograph mCT (Siemens Medical Solutions USA, Inc.). This method assumes that, both, scattered events and the sum of scattered and true events-estimated as 'prompts minus delayed'-have a similar time distribution in each LOR. It is relevant in this context that the mCT scanner uses an axial span of 11 and twofold angular mashing even in list-mode format due to intrinsic data format limitation which heavily reduces the total number of LORs and correspondingly increases the number of events in each LOR. In combination with a relatively small number of TOF bins (13) excessive noise in the 'prompts minus delayed' time distribution can thus be avoided in this specific setting. However, in the practically more relevant case of small event numbers per LOR and TOF bin this approach leads to significant scatter overcorrection by wrongly assigning the estimated scatter to just a few TOF bins. Moreover, the initial assumption of the method will not hold for LORs in the vicinity of large structures like the bladder exhibiting high activity concentrations. In this case, the scatter TOF distribution will be significantly affected by neighboring activity while the trues TOF distribution is expected to remain relatively flat.
To overcome the limitations of the available approaches, we have developed a new time efficient TOF-SC method. Our method also relies on the principal idea of separate estimation of scatter spatial distribution (via SSS) and scatter time distribution via a dedicated fast algorithm. The key difference to existing approaches is that the proposed scatter time distribution algorithm works reliably by design for all practically relevant event numbers and subdivisions into TOF bins. Moreover, the activity distribution dependence of the scatter TOF profiles was addressed.
The proposed approach was evaluated in dedicated phantom measurements providing challenging high activity contrast conditions and compared against several other reconstruction schemes including TOF-SSS (which we consider as the accepted reference method). A further evaluation was performed in representative clinical patient data sets.

Methods
Image reconstruction was performed with the Tube of Response High Resolution OSEM Reconstruction tool (THOR), previously developed by our group (Lougovski et al 2014). THOR uses the ordered subsets accelerated listmode version of maximum-likelihood estimation-maximization algorithm (MLEM), see Shepp and Vardi (1982), Hudson and Larkin (1994) and Levkovitz et al (2001). The original version of THOR as presented in Lougovski et al (2014) features a unique system matrix calculation approach and includes all corrections necessary for quantitative image reconstruction, i.e. proper normalization as well as attenuation, dead time, randoms, and (non-TOF) SSS-based scatter correction. So far, however, THOR did not utilize any TOF information. We, therefore, have extended THOR for the present investigation accordingly in order to handle and use TOF data as well.
In the TOF case, the iterative formula for listmode MLEM can be written in the form where λ (n) j denotes estimated activity concentration in voxel j after the iteration n, i enumerates the different LORs, c ik,δ is a system matrix element for LOR i and voxel k with attenuation and normalization corrections factored in. The index δ explicitly denotes the given dependency of the respective quantity on the difference of photon arrival times at both detectors. The notation i(e) and δ(e) indicates that these indices are depending on the considered specific event e. Accordingly, s i,δ and r i,δ are TOF weighted additive contributions of scattered and random events in LOR i. In the following, we use the convention that omission of index δ designates the non-TOF value of the respective quantity.
Since random events have a uniform time distribution the non-TOF random estimate r i,δ ≡ r i can be used directly in (1). Scattered events, on the other hand, have a non-uniform time distribution which has to be accounted for. In order to achieve this, we express the TOF scatter estimate s i,δ for LOR i as a product of the non-TOF scatter value s i and TOF weight factor w i,δ : For calculation of w i,δ we use a two step procedure. First, we calculate 'true' scatter time distributions m iτ which correspond to a scanner with infinite time resolution. For practical reasons, we store a discrete approximation of this distribution using a sufficiently large number N τ of TOF bins and normalize it as follows: The required w i,δ coefficients are then computed for each event individually during the MLEM update as a scaled and weighted sum of the m iτ where σ TOF is the scanner's time resolution and A δ is a normalization constant. The integrals in (4) and (5) extend over the time intervals ∆ τ corresponding to the individual time bins τ . The sum of these integrals in (5) is thus equal to the integral over all possible TOF time differences along the whole LOR. Equation (4) serves as our model to describe the influence of finite time resolution being a convolution of a truncated Gaussian kernel and the discretized true scatter TOF distribution m iτ with a proper normalization enforced by (2) and (3). Obviously, in the non-TOF limit all w i,δ ≡ 1 whereas for infinitely increasing time resolution w i,δ becomes asymptotically identical to the product of the respective m iτ and N τ . TOF-SSS estimates s i and m iτ simultaneously and serves as the reference method in this study. The newly proposed SC method, however, is combining a non-TOF SSS-based computation of the scatter distribution s i with a dedicated fast algorithm to estimate the scatter time distribution m iτ as explained in the following.

Single scatter simulation
The original SSS according to Watson (2000) was integrated into THOR and used for this study. As mentioned in the original paper, due to the smoothness of the scatter distribution it is sufficient to estimate the scatter only for a small subset of LORs and then to interpolate between these. We used an overall undersampling by a factor of 128 (2-fold angular × 8-fold radial × 8-fold axial) resulting in a corresponding acceleration of the computations. The same interpolation procedure was used for all scatter estimation procedures within this work.
SSS does not account for multiple scatter events and out-of-FOV scatter. To correct for these effects, a scaling of the SSS scatter distribution can be performed (Polycarpou et al 2011). However, accuracy of this correction depends on the chosen scaling method, e.g. variants of the tail fitting approach, maximum-likelihood scatter scaling (Rezaei et al 2017), or Monte Carlo simulation based scaling (Ye et al 2014). According to our experience, a good combination of simplicity, accuracy, and reliability can be achieved with plane-wise scatter scaling using weighting coefficients of (1/a i ) 4 for the LORs in the tails region. This method was used in the present investigation. The activity support was determined by suitable thresholding of the attenuation map excluding the patient bed. A safety margin of 10 mm was added to the derived body outline in order to avoid any contamination of the LORs used for tail fitting by residual activity contributions due to finite resolution or residual patient motion.
The out-of-FOV scatter estimate can be further improved if information regarding activity and attenuation from adjacent bed positions is available, see Iatrou et al (2009). THOR takes advantage of this fact for whole-body reconstruction by performing the scatter simulation for a 'virtual total body scanner' whose axial FOV covers all bed positions of the respective scan.
Formally, SSS yields scatter within the considered LOR after integration over all photon arrival time differences, i.e. it does not provide the probability distribution of the scatter over the different TOF bins. Nevertheless, one might consider to use SSS in a TOF reconstruction by assuming uniform distribution of the scatter over all TOF bins. This is equivalent to assuming that all TOF weights w i,δ are equal to 1. We have included this hybrid approach in our investigation in order to evaluate the importance of adequately modeling the time structure of the scatter.

Immediate scatter approximation
Our objective is to adequately model the dependency of scatter TOF profiles on activity distribution (while neglecting the additional influence of spatial variations of the attenuation coefficient). The relevant geometry is shown in figure 1. The method utilizes the fact that photon pairs emitted in point E and detected in LOR D 1 D 2 after single scatter at some point S in most cases will exhibit an arrival time difference close to (ED 1 − ED 2 ) (we are using units such that the speed of light c = 1). For illustration, figures 2 and 3 demonstrate this behavior for a 2D acquisition of a point source in the center of a cylindrical water phantom. In this example, the TOF difference should thus be close to zero for most events. The simulation underlying these plots scans a sufficiently fine grid of scatter locations and scatter angles, computing for each grid point the resulting LOR and TOF difference as well the relative scatter and detection probability in the given energy window. As can be seen in figure 2, the probability of detecting a scattered event rapidly decreases with increasing LOR distance as well as TOF difference. This behavior arises due to a combination of two factors. First, the actual arrival time difference is (ED 1 − (ES + SD 2 )) and for LORs with a sufficiently small orthogonal distance from E (corresponding to small angle forward scattering) we have ES + SD 2 ≈ ED 2 even for relatively distant scatter points S (for which ES is not small compared to ED 2 ). Second, although for more distant LORs the dependency of measured TOF information on the distance ES becomes more prominent, the relative contribution of the considered emission point to an LOR rapidly decreases with increasing distance of the LOR from E. The reason for this is lower scatter probability at large angles in combination with lower detection/acceptance probability of the scattered photon by D 2 . The latter is a consequence of the reduced energy of the scattered photon which reduces overlap of the detector signal with the given energy window of 460-665 keV for the considered scanner. Therefore, for the overwhelming number of events contributing to a given LOR the TOF difference of the registered scatter events essentially only depends on the position of the emission point relative to the LOR while the dependence on the actual scatter point position might be neglected.
To a lesser extent, approximate independence of scatter point location is also valid for the scatter angle: for most events the actual scatter angle is not much larger than for a scatter occurring immediately after the emission at E. So to a rough approximation one can compute the scatter angle for that limiting case only and neglect here, too, the dependence on actual scatter position.
Altogether, we surmise that it should be sufficient to estimate scatter profiles m iτ based on the assumption that for all scatter events only one of the photons scatters and does so immediately after emission. This is the key assumption of our new immediate scatter approximation (ISA) which we propose and evaluate in the present study. Our ISA implementation models the object as a relatively small set of emission points E with activities I E taken from the current image estimate. The locations of the emission points are chosen to coincide with those of the scatter points used for SSS according to Watson (2000).
For each LOR D 1 D 2 and emission point E the arrival time difference of any event originating from that point is estimated as δ i,E = ED 1 − ED 2 (neglecting the residual scatter point dependency as explained above). By repeating this procedure for all emission points one arrives at the scatter time distribution where θ is scattering angle (for scatter happening directly at E), [dσ c /dΩ](θ) is the differential Compton cross section at 511 keV, and ε(θ) is the energy dependent detection probability for the scattered photon. ∆Ω 1 and Figure 1. Illustration of ISA and TOF-SSS. D 1 , D 2 : detectors, S: scatter point, E: emission point. E is the corresponding apparent source position as resulting from TOF-SSS while E is the position resulting from ISA. Note that the difference between E and E is small despite considering an LOR with very large radial offset from E and scatter at a rather distant point.
∆Ω 2 are the solid angles extended by both detectors at E. N i is a normalization coefficient ensuring satisfaction of (3). Note that the quantities [dσ c /dΩ](θ), ε(θ), ∆Ω 1 , ∆Ω 2 in (6) are already required for SSS and can therefore be recycled by ISA if the emission points are chosen to coincide with the SSS scatter points.
TOF single scatter simulation TOF-SSS was implemented similar to Werner et al (2006). The LOR undersampling factor and scatter point locations were the same as those used for SSS.

Patient and phantom data
The different TOF-SC algorithms were evaluated on a set of phantom and patient scans acquired with the Ingenuity PET/MR scanner (Philips Healthcare, Best, The Netherlands) in list-mode format. This scanner has a count rate dependent time resolution of 600-700 ps, a coincidence time window of 6 ns, and a list-mode TOF bin width of 25 ps. Attenuation maps were generated from MR image segmentation into 3 classes (air, lungs, soft Figure 2. Illustration of the TOF difference versus LOR distance distribution of the scattered events produced by a point source located at the center of a cylindrical water phantom of 40 cm diameter centered in the field of view. The chosen scanner geometry and energy acquisition window are that of the Philips Ingenuity PET/MR. The 2D case is considered. The influence of energy dependent attenuation (different before and after scatter), anisotropic scatter probability, energy dependent detection probability, and position dependent solid angles are included while multiple scatter was neglected. The plot shows the relative probability of detecting a scattered event at TOF difference δ in an LOR with radial offset h. The distribution rapidly approaches zero with increasing LOR distance and TOF difference.  figure 2 along the corresponding axis. Shown in blue are the predictions resulting from ISA for this configuration (zero TOF difference for all scattered events, modest overemphasis of larger LOR distances). Note the small range of actually occuring TOF differences (60 ps corresponding to 9 mm shift in apparent event position along the LOR relative to ISA prediction).
tissue) with the default vendor provided software, see Hu et al (2009). Additional truncation compensation was performed when necessary (Schramm et al 2013).
The phantom modeled extreme high contrast conditions in the pelvic region in order to evaluate scatter correction in a worst case scenario. We used the whole-body phantom L981602 (PTW-Freiburg, Freiburg, Germany) with two spherical inserts. The bigger sphere (R = 36 mm) is located in the center and simulates a hot bladder. The second sphere (R = 13.8 mm) is located 5 cm away from the other one to describe a tumor lesion near the bladder. We chose a 'bladder:lesion:background' activity concentration ratio of 40:5:1 which is comparable to concentration ratios frequently observed in the pelvic region in clinical 18 F-FDG PET scans. A homogeneous attenuation map with the attenuation coefficient set to that of water was used since this is a good approximation for this setup of a water-filled phantom with water-filled inserts.
Accuracy of the SC algorithms was quantitatively assessed by comparing image-derived activity ratios of relevant regions with their true values. Of special interest in this context is the concentric neighbourhood of the bladder since it is most strongly affected by inaccuracies of the scatter correction. A frequently observed problem in this region is severe scatter overcorrection (known as photopenic or halo artifact) which thus provides a critical test of SC accuracy. The activity in the background was determined from a set of 14 cylindrical regions-ofinterest (ROI) with a total volume of 805 cm 3 equally spaced along the rim approximately 3 cm inward from the edge of the phantom. Activities in the bladder and lesion were measured in spherical ROIs of 59 cm 3 and 2 cm 3 , respectively, centered within the respective region ( figure 5).
Performance of the different TOF-SC algorithms was also assessed in typical whole-body and brain scans acquired on the Philips Ingenuity PET/MR scanner at the University Hospital Carl Gustav Carus, Dresden. Informed consent was obtained from all subjects. Results for two patients are presented in the following. The first patient had a cervical lymph node metastases of a squamous-cell carcinoma. He received 333 MBq dose of 18 F-FDG and underwent 10 bed positions whole-body PET scan starting 65 min post-injection with 2 min scan time per bed position. The second patient exhibited a grey matter heterotopia in the left mesial temporal lobe. A brain scan of 20 min was performed starting 30 min after injection of 210 MBq of 18 F-FDG.
The data were reconstructed with THOR in TOF and non-TOF mode. Reconstruction settings are listed in table 1. 25 TOF bins yield a bin width of 240 ps and 50 TOF bins yield a bin width of 120 ps. The same iteration schemes were used for both TOF and non-TOF reconstruction. No post-smoothing was applied.
Initial scatter estimation was performed in TOF mode with respective algorithm and the scatter time distributions were derived (i.e. the time demanding procedure of scatter time distribution calculation was performed only once per reconstruction). Non-TOF SSS was used for subsequent scatter updates during the iterative reconstruction process. A simplified scheme of the reconstruction process is shown in figure 4. Two scatter updates were done for the whole-body scans while three scatter updates were done for the brain and phantom scans since these were reconstructed with a protocol using an increased number of iterations.

Run time performance
All THOR reconstructions were performed on a cluster of ten servers (8 × 28 cores Intel Xeon E5-2690 v4 + 1 × 24 cores Intel Xeon E5-2690 v3 + 1 × 24 cores Intel Xeon E5-2697 v2, 128 GB RAM each), distributing the computational load according to the given performance differences of the available machines. Reconstruction times were averaged over 10 runs.

Results
Reconstructed phantom images are shown in figure 6. The corresponding quantitative evaluation is presented in table 2. The presented results are accompanied by sample phantom image profiles comparison given in figure 7. The vendor (TOF BLOB-OS-TF) reconstruction cannot handle this extreme high contrast case. Scatter is heavily overcorrected near the 'bladder', resulting in a massive halo artifact with a spurious 65.8% signal drop. Activity in the background remote from the central sphere approaches its true value (=1 in the units used in figure 6) only at the very edge of the phantom.
In comparison, non-TOF THOR utilizing SSS exhibits less severe but non-negligible scatter correction artifacts. For instance, one observes 8.4% activity underestimate (overcorrection of scatter) in the halo region. Moreover, activity distribution around the central sphere is visibly nonuniform and introduces spurious structures in the reconstructed images.
Switching from non-TOF THOR to TOF-THOR while keeping simple SSS is not a viable approach as can be seen in the 3rd column of figure 6: here, scatter is heavily underestimated resulting in 138.3% higher activity around the 'bladder' compared to the ground truth.   Utilizing TOF reconstruction together with TOF-aware scatter correction improves the situation distinctly as can be seen in the last two columns of figure 6. In fact, ISA and TOF-SSS produce near identical results. Both are able to suppress the halo artifact nearly completely and adequately reproduce the correct activity ratios between the different regions. The remaining small differences between ISA and TOF-SSS as shown in the last two columns in the bottom half of figure 6 result in a small overestimate of mean background by 1.5%. All other measures are virtually identical, see table 2.
The difference in scatter time profiles estimated with ISA and TOF-SSS is demonstrated in figure 8. The left plot shows the scatter TOF distribution obtained with TOF-SSS and its ISA approximation for the horizontal LOR shown in figure 5. The right plot shows the profiles-the actually applied TOF weights w i,δ according to (4)-when taking into account the finite time resolution of the PET scanner. It is apparent from these plots that the deviations of ISA from TOF-SSS are mostly very small and larger deviations are sharply localized (notably at the edges of the central peak). The smoothing effect of the finite time resolution eliminates these differences essentially completely as is obvious on the right hand side of figure 8.  The different reconstruction and scatter correction methods were also evaluated in representative clinical PET scans. A whole-body investigation is shown in figure 9 and a long duration/high statistics brain investigation in figure 10 (as mentioned, in the latter case the vendor uses a different reconstruction than for whole-body studies). The results are qualitatively in agreement with those from the phantom study. Due to the less extreme contrast conditions deviations from the TOF-SSS results are visually less obvious, but quantitatively relevant (except for ISA which again yields virtually the same results as TOF-SSS) as can be appreciated from the difference images in the bottom half of both figures.
The runtime performance of the implemented approaches for these clinical data is summarized in tables 3 and 4. As can be seen, ISA outperforms TOF-SSS by a factor of 2-5 and is only 16%-24% slower than non-TOF SSS. The performance gain of ISA over TOF-SSS translates into nearly a factor of two acceleration of the wholebody reconstruction, resulting in a reconstruction time comparable to the one required when using ordinary SSS which indicates that scatter correction is responsible for a substantial part of total reconstruction time. Reconstructing a single bed position of the the same study, we observe acceleration of the reconstruction by a factor 1.4 in comparison to TOF-SSS when using ISA and basically identical performance compared to SSS. For the longduration (large number of events) brain investigation the overall time gain achieved with ISA versus TOF-SSS is only minimal (<4%).

Discussion
The main result of the present investigation is that ISA-based scatter correction yields reconstructed images that are virtually identical to those obtained using full TOF-SSS while accelerating the computation of the required TOF scatter by a factor of 2-5. Our phantom results demonstrate that ISA, like TOF-SSS, is capable of minimizing scatter related image artifacts and deviations from true contrast between different regions even under extreme high contrast conditions. According to these results, ISA could be used as a drop-in replacement for TOF-SSS in clinical PET image reconstruction without compromising quantitative accuracy.
The very similar quantitative performance of ISA and TOF-SSS can be traced back to the fact that the resulting scatter TOF weights w i,δ are almost identical, see figure 8. The TOF-weights are computed from an estimate of the scatter time distributions (corresponding to infinite TOF resolution) by convolution with a scanner-specific TOF kernel. The scatter in the LOR considered for illustration has a non-trivial time profile shape (superposition of a broad and a sharply peaked distribution) caused by immediate vicinity of the central sphere of the phantom. Despite of these challenging conditions, the scatter time distribution estimated by TOF-SSS is remarkably similar to its ISA approximation. Remaining differences are negligible in view of the scanner's limited TOF resolution. The scatter TOF weights in figure 8 were calculated for the modest TOF resolution of the Philips Ingenuity PET/MR scanner used in the present study, but the same reasoning will remain valid for any realistic TOF resolution (say 200 ps). Altogether, it is the good approximation of the scatter time distributions which justifies the ISA approach.
In comparison to ISA and TOF-SSS, the other reconstruction and scatter correction schemes have to be considered unsatisfactory to a varying extent. The discrepancies are most pronounced in the phantom measurements shown in figure 6 due to the chosen very high contrast but are also visible in the clinical data sets in figures 9 and 10. Regarding the vendor provided reconstruction, no implementation details are available and we are not in a position to explain the reason for the observed overcorrection of scatter in most areas. The effect becomes severe at high contrasts and generates notable artifacts not only in the phantom but in the whole-body scan in figure 9, too.
Compared to non-TOF THOR, ISA and TOF-SSS reconstructions clearly show the benefits of fully utilizing available TOF information in order to improve quantitative accuracy and avoidance of image artifacts. Regarding specific reasons for inferior performance of non-TOF THOR, we verified that insufficient convergence of non-TOF versus TOF reconstrcution is not the cause (increasing the iteration count did not lead to relevant changes). Rather, the differences seemingly have to be attributed to principally inferior performance of non-TOF reconstruction and scatter correction. Regarding the clinical data sets, they might also be partly explained by the known sensitivity of non-TOF reconstruction to possible inconsistencies/errors in the MR-derived attenuation maps (which are much better tolerated by TOF reconstruction).
Compared to the non-TOF reconstruction, the attempt to combine TOF reconstruction with non-TOF scatter correction (by distributing the scatter equally over all TOF bins) leads to distinctly more severe image artifacts and increased quantitative errors resulting from a rather pronounced under-correction of the scatter (and thus a positive bias in the resulting images). Explanation of this observation is rather straightforward if one considers the scatter TOF profiles in figure 8 which reveal that the scatter is by no means uniformly distributed over the TOF bins but rather localized within the imaged object boundary. Disregarding this behavior results, in our case, in scatter underestimation by a factor of up to 6.
The speed advantage of ISA over TOF-SSS arises from the fact that it requires neither TOF forward projections in the SSS caching step nor heavy cached arrays manipulations and merging for each scatter point and LOR in the combination step, effectively replacing them by very few (≈5) floating point operations. The reduction in overall image reconstruction time achieved when replacing full TOF-SSS by ISA depends of course on the considered specific implementation of iterative image reconstruction as well as on the input data.
In this study, we have used THOR, our tube-of-response listmode OSEM reconstruction. This type of reconstruction, notably the PSF modeling achieved by the tube-of-response approach, is computationally quite intensive and THOR thus is not as fast as highly optimized scanner specific sinogram based reconstructions. Nevertheless, we observe nearly a factor of two reduction in overall computation time for whole-body studies when performing the initial scatter time distribution computation with ISA rather than with TOF-SSS. This demonstrates that in this case the scatter estimation procedure is responsible for a large fraction of the computation time. The acceleration factor reduces to 1.4 if a single bed position of a whole-body study is being recon- Figure 10. Comparison of scatter correction and reconstruction methods in a clinical brain study. The abbreviations used are clarified in the text. The top row shows a transaxial view of the activity images reconstructed with the different algorithms, and the bottom row the absolute differences between the indicated respective reconstruction method and TOF-SSS. In TOF reconstructions red regions reflect undercorrection of scatter and blue ones overcorrection relative to TOF-SSS. structed whereas the scatter estimation alone is accelerated by a factor of about five, demonstrating the reduced contrib ution of scatter correction to total reconstruction time. For the long-duration/high-statistics brain scan, the overall reconstruction time-which in this case is completely dominated by the actual image reconstruction process-is only reduced by about 4%.
To better understand the rather pronounced differences between the different scenarios (whole-body, low statistics/large object single-bed, high statistics/small object single-bed scan) it first should be noted that the computation times for MLEM update and scatter simulation, respectively, scale differently with scan statistics, number of bed positions, and dimensions of the reconstruction grid. Second, the number of emission/scatter points over which the SC algorithm iterates is much higher for scans of thorax or abdomen than in brain investigations due to distinctly larger transaxial object extension. The number of emission/scatter points is further increased in multi-bed studies by utilizing adjacent bed positions to also account for out-of-FOV scatter which explains the discrepancy between one-bed and multi-bed whole-body reconstruction acceleration factors.
Overall, it is thus clear that the benefits of switching from TOF-SSS to ISA are most important for whole-body studies but also are notable for single-bed scans of the thorax or abdomen. In view of the fact that oncological whole-body studies are currently by far the most frequent application of clinical PET, an achievable factor of up to two reduction in overall reconstruction time for this type of study when using TOF-aware scatter correction seems highly relevant.
The influence of small variations/errors in the underlying attenuation map (e.g. differences between attenuation maps derived in PET/MR and PET/CT, respectively) has been found to be negligible for non-TOF SSS in brain PET (Burgos et al 2014, Teuho et al 2017) and the situation should be similar for whole body PET. ISA does not utilize attenuation information beyond the non-TOF SSS calculation effectivly used for correct scaling of the TOF profiles. It might thus be concluded that ISA is not affected any more by uncertainties of the attenuation map than non-TOF SSS.
ISA also has some potential shortcomings. Unlike TOF-SSS, ISA cannot account for the influence of spatially variant attenuation along ES on the resulting TOF scatter profiles. However, in our data, we were so far not able to see any adverse consequences of this fact, e.g. in the lungs in figure 9. Nevertheless, a more comprehensive investigation of this question might be worthwhile.
A related question concerns a possible adverse effect of ISA on joint attenuation and activity reconstructions such as MLAA and MLACF since these algorithms are known to be sensitive to any kind of data inconsistencies, see Nuyts et al (2018). Another potential shortcoming of ISA is that it does not explicitly account for the influence of the actual scatter medium geometry on the contributions from different emission points. This could be modeled in (6) by an additional emission point position and scatter angle dependent factor accounting for the difference between ISA and TOF-SSS predictions such as those shown in figure 3. However, the presented results and preliminary testing of tentative further modifications to ISA do not indicate that this effect is a practically relevant issue, presumably, because the averaging of the contributions from a huge number of emission points in (6) leads to canceling of the residual differences present for the individual emission points.

Conclusion
ISA is a viable alternative to TOF-SSS accelerating TOF scatter estimation by a factor of up to five. Images reconstructed using ISA are in excellent quantitative agreement with those obtained when using TOF-SSS while overall reconstruction time is reduced by a factor of two in whole-body studies.