Time-of-flight (TOF) implementation for PET reconstruction in practice

The time-of-flight (TOF) feature of PET scanners has been used for a long time in PET reconstruction, but many implementational aspects are still incomplete or ambiguous in the literature. Here we formalize and present theoretical and practical implementation details for the reconstruction of clinical TOF histogram and list-mode data using ML-EM. Relevant aspects include the computation of the TOF component of the system matrix, the processing of TOF bins, the use of estimations of random and scattered coincidences, and differences between histogram and list-mode ML-EM TOF reconstruction. Several approaches and approximations have been implemented in the CASToR platform and compared for OSEM reconstructions of patient data from the GE Signa PET/MR scanner. Differences between implementations are not larger than the typical bias in clinical data reconstruction. The largest difference and contrast loss occur when the processing of histogram TOF bins is simplified, and list-mode reconstruction is most sensitive to the truncation of the Gaussian TOF probability distribution.


Introduction
There is a large amount of literature on different aspects of using time-of-flight (TOF) measurements in PET scanners. However, implementation details for reconstruction algorithms are often incomplete and scattered in various papers dealing either with reconstruction or with related topics (i.e. TOF PET scanners, estimations of random and scattered coincidences). Also, there is no formalization nor comparison of different implementation approaches and approximations. Here we focus on standard quantitative reconstructions from clinical data and on usual primary data formats (list-mode and histogram or raw sinogram), but the considerations are more general and apply to a wide family of algorithms and data formats. The optimization of LOR and TOF binning is not relevant for our purposes and we do not consider the family of TOF reconstruction methods based on preprocessed data formats (i.e. methods performing a single TOF backprojection and iterative operations in the image domain, Snyder and Politte (1983), Matej et al (2009)). First, we formalize the implementation of TOF and investigate some details and approximations. For each aspect, we mention articles that provided some theoretical or practical details. Then, we compare several implementations to investigate the impacts of approximations, using the CASToR platform, Merlin et al (2018), CAS (2017). This paper is also in support of the implementation of TOF in the CASToR platform.

Theory
The TOF measurement for a pair of detected coincident photons is defined as the difference of photon arrival times ∆t, but in practice this time difference is converted into v, the presumed spatial position of the annihilation on the LOR. Here v is defined as the position with respect to the center of the LOR, and thus is expressed as v = c∆t 2 , c being the speed of light.

Time-of-flight (TOF) implementation for PET reconstruction in practice
The uncertainty of TOF measurements is usually modelled with a normalized (integral =1) Gaussian function centered at v (Gaussian probability distribution). The Gaussian model and its parameters may be open to consideration because of imperfections of the imaging system, as for instance the spatial inhomogeneity of the timing resolution, the issue of timing alignment, Clementel et al (2013), the dependence on the count rate, Levin et al (2016), or the accuracy of the estimation of the standard deviation σ, Daube-Witherspoon et al (2006).

Data formats
TOF measurements represent essentially additional information and do not modify other aspects in the acquired data. TOF histogram data contain an additional data matrix dimension for TOF bins: all the counts acquired in a LOR with a certain range of associated TOF measurements are summed into the corresponding TOF bin for that LOR. This implies the sum consistency property: the counts y ib for all TOF bins b in the LOR i must sum to the total counts y i detected in the LOR, b y ib = y i . (1) TOF list-mode data contain an additional TOF measurement value for each detected count, where the measurement can be continuous, though in practice it is necessarily saved with a limited precision, using a certain quantization step.

Maximum-likelihood expectation-maximization (ML-EM) equations
Let A be the system matrix, λ the radioactive concentration, r the expected count rate of random coincidences, s the expected count rate of scattered coincidences and j the voxel index. The multiplicative update part u of the ML-EM equation λ t+1 j = u t j λ t j for voxel j and iteration t differs for each data format. For TOF histogram data, the measured data are modelled with a Poisson distribution with expectation and thus the update term for histogram data is This equation is used in the papers providing details about TOF histogram reconstruction and in the CASToR platform. The equation for list-mode data with quantized TOF measurements can be derived from the histogram equation (3) as One could also write an equation for list-mode data with continuous TOF measurements, where the system matrix elements and background contributions are expressed as continuous functions of the TOF measurement v n, This equation is used in the CASToR version 2. The equation for list-mode data with quantized TOF measurements can also be derived directly from this list-mode equation with continuous TOF measurements (5) as without going through the histogram equation (3) (the superscript q stands for 'quantized').

System matrix elements and TOF weights
Ideally, the system matrix should be modelled by taking into account all the aspects of the acquisition process at once, but in practice it is convenient and reasonably accurate to approximate the system matrix elements as a multiplication of independent terms. Let A ij be the complete explicit system matrix elements for non TOF data, containing the geometric projection and other components. Let w be the TOF weights, independent of A ij , such that the system matrix elements for TOF histogram data are A ijb = A ij w ijb , for TOF list-mode data . From the sum consistency (1) and from the model (2) follows an important sum consistency property for TOF weights w: the weights for a single LOR, a single voxel, and all TOF bins b or all TOF measurement values v must sum or integrate to 1, and thus represent a probability distribution.
There is a practical limitation to this property, because of the finite range of TOF measurements. This range can be called the TOF field-of-view and usually matches the size of the coincidence timing window. For voxels close to extreme TOF bins or measurements, the sum/integral (7) decreases. An illustration of this effect is shown in figure 1. The importance of this issue depends on the size of the imaged subject with respect to the scanner. For instance, we may say that the issue becomes of little importance if the distance between the edges of the patient body and of the TOF field-of-view is larger than 3 standard deviations of the TOF uncertainty function (75 mm for the Signa PET/MR). It should be noted that theoretically the sum can never reach exactly 1, because the TOF Gaussian uncertainty function is infinite and the range of TOF measurements finite, but in practice, for instance in the test in figure 1, the sum is equal to 1 in the scanner center.

Voxel sensitivity
As a consequence of the sum consistency property (7), the voxel sensitivity factors are always the same, regardless of the use of TOF and regardless of the data format (the only paper that mentions this property is Groiselle and Glick (2004)),

TOF bin definition
The notion of TOF bin for TOF histogram data is frequently used and straightforward. However, it does not account for the approach (6) for list-mode data with quantized TOF measurements, for the following reasons: this approach is not derived from the TOF histogram equation (3) and it deals with quantization instead of histogramming. In order to provide expressions for TOF weights and estimations of random and scattered coincidences for all types of ML-EM equations using the same formalism, we redefine the term 'TOF bin', so that it refers either to (a) the quantization bin, which only makes individual TOF measurement values less precise, without implying any summing of counts in the raw data, or to (b) the cumulative bin, which implies the accumulation or summing of counts with a certain range of associated TOF measurement values. Hence, the quantization bin is used in (6) and the cumulative bin is used in (3) and (4). The distinction is relevant because the TOF weights and the estimations of expected count rates of random and scattered coincidences have to be computed accordingly. Both types of TOF bin can be represented with a box function Π b , with width ∆v b , center v b , and value h, as (9). The cumulative TOF bin has h = 1 and its integral equal to the bin width ∆v b , whereas the quantization TOF bin has h = 1/∆v b and its integral equal to 1.

Computing TOF weights
Here we investigate different approaches to computing the TOF weights, under the assumption of a Gaussian TOF uncertainty model. TOF weights can be expressed from two equivalent points of view: (a) by considering the contribution of each voxel to a TOF bin or measurement, or (b) by considering an uncertainty function associated to a TOF bin or measurement, and sampled for each voxel. Usually, TOF bins are equidistant and have the same width, which simplifies computations greatly. When computing the TOF weight for a list-mode detected coincidence n with its associated continuous TOF measurement v n , it is equivalent to (a) center the Gaussian function at the voxel projection onto the LOR v j and sample it at v n and (b) center the Gaussian function in v n and sample it at v j (see figure 2 left).
TOF weights for TOF bins can be formulated in the same manner for both cumulative (w ijb ) and quantization (w q ij (v b )) bins, based on several approaches: The contribution of a voxel j to a bin b for a LOR can be expressed using (10a) (see figure 2 center). This was explicited in Mehranian et al (2016) for cumulative TOF bins. Equivalently, it can be expressed by first computing the convolution (10b), which is a function of the distance between the TOF bin center and the voxel projection onto the LOR, and then sampling it for the given voxel and TOF bin (see figure 2 right). These two approaches are equivalent in theory, but in practice (10a) is more suited for computation on the fly and (10b) is more suited for precomputation (slight differences may occur depending on the precision of the precomputed function, see figure 3). It should be noted that the difference between weights for cumulative bins w ijb and for quantization bins w q ij (v b ) boils down to a single multiplicative factor equal to the TOF bin width: w q ij (v b ) = w ijb /∆v b . Therefore, the two approaches to list-mode reconstruction with quantized TOF measurements (4) and (6) are equivalent, the difference being conceptual.
Various approximations can be used for simplifying computations. TOF quantization bins can be simply neglected, which implies reverting to the continuous list-mode (5). The weights for cumulative TOF bins can be approximated by simple multiplication of a Gaussian sample with the width of the TOF bin (10c). This is equivalent to the approach in Lois et al (2010), Watson (2007), and in the CASToR version 2, where there is no concern about absolute value: the TOF weights are computed for all TOF bins for a voxel and a LOR, and then normalized by their sum to satisfy the sum consistency (7). Equation (10c) is a crude approximation to (10a) or (10b), and can be viewed as a decrease of the nominal σ. There may be other approximations in between, such as (10c) with a larger σ that accounts for the TOF bin width.
The nominal Gaussian function can be truncated, for instance to a certain number of standard deviations. Truncation implies setting some TOF weights and consequently system matrix elements to zero. This impairs the sum consistency (7), which may have the following consequences: (a) if the voxel sensitivity factors are computed using the system matrix with truncation, they will differ from non TOF sensitivity factors and (8) will not be satisfied anymore, (b) if the voxel sensitivity factors are computed using the system matrix without truncation, (8) will be satisfied but discrepancies will occur between the system matrix elements used Figure 2. Illustration of the computation of TOF weights for continuous TOF measurements (left), for TOF bins using approach (10a) (center) and for TOF bins using approach (10b) (right).
for projection/backprojection operations and the voxel sensitivity factors. To ensure that (8) is satisfied, the truncated Gaussian function may be normalized so that its integral equals 1. However, such a function would represent a different TOF uncertainty model with a different standard deviation.
Building the TOF weights independently from the projection component p of the system matrix implies that the TOF component has a constant value over the volume of a single voxel. In theory, it may be more accurate to account for TOF directly in the formulas for computing the projection component of the system matrix p ijb or p ij (v). In this case, a more general sum consistency property must be satisfied: b p ijb = p ij and p ij (v)dv = p ij . However, in practice, this question has not been studied and may become relevant if the TOF uncertainty is reduced to the order of the voxel size.

Random coincidences
The random coincidences depend on the detection process and their associated TOF measurements may be considered random as well. In view of the sum consistency (1), we want a model (2) such that br ib =r i and r i (v)dv =r i . Let us consider a widespread method for estimating the expected count rate of random coincidences for a LOR i, r i = S i1 S i2 2τ , where S i1 and S i2 are the rates of single counts detected at LOR end points, and τ is the timing window for coincidence detection. The expected random count rate per cumulative TOF bin r ib can be estimated by replacing 2τ with the temporal size of the TOF bin ∆t b . This is equivalent to dividing the standard estimation r i by the number of TOF bins. This estimation was mentioned in Haynor et al (1988), and used for list-mode reconstruction with quantized TOF measurements (4) in Wang et al (2006), Zhang et al (2017). The continuous expected random count rate r i (v) can be regarded as a limit case when the width of the TOF bin approaches 0. As it is a function of the distance v along the LOR, it can be estimated by dividing the estimation r i with the spatial total span of TOF measurements (∼spatial equivalent of the coincidence timing window), which gives a constant 2S i1 S i2 /c. The estimation r q i (v b ) for quantized TOF measurements in (6) also equals the same constant. Hence, for a LOR, the estimation of the random count rate is the same for all TOF bins, or equivalently for all continuous or quantized TOF measurements.

Scattered coincidences
The estimation of scattered coincidences must also include the estimation of their associated TOF measurements. Two main methods have been developed for the estimation of the expected count rate of scattered coincidences: Watson (2007) provides scatter estimation s ib for cumulative TOF bins and Werner et al (2006) provides a continuous scatter estimation s i (v). As for random coincidences, we want a model such that bs ib =s i and s i (v)dv =s i . Depending on the availability of estimation methods, the cumulative and quantization bin estimations s ib and s q i (v b ) can be obtained from continuous estimations as and continuous estimations can be approximated from cumulative bin estimations as Figure 3. b w ijb for all voxels j along the horizontal LOR passing through the center of the scanner, for different implementation approaches with truncation of the Gaussian function at 3σ.

Method
All the approaches mentioned in this paper have been implemented in the CASToR platform and will be included in the CASToR version 3. Hence, TOF weights for TOF bins can be computed using (10a) (integral computation on the fly using erf), (10b) (precomputed convolution sampled on the fly) and (10c) (either precomputed function sampled on the fly or computation on the fly using exp). TOF weights for continuous TOF measurements can be computed on the fly or precomputed. The reconstruction of list-mode data with quantized TOF measurements is implemented using (6) only, which is equivalent to (4). The Gaussian distribution can be truncated to a given number of standard deviations σ. Different implementations are referred to as equation (10) t number of σ for truncation if present, for instance Ha_t3 for histogram reconstruction with (10a) and truncation to 3σ.
Whole body scans of two patients with high body mass index (BMI) were selected from the GE Signa PET/ MR scanner (GE Healthcare, Milwaukee, WI, USA). The scanner TOF characteristics were: 13.02 ps (1.95 mm) quantization TOF bin, 169.26 ps (25.37 mm) cumulative TOF bin, and 420 ps (62.96 mm) nominal TOF resolution (FWHM). The data were formatted into a TOF histogram and a list-mode with quantized TOF measurements. The reconstruction parameters were similar to parameters in clinical routine: OSEM reconstruction with two iterations and 27 subsets, voxel size 2.34 × 2.34 × 2.78 mm 3 , sieve reconstruction with kernel equal to the scanner resolution PSF (Snyder et al 1987, Stute andComtat 2013).
It should be noted that there are differences in computation of voxel sensitivities for different data formats. For list-mode data, sensitivity is precomputed using all the LORs and using system matrix elements without TOF weights, because of (8). Hence, inconsistencies may occur between the sensitivity and the system matrix elements when approximations to TOF weights are used. For histogram data, there are no such inconcistencies because the sensitivity is computed on the fly for each OSEM subset using approximated TOF weights.
We tested all the approaches to the computation of TOF bin weights (10a)-(10c), without and with Gaussian truncation at 3σ. All the approaches for a data format were compared to the reference approach for that data format, which was (10a) without truncation. Differences were evaluated using the root-mean-square error presented as a % of the mean value in the image, using difference images and using SUV quantification in tumoral regions of interest (SUV mean, max, peak). Figure 3 shows an illustration of the effect of Gaussian truncation on the sum consistency (7). The implementation Ha_t3 presents a sum lower than 1 (∼0.9973), stable along the LOR. The implementation Hb_t3 presents slight fluctuations due to the precision of the precomputed convolution (here set to 0.01 mm to make the effect visible), and Hc_t3 presents larger fluctuations due to the approximations.

Results
Here we expose only some noteworthy differences resulting from the comparison of images reconstructed using different implementations. The order of magnitude of differences is presented in table 1 and relevant images of differences are presented in figure 4.
For histogram data, the approximation of (cumulative) TOF bin weights (Hc versus Ha) resulted in the largest differences and in contrast loss (see figure 4(b)). Truncating the Gaussian (Ha_t3 versus Ha) produced lower differences, partially due to the fact that the sensitivity is consistent with the computed truncated TOF weights. For list-mode data, the approximation of (quantization) TOF bins (Lc versus La) produced the same effect as for histogram data, but with much lower differences, because the TOF bin is one order of magnitude smaller. Truncating the Gaussian (La_t3 versus La) resulted in comparatively larger differences and in a general loss of intensity (see figure 4(d)), due to the inconsistencies between the voxel sensitivity and the computed truncated TOF weights.
Direct comparison between the histogram and the list-mode reconstruction is not shown, because the difference depends on several factors, not necessarily related to TOF. In order to perform a fair comparison, the number of subsets should be 1, the voxel sensitivity factors strictly identical, the TOF bin identical, and (4) used for list-mode reconstruction. The main reason for using approximations is reducing the implementation complexity and the computation time. These depend strongly on the software framework and on computer architecture. In CASToR implementation, the approximations do accelerate computations up to a factor of 4, but a detailed analysis of computation time was not an aim of this paper. For the scanner and setup used in this study, the order of magnitude of differences between implementations (less than a couple of percents) is not larger than the order of magnitude of typical bias in standard clinical OSEM reconstructions (Schöder et al 2004). The differences are in accordance with the presented theoretical and practical implementation details.

Conclusion
We provide theoretical and practical details for implementing TOF in PET reconstruction. Several approaches have been implemented in the CASToR platform and shall be available in the next release (version 3). The differences between approaches are shown to be of low importance for standard reconstructions of clinical data.