Time-based position estimation in monolithic scintillator detectors

Gamma-ray detectors based on bright monolithic scintillation crystals coupled to pixelated photodetectors are currently being considered for several applications in the medical imaging field. In a typical monolithic detector, both the light intensity and the time of arrival of the earliest scintillation photons can be recorded by each of the photosensor pixels every time a gamma interaction occurs. Generally, the time stamps are used to determine the gamma interaction time while the light intensities are used to estimate the 3D position of the interaction point. In this work we show that the spatio-temporal distribution of the time stamps also carries information on the location of the gamma interaction point and thus the time stamps can be used as explanatory variables for position estimation. We present a model for the spatial resolution obtainable when the interaction position is estimated using exclusively the time stamp of the first photon detected on each of the photosensor pixels. The model is shown to be in agreement with experimental measurements on a 16 mm  ×  16 mm  ×  10 mm LSO : Ce,0.2%Ca crystal coupled to a digital photon counter (DPC) array where a spatial resolution of 3 mm (root mean squared error) is obtained. Finally we discuss the effects of the main parameters such as scintillator rise and decay time, light output and photosensor single photon time resolution and pixel size.


Introduction
Gamma-ray detectors based on bright monolithic scintillation crystals coupled to pixelated photosensors are currently being considered for several applications in the medical imaging field such as clinical time-of-flight positron emission tomography (TOF-PET) (Moliner et al 2012), small animal PET , Miyaoka et al 2010, Sánchez et al 2012 and Compton cameras for dose monitoring during hadron therapy treatments (Llosá et al 2013). Typically, the position of interaction of the gamma ray inside the scintillator is estimated by the light distribution (measured on the photosensor pixels) using statistical learning (Bruyndonckx et al 2004, van Dam et al 2011a or model-based methods (Li et al 2010). Such methods have shown to provide excellent spatial resolution in the lateral dimensions  while many promising algorithms are available for the estimation of the depth-ofinteraction (DOI) (Lerche et al 2005, Ling et al 2007, van Dam et al 2011b. One advantage of monolithic detectors (coupled to fast, low-noise photosensors) in applications involving timing (e.g. TOF-PET) is the fact that each scintillation pulse is sampled by many photosensor pixels. For this reason it is possible to measure multiple time stamps corresponding to the first scintillation photon(s) detected by each of the photosensor pixels. It was demonstrated (van Dam et al 2013) that an improvement in the timing performance can be obtained when the information provided by all these time stamps is exploited. Moreover, it was shown that the degradation in the timing performance due to the optical transport of scintillation photons inside the crystal can be partially corrected for by additionally utilizing the information of the position of the gamma interaction provided by the light intensities. In other words, in a monolithic detector the timing information is correlated with the spatial information via the optical transport of the scintillation photons. In this work, we further study this correlation and investigate the possibility of using the spatio-temporal distribution of the time stamps registered in a scintillation event to predict the position of the interaction point. From this perspective, statistical fluctuations in the registration times of individual scintillation photons translate into an uncertainty in the spatial localization of the interaction point. The ability that time stamps have to estimate the interaction position of a gamma ray, defined as the spatial resolution, will therefore depend on temporal parameters (e.g. rise and decay time of the scintillator and photosensor single photon time resolution) as well as on the photosensor pixel size and the scintillation light yield. In the following sections, a model predicting the spatial resolution as a function of each of these parameters is presented and compared to experimental data acquired on a 16 mm × 16 mm × 10 mm Ca-codoped LSO : Ce crystal coupled to a DPC array. Finally, the effect of each of the parameters on the spatial resolution is discussed.

Theory
In this section we derive a model for the spatial resolution that is obtained when the coordinates of the interaction point are estimated using the time stamps only (i.e. we ignore the light intensities). We start by assuming that P photosensor pixels detect the optical photons emitted by a scintillation event in which a gamma ray deposits its full energy γ E at position x = (x,y,z) T (or within a volume around x whose dimensions are small compared to the spatial resolution of the system) and time = t 0 0 inside a scintillation crystal with light yield Y. Each photosensor pixel i detects a fraction f i of the optical photons so that the number of detected photons n i on each pixel will follow a Poisson distribution (assuming f i is small) with expected value [ ] = γ E n E Yf i i (Barret and Myers 2004). Let t i delay be the delay between the gamma interaction time and the time a scintillation photon is actually registered in photosensor pixel i. In general, the t i delay of the scintillation photons can be considered to be statistically independent and identically distributed (i.i.d.) according to the probability density function (pdf) ( ) p t t i delay . The three main contributions to t i delay are the time of emission t e of the scintillation photon relative to the gamma interaction time, its optical propagation time t i op and the signal transport delay t s : In the following we discuss the three contribution separately. The pdf of the scintillation photon emission times t e is described as the sum of two exponential functions of the form: In this work we assume a single scintillation process, with rise time constant t r and decay time constant t d . However, equation (2) can be easily generalized to the case of multiple scintillation processes (ter Weele et al 2014).
The optical propagation delay t i op > 0 is the time between the actual emission time at point x and the time at which the photon is absorbed in pixel i. In order to emphasize the dependence . In practice t i op represents the transport of scintillation photons inside the crystal. Its distribution depends on several factors such as the dimensions of the scintillator, its refractive index, the reflective wrapping and the roughness of the crystal surface. The dependence on x can be intuitively explained by geometrical considerations, that is, smaller t op are expected when the interaction point x is closer to pixel i.
The signal transport time t s represents the total electronic delay and consists of components that find their origin within the photosensor as well as in the readout electronics. In general, the signal transport time t s and its pdf ( ) p t ts will depend on the type of photosensor and readout electronics used. For the digital photon counter considered in this work, which is described in the next section, ( ) p t ts can be modeled as a normal distribution with mean μ s and standard deviation σ μ ≪ s s , hereinafter referred to as the single photon time resolution (SPTR) (in the literature sometimes also referred to as photosensor transit time spread).
Since the three contributions to t i delay expressed by equation (1) are statistically independent, the pdf of t i delay is given by: where * denotes the convolution operator. In the following we restrict ourselves to an example setting in which a photodetector pixel can only register the time of arrival of the first scintillation photon (while the position of the absorption point of the scintillation photon within the pixel is unknown). That is to say, among all n i photons detected by pixel i, only the time stamp t i first of the photon with the smallest t i delay is registered. The pdf of t i first is given by the first order statistics (or sample minimum) (David 1989): is the cumulative distribution function (cdf) of t i delay . For simplicity, in equation (4) we replaced n i by its expected value [ ] E n i as is also done in (Seifert et al 2012), which is justified if [ ] ≫ E n 1 i . Since the time stamps of the first photons registered by each photosensor are statistically independent it follows that the joint distribution (or likelihood function) of , , P first first 1 first 2 first is given by: The likelihood function can then be used to find the posterior probability ( ) x t p | x first using Bayes theorem: where ( ) x p x is the prior, representing the probability of a gamma interaction at a given point in the crystal, while ( ) t p tfirst can be considered as a normalization factor since it does not depend on x. The prior probability depends in general on the irradiation setup. An accurate model for ( ) x p x can be obtained for instance from Monte Carlo simulations of the transport of gamma rays within the scintillation crystal. Here we approximate ( ) x p x by a uniform distribution, i.e. it is constant for x inside the crystal volume, and zero outside.
We recall that the posterior probability ( ) x t p | x first represents the probability that the gamma interaction occurred at point x, having measured the set of timestamps t first . The definition of a single point estimator x will depend on our choice for the spatial resolution. In this work we define the spatial resolution for each of the three coordinates σ σ σ σ = ( ) , , x x y z as the root mean squared error (RMSE): where E[] denotes expectation with respect to x and t first . The mean square estimate which minimizes (7) is then defined as (Lehmann and Casella 1998): It is worth noting that if the timestamps of the first photons had no correlation with the position of the gamma interaction, i.e.
, the posterior of equation (6) would be equal to the prior (i.e. we gain no information on x by measuring the time stamps). In this worst case scenario, the spatial resolution for a flat prior would be equal to the RMSE of the quantization error: , , x y z MAX (9) where L L L , , x y z are the crystal dimensions. In other words, σ MAX represents the spatial resolution we obtain if we always estimate the interaction position as the center of the crystal, regardless of the time stamps we measure.
It must be noted that the model assumes that the t first are defined with respect to the gamma interaction time = t 0 0 . In general, the gamma interaction time is not known and the time stamps are measured with respect to a reference clock that can be arbitrarily chosen. Thus, assuming  (10) It is emphasized that t start is not correlated with the time of interaction of the gamma photon and thus carries no information on the interaction point. To get rid of the dependency on t start we subtract the average value Geometrically this is equivalent to projecting the vector t obs along the direction ( … ) 1, 1, , 1 . The transformed vector: where P ( … ) 1, ,1 denotes the projection matrix, lies on the hyperplane ∑ = t 0 i obs and has − P 1 degrees of freedom. The likelihood of the linearly transformed time stamps is given by:

Materials and methods
The spatial resolution of a 16 mm × 16 mm × 10 mm detector obtained using time-based position estimation was obtained both experimentally and from the model of equation (7). We start by describing the experiment, which is also the starting point for the model calculation, as detailed in the following. The array consists of 16 autonomous sensors (DPC dies) placed at a pitch of 8 mm into a 4 × 4 array. A die measures 7.875 mm × 7.15 mm and is subdivided into 2 x 2 so-called DPC pixels, whose dimensions are 3.875 mm × 3.2 mm. The four DPC pixels share a single time-to-digital converter (TDC). Each DPC pixel comprises a total of 3200 microcells. One microcell is a single photon avalanche photodiode (SPAD) connected to a dedicated CMOS quenching and digitization circuit. The DPC sensor operates on the basis of an asynchronous, configurable-length acquisition cycle. The acquisition sequence of a die is started by a trigger, whose threshold can be set by the user. In this work the trigger level MT =1 was used, i.e. a trigger is generated as soon as a microcell discharges on the die. Whenever a trigger is generated, a time stamp is registered and the die starts the acquisition cycle. At the end of an acquisition cycle a die outputs one time stamp and the number of fired cells on each of the DPC pixels. It must be noted that the DPC sensor is unable to identify which microcell within the die generated the trigger. Because the DPC array is larger than the crystal used in this work, only the four central dies of the DPC were enabled and optically coupled to the crystal using a transparent silicone material (Sylgard 527, Dow Corning). The sides of the crystal were wrapped with a specular reflector foil (Vikuiti ESR, 3M) whereas the top surface was covered with Teflon tape. The detector was operated at a temperature of −25 °C.

Irradiation setup and data acquisition.
The monolithic detector under investigation was irradiated with a perpendicularly incident pencil beam of 511 keV annihilation photons from a 22 Na point-source (Ø 0.5 mm). A detailed description of the collimator setup used for this measurement is reported in (Borghi et al 2015). The x-y surface of the crystal (measuring 16 mm by 16 mm) was irradiated at a grid of 64 × 64 reference positions, with a pitch of 0.25 mm (figure 1). At each irradiation point, 50 events in the energy photopeak were registered. Each event of the data set consists of a set of four time stamps (one time stamp per DPC die) corresponding to the time stamps of the first detected photon on the die.

Data processing and analysis
First the data were linearly transformed as described in equation (11). Then the data set was split into two subsets. One was used for testing (test data set) and the other one for training (training data set).
The spatial resolution was computed on the test dataset separately for x and y as follows: Here, N test is the number of events in the test data set, test are the coordinates of the entry points of the gamma photons on the crystal front surface and x i test are the corresponding estimated coordinates. The estimation was done non-parametrically using the 'kknn' package of the R environment (R-Core-Team 2015). This package implements k-nearest neighbors regression (Hechenbichler and Schliep 2004). The method works as follows: first, for each event i with timestamps t i test in the test dataset, the Euclidean distances of t i test to all of the events in the training set are computed and ordered. Then the subset N i of the training set containing the k events with the smallest distances (nearest neighbors) is selected. The position of interaction x i test subsequently is estimated as the average interaction position of the k-nearest neighbors: where x ji train are the positions of irradiation corresponding to the k-nearest neighbors of t i test in the subset N i . In practice k-nn represents a direct implementation of equation (8) with the two approximations that (i) expectation is approximated by averaging over sample data, (ii) conditioning on a point is replaced by conditioning on a neighborhood of that point.

Model calculation
In order to compare the model predictions to the experimental results, we chose the model parameters to represent our experimental set-up (described previously) as accurately as possible. Calculation of equation (3) requires a model of the optical transport inside the crystal for which no analytical expression is readily available. Therefore a Monte Carlo simulation of a 3D grid (with a pitch of 1.6 mm) of optical-photon point sources inside the crystal was performed using GATE (Jan et al 2004). The crystal (sized 16 mm × 16 mm × 10 mm) was wrapped with diffuse reflector material with reflectivity 0.99. Although this choice of reflector does not exactly match our experimental setup, we believe such discrepancy has negligible effects as in our experience the behavior at the boundary (determined by the surface finish and refractive indexes) is the most relevant. The geometry of the photosensor array resembled that of the DPC array used in the experiment, as shown in figure 2. It must be noted that one photosensor pixel as defined in the model of section 2 corresponds to one DPC die (and not to a DPC pixel) since the DPC used in this work has one TDC per die. A photon detection efficiency (PDE) of 0.3 was assumed for the active area of each photosensor pixel, while the gaps between them were modeled as black absorbers. The main contribution to the single photon time resolution in a DPC array is given by the skew of the trigger network (Frach et al 2009). A value of σ = 71 ps s (168 ps FWHM) is assumed for our model (Brunner 2014). A summary of the main parameters used in the model can be found in table 1.
For every optical point source 0.5 million optical photons (generated at time = t 0 0 ) were tracked and the time of arrival was recorded for each of the photons detected by pixel i. The fraction of detected photons per pixel f i was calculated as the ratio between the number of detected photons on pixel i and the total number of photons generated. The scintillator light yield Y then was estimated in such a way that the expected total number of photons detected for a 511 keV gamma would match the experimental mean number of photons measured in the photopeak (≈ 3000 photons).
The probability density of the optical propagation time ( ) x p t t i op was computed non-parametrically with kernel density estimation using a rectangular window and a bandwidth of 2 ps. Equations (3) and (4) were solved numerically while the integral of equation (7) was estimated by Monte Carlo integration. Additionally, we considered the case where the gamma interaction time is unknown and the time stamps are linearly transformed according to (11). Finally, the effects of each of the main physical parameters on the spatial resolutions were studied. In particular the influences of the decay time, the rise time, the scintillator light output and the single photon time resolution were investigated by sweeping the value of the parameter under investigation while keeping all the others constant at the values reported in table 1. To study the effect of the photosensor pixel size the model was recomputed using a hypothetical pixel size of 3.2 mm × 3.9 mm (leading to 16 time stamps measured for each gamma interaction).

Results and discussion
In figure 3, the likelihood function calculated using equation (4) with parameters reported in table 1 is plotted for two different positions of irradiation about 9 mm apart. Noticeable changes in the likelihood function can be observed, which indicate a correlation between the position of interaction and the registered time stamps. Computation of the spatial resolution for this setting, averaged over the entire crystal, resulted in 3.07 mm, 3.12 mm and 2.15 mm for σ x , σ y and σ z respectively. It can be noted that the time stamps indeed carry information on the position of interaction. In fact, ignoring such information would result in a spatial resolution σ MAX of 4.62 mm, for x and y, and 2.89 for z according to equation (9), which represents the RMSE of the quantization error. Table 2. Experimental spatial resolutions averaged over the entire crystal measured for the two transversal coordinates. Position estimation was done using non-parametric k-nn regression as described in section 3. The value of k was found using 5-fold cross validation.

Coordinate
σ exp k N train x 3.04 mm 164 50 000 y 3.12 mm 134 We emphasize that the above results were calculated assuming that the gamma interaction time is known and equal to zero. To study the more realistic case of an unknown gamma interaction time, the spatial resolution was also calculated for the linearly transformed time stamps ′ t first using the likelihood function of equation (12). It was found that using the transformed ′ t first instead of the t first results in a negligible degradation (<1%) of the spatial resolution. In other words, the information on the position of the gamma interaction is mainly in the relative differences of the time stamps measured across the photo sensor array and not in their absolute values.
The experimental results for the spatial resolution, averaged over the entire crystal, are reported in table 2 together with the other parameters used for the position estimation. Spatial resolutions of 3.04 mm and 3.12 mm were found for the x and y coordinates respectively, which is in good agreement with the model. The difference in performance for the two coordinates can be explained by the asymmetries of the photosensor die in the 2D (see section 3.1.1). In order to study the degradation in the spatial resolution due to the crystal edges, typical of monolithic detectors (Seifert et al 2013, Vinke andLevin 2014), we plot in figure 4 the experimental and the modeled spatial resolution in the x coordinate as a function of x. A significant worsening of the spatial resolution at the border of the crystal can be noticed, which is attributed to surface reflections and the relatively large size of the photosensor pixel (DPC die) compared to the crystal dimensions.
The results of the analysis of the influence of the model parameters are shown in figure 5 for the x coordinate where the blue (solid) and red (dashed) lines indicate model calculations assuming a photosensor pixel size of 7.15 mm × 7.9 mm and a hypothetical size of 3.2 mm × 3.9 mm, respectively. It can be seen that the main limiting factor of our setup is the photosensor pixel size (corresponding to a DPC die), while the scintillation rise time (figure 5(a)) and photosensor SPTR (figure 5(b)) have a negligible effect over the range of values considered. On the other hand, light output (more generally the number of detected photons, figure 5(c)) and decay time (figure 5(d)) have a stronger influence on the spatial resolution. We emphasize that the dependence on timing parameters such as the scintillator decay time is particularly interesting for instance in applications like time-of-flight PET, in which the use of a faster scintillator would simultaneously improve spatial resolution and coincidence resolving time. For instance, the model indicates that with a scintillator as bright and fast as LaBr 3 : Ce (t d < 16 ns, Y > 70 000 ph MeV −1 ) a spatial resolution of less than 2 mm in principle is possible from the time stamps only. It must be noticed, however, that other performance parameters such as energy resolution and sensitivity would be affected by the choice of the crystal material and therefore should be also taken into account when designing a gamma ray detector.
In figure 6, we plot the spatial resolution calculated with the model of section 2, for hypothetical fast detectors (τ d < 15 ns) with relatively low light output. In this case, traditional positioning methods based exclusively on the light intensities are expected to perform poorly since their RMSE depends mainly on the number of detected photons . On the other hand, it appears that with time-based position estimation the fast detector response can compensate for the low light signals so that spatial resolutions of less than 2 mm are achievable with less than 1000 detected photons provided that the scintillator decay time τ d < 5 ns.
In general, both the light intensities and the time stamps will be available. Given the encouraging results obtained here, in principle one may think of combining them to increase the accuracy of the positioning algorithms. For instance, this can be done by (i) fusing the two types of variables into a single composite vector and then train a statistical algorithm on this composite vector as its input; (ii) combining the two types of information on a higher level, that is, the vector of time stamps and the vector of light intensities are used as separate inputs to some statistical algorithm. The two outputs are then combined through some integration technique. An investigation on the best strategy and the potential improvement obtainable from such an approach are however outside the scope of the present paper and are left for future work.

Conclusions
In this work we showed both theoretically and experimentally that in state-of-the-art monolithic scintillator detectors the first-photon time stamps recorded by each photosensor pixel are correlated with the position of interaction of the gamma ray. Such correlation suggests that the time stamps can be used as explanatory variables to estimate the interaction position. We presented and discussed a model that relates the spatial resolution obtainable from the timebased position estimation to the main parameters such as scintillator rise and decay time, light output, single photon time resolution and photosensor pixel size. Further work is needed to investigate the potential improvement in the spatial resolution that could be derived from the combined use of the time stamps and the light intensities when both are available.