Towards monolithic scintillator based TOF-PET systems: practical methods for detector calibration and operation

Gamma-ray detectors based on thick monolithic scintillator crystals can achieve spatial resolutions  <2 mm full-width-at-half-maximum (FWHM) and coincidence resolving times (CRTs) better than 200 ps FWHM. Moreover, they provide high sensitivity and depth-of-interaction (DOI) information. While these are excellent characteristics for clinical time-of-flight (TOF) positron emission tomography (PET), the application of monolithic scintillators has so far been hampered by the lengthy and complex procedures needed for position- and time-of-interaction estimation. Here, the algorithms previously developed in our group are revised to make the calibration and operation of a large number of monolithic scintillator detectors in a TOF-PET system practical. In particular, the k-nearest neighbor (k-NN) classification method for x,y-position estimation is accelerated with an algorithm that quickly preselects only the most useful reference events, reducing the computation time for position estimation by a factor of ~200 compared to the previously published k-NN 1D method. Also, the procedures for estimating the DOI and time of interaction are revised to enable full detector calibration by means of fan-beam or flood irradiations only. Moreover, a new technique is presented to allow the use of events in which some of the photosensor pixel values and/or timestamps are missing (e.g. due to dead time), so as to further increase system sensitivity. The accelerated methods were tested on a monolithic scintillator detector specifically developed for clinical PET applications, consisting of a 32 mm  ×  32 mm  ×  22 mm LYSO : Ce crystal coupled to a digital photon counter (DPC) array. This resulted in a spatial resolution of 1.7 mm FWHM, an average DOI resolution of 3.7 mm FWHM, and a CRT of 214 ps. Moreover, the possibility of using events missing the information of up to 16 out of 64 photosensor pixels is shown. This results in only a small deterioration of the detector performance.

Gamma-ray detectors based on thick monolithic scintillator crystals can achieve spatial resolutions <2 mm full-width-at-half-maximum (FWHM) and coincidence resolving times (CRTs) better than 200 ps FWHM. Moreover, they provide high sensitivity and depth-of-interaction (DOI) information. While these are excellent characteristics for clinical time-of-flight (TOF) positron emission tomography (PET), the application of monolithic scintillators has so far been hampered by the lengthy and complex procedures needed for position-and time-of-interaction estimation. Here, the algorithms previously developed in our group are revised to make the calibration and operation of a large number of monolithic scintillator detectors in a TOF-PET system practical. In particular, the k-nearest neighbor (k-NN) classification method for x,y-position estimation is accelerated with an algorithm that quickly preselects only the most useful reference events, reducing the computation time for position estimation by a factor of ~200 compared to the previously published k-NN 1D method. Also, the procedures for estimating the DOI and time of interaction are revised to enable full detector calibration by means of fan-beam or flood irradiations only. Moreover, a new technique is presented to allow the use of events in which some of the photosensor pixel values and/or timestamps are missing (e.g. due to dead time), so as to further increase system sensitivity. The accelerated methods were tested on a monolithic scintillator detector specifically developed for clinical PET applications, consisting of a 32 mm × 32 mm × 22 mm LYSO : Ce crystal coupled to a digital photon counter (DPC) array. This resulted in a spatial resolution of 1.7 mm FWHM, an average DOI resolution of 3.7 mm FWHM, and a CRT of 214 ps. Moreover, the possibility of using events missing the information of up to 16 out of 64 photosensor pixels is shown. This results in only a small deterioration of the detector performance.
Keywords: monolithic scintillator detector, positron emission tomography (PET), practical calibration, accelerated k-NN position estimation, maximum likelihood interaction time estimation (MLITE), depth of interaction (DOI), incomplete light distributions (Some figures may appear in colour only in the online journal)

Introduction
Monolithic scintillator detectors are based on a single-crystal scintillator with typical edge dimensions of 15-50 mm and a thickness of 10-25 mm, coupled to a pixelated photosensor (figure 1). For each scintillation event, the photosensor registers the light intensity on each of its pixels and acquires one or more timestamps, from which the position-, time-and energyof-interaction are derived. These detectors can achieve a spatial resolution better than 2 mm full-width-at-half-maximum (FWHM) (Li et al 2012, Cabello et al 2013, Borghi et al 2015 in combination with a coincidence resolving time (CRT) below 200 ps FWHM . In addition, they can provide depth-of-interaction (DOI) information (Ling et al 2007, Hunter et al 2009, Li et al 2010, van Dam et al 2011b, high sensitivity, and good energy resolution (Seifert et al 2012a. Monolithic scintillator detectors thus offer a unique combination of characteristics with an excellent match to the requirements of clinical time-of-flight (TOF) positron emission tomography (PET).
The position of interaction inside the crystal is estimated using analytical or statistical algorithms. Analytical algorithms range from simple center-of-gravity (COG) methods to more complex models of the expected light distribution as a function of the 3D position of interaction (Ling et al 2008, Li et al 2010. Examples of statistical models are maximum-likelihood positioning (Ling et al 2007, Hunter et al 2009, neural networks (Bruyndonckx et al 2003(Bruyndonckx et al , 2004, and k-nearest neighbor (k-NN) methods (Maas et al 2006, van Dam et al 2011a.
Statistical methods often provide the best spatial resolution. However, they require meticulous calibration procedures based on a fine sampling of the detector response over a large number of positions, for example by means of pencil beam (PB) irradiations, which translates into long calibration times. Moreover, these algorithms can be computationally demanding. Similar issues affect the methods used to estimate the DOI and time of interaction, both of which may require position-dependent calibration datasets (van Dam et al 2011b. These practical difficulties add up if multiple detectors have to be used and have so far hampered the development of clinical TOF-PET systems based on monolithic scintillator detectors. In this paper we aim to make monolithic scintillator detectors more practical for largescale application by revising the calibration procedures and estimators for position, time-ofinteraction and energy previously developed in our group, in order to reduce their complexity and time consumption and make them more robust with respect to missing data and noise. Accelerated calibration procedures and estimators for position and energy are presented in section 2, while an improved version of the time-of-interaction estimation technique is described in section 3. In section 4 a new approach is introduced to use events with incomplete information (i.e. with missing photosensor pixel values and/or timestamps), so as to improve sensitivity while softening the requirements on the detector operating parameters (dead time, temperature, etc). In these sections the problems that affected the previously used calibration procedures and estimators are briefly reviewed, after which the general description of the techniques developed to solve these issues is given. The applicability and effectiveness of the new methods are then demonstrated on a monolithic scintillator detector developed for clinical TOF-PET applications, consisting of a 32 mm × 32 mm × 22 mm LYSO : Ce crystal coupled to a digital photon counter (DPC) array. Complete descriptions of the detector and of the measurements performed to calibrate and test it are presented in section 5. The practical implementation of the calibration procedures and estimators and the performance that can be achieved with the detector are reported in section 6.

Improved calibration procedures and estimators for position and energy
2.1. x,y-position estimation Maas et al (2006) introduced the k-NN algorithm to estimate the position of interaction in monolithic scintillator detectors. This method compares the light distribution of an unknown event with the light distributions of a large dataset of reference events, i.e. events for which the position of interaction is known. The reference events most similar to the unknown events are selected and their known position of interaction is used to estimate the position of interaction of the unknown event (see section 2.1.2). For sufficiently large reference datasets, this method should approach the theoretical minimum of the misclassification probability. However, the implementation by Maas et al required very large reference datasets acquired with a PB at a finely spaced grid on the crystal surface for hundreds of different angles of incidence. A modified version was therefore introduced by van Dam et al (2011a), which required reference data acquired with a perpendicularly incident PB only, decreasing the calibration time by at least two orders of magnitude. The next implementation of the k-NN method, by Borghi et al (2015), utilized independent reference datasets along the x-and y-directions obtained with a fan beam (FB) of annihilation photons, decreasing the calibration time by at least two additional orders of magnitude. Reasonable calibration times of 2-3 h per detector were thus achieved for monolithic scintillator crystals with edge dimensions of ~2 cm.
Nevertheless, these k-NN implementations still required the calculation of the Euclidean distance of the light distribution of the unknown event to the light distributions of all reference events. Since the number of reference events may be in the order of 10 5 -10 6 and modern photosensors offer a large number of channels per detector, the resulting computational requirements can become prohibitive. Therefore, a new approach is introduced here, which is based on a coarse but fast position pre-estimation for reference and unknown events. Using this pre-estimation, the final, more accurate position estimation can be performed with the k-NN algorithm using only a small fraction of the reference light distributions, viz. those that are estimated to originate in the same region of interaction as the unknown event. Reducing the amount of reference events used in the k-NN algorithm, the method is significantly accelerated.
2.1.1. Position pre-estimation algorithms. Position pre-estimation is performed using a clustering method derived from the DOI estimation algorithm introduced by van Dam et al (2011b), which in turn was inspired by Ling et al (2007). The method uses a parameter correlated to the coordinate to be estimated and a lookup table (LUT). For the detector under study (figure 1), the x-and y-coordinates of the center of gravity of the measured light distribution (COG x and COG y ) were chosen as parameters to perform the classification in the x-and y-directions, whereas the sum of the squared pixel intensity S (SSPI) was used for the DOI: Here, n i is the light intensity on each of the light sensor pixels in the normalized light distribu- , with N pxls the total number of pixels. The parameter S can be considered as a measure of how much the light distribution is peaked. The closer the event occurs to the photosensor, the narrower the light distribution and, therefore, the higher is S, for events having the same x,y-coordinates. The COG and SSPI were selected because they are simple, robust, and computationally inexpensive parameters for the 3D position pre-estimation.
A calibration dataset containing a sufficient number of events distributed over the whole crystal volume is required to build the LUTs. Knowledge of the real positions of irradiation on the crystal of the events belonging to the calibration dataset is not necessary, therefore a flood irradiation with 511 keV photons could be used. However, given the irradiation condition used to acquire the dataset, it is necessary to know the resulting probability distribution of the energy deposition centroids, i.e. the center of gravity of the different energy deposition points for each gamma-ray interaction, over the crystal volume. This information can be obtained by means of a simple Monte Carlo (MC) simulation.
To obtain the LUT for pre-estimating the x-coordinate (LUT x PRE ) the crystal is subdivided into n x PRE regions of equal width along the x-axis (covering the whole crystal in the y-and z-directions) and the expected fractions f x x ( ) of events interacting in each region are derived from the simulated probability distribution of the energy deposition centroids. Next, the COG x values of the calibration events are calculated and sorted in ascending order, after which the n 1 x PRE − COG x values demarcating the fractions f x x ( ) of the sorted series are selected and stored in LUT x PRE . To pre-estimate the x-position of an unknown event, first its COG x value is calculated. Then, the x-region delimited by the LUT x PRE values that are below and above COG x is selected and its x-position is taken as the pre-estimated position of the unknown event. An equivalent procedure is followed for position pre-estimation along the y-axis. For DOI pre-estimation, the detector x,y area is divided into n n x y PRE-DOI PRE-DOI × regions in the x,y-plane and in n z PRE-DOI depth-ranges in the DOI-direction. For each x,y-bin, the fractions f z x y , z ( ) of events interacting in each DOI layer relative to the total number of events interacting in that x,y-bin are derived from the MC simulation. The COG x , COG y and SSPI values are calculated for each calibration event, which is then assigned to one of the x,y-bins using their pre-estimated x,y position, obtained with LUT To pre-estimate the DOI of an unknown event, the appropriate x,y-bin of the LUT DOI PRE is first selected using its pre-estimated x,y-position and then the pre-estimation of the DOI value is derived by comparing the SSPI value of the event to the SSPI values in the LUT DOI PRE for that x,y-bin.
2.1.2. Accelerated k-NN 1D method for x,y-position estimation. The more accurate, definitive position estimation is performed using an accelerated version of the smoothed k-NN 1D method described by Borghi et al (2015). This method requires two reference datasets acquired with a collimated FB of 511 keV photons, one with the FB aligned perpendicularly to the crystal x-axis (x-subset) and one with the FB rotated by 90 degrees (y-subset). For each dataset, a fixed number of full-energy events is acquired at a series of equidistant, accurately known, positions along the x-or y-axis, covering the whole crystal surface. The positions of all reference events are subsequently pre-estimated as described in section 2.1.1.
To estimate the x-coordinate of the interaction of an unknown event using the accelerated k-NN 1D algorithm, the event position is first pre-estimated with the method described in section 2.1.1. Using the pre-estimated positions, the events of the x-subset having a distance to the unknown event in the x,y-plane smaller than an optimized value (defined r xy PRE ) and a distance in the DOI-direction smaller than another value (defined r z PRE ) are selected as the reference dataset for the k-NN algorithm. Then, all the Euclidean distances of the light distribution of the unknown event to the light distributions of all the selected reference events are calculated and the k closest light distributions (nearest neighbors) are found. A 1D histogram of the x-position of the nearest neighbor events is then built and smoothed with a moving average filter n bins wide, whose dimensions are asymmetrically reduced when it approaches the edges. The x coordinate of the unknown event is estimated as the position of the maximum of the histogram. If multiple maxima are present, one of them is randomly chosen. An analogous procedure is used to estimate the y-coordinate.

DOI estimation
The more accurate, definitive DOI estimation is performed with a method similar to that used for position pre-estimation. The only difference is that the accelerated k-NN 1D method is used to estimate the x,y-positions of the events of the calibration dataset used to build the final LUT, LUT DOI , as well as the x,y-position of the unknown event whose DOI has to be estimated.
The advantage of this method over the one described by van Dam et al (2011b) is that the estimated position of the calibration events instead of their irradiation positions is used to build LUT DOI and therefore it is not necessary to perform PB irradiations to acquire the calibration events.

Energy correction for position-dependent detector response
Monolithic scintillator detectors usually show a rather homogeneous energy resolution across the whole crystal volume if they are completely covered with a highly reflective material due to the favorable light-collection conditions (Maas et al 2009, Seifert et al 2012a. Nevertheless, a simple method to estimate and correct for any remaining energy response variations across the detector volume has been developed. First, the crystal volume is divided into n n n voxels. Next, each event of a large calibration dataset distributed over the whole crystal volume is assigned to the voxel containing its position of interaction estimated using the methods described in sections 2.1 and 2.2. An energy spectrum is created for each voxel and the center position of its 511 keV photopeak is estimated using a Gaussian fit. A correction factor for each voxel is subsequently calculated as the ratio between the photopeak position of the whole detector and the photopeak position of that voxel and stored in a LUT (LUT EN ). When an unknown event is registered, first its position of interaction is estimated and then its energy is corrected using the factor of the corresponding voxel of LUT EN .

Improved calibration procedures and estimators for time of interaction
Pixelated photosensors may register a timestamp for each pixel or group of pixels, providing a maximum of N ts timestamps per event. To distinguish the N pxls photosensor regions that independently measure the amount of incident photons from the N ts regions that provide a timestamp, the former will be referred to as light-pixels and the latter as time-pixels. If each time-pixel accurately measures the time of arrival of the first scintillation photon(s) detected, the spatio-temporal distribution of the N ts timestamps can be used to correct for the optical transport times of the photons inside the crystal and to calculate more precisely the time of interaction of each gamma photon. In order to optimally perform this correction, a method for maximum-likelihood interaction time estimation (MLITE) was proposed by van . Here, this method has been modified in order to make the calibration procedure faster and more practical, as well as to improve the timing performance by selecting only the timestamps with reliable information and using a non-parametric fit for the probability density functions.
The calibration procedure is based on the acquisition of a large calibration dataset of events distributed over the whole crystal volume in coincidence with a fast reference detector. No prior knowledge on their positions of irradiation is needed. We first attempt to remove premature timestamps (e.g. generated by dark counts), as well as late timestamps that do not contain relevant information (Seifert et al 2012b) from the set of N ts timestamps acquired for each scintillation event. This operation makes use of the fact that timestamps acquired in the first part of the scintillation pulse (e.g. in the first ~5 ns for a LYSO crystal) should occur close to each other in time. To perform this operation it is important that no significant time skews are present between the different time-pixels, e.g. due to electronic jitter or readout delays. First, all timestamps of each scintillation event are sorted in order of acquisition (i.e. earliest timestamp in first position, etc) and the timestamps which are followed by another timestamp in a time period t ts-seq start ∆ are selected. The first valid timestamp is determined as the earliest selected one which is followed by another selected timestamp. Similarly, the last valid timestamp is determined sorting the timestamps in inverse order of acquisition, selecting the timestamps preceded by another timestamp in a time period t ts-seq end ∆ and finding the last selected timestamp preceded by another selected timestamp. All timestamps registered in between the first and last valid ones are also considered valid. To make sure that no useful timestamps are discarded, it is possible to define an additional time window before the first valid timestamp ( t t ts start ts-seq start ∆ > ∆ ) and after the last one ( t t ts end ts-seq end ∆ > ∆ ), accepting also the timestamps acquired in these windows. The length of the four time windows has to be optimized for each type of detector.
Once the valid timestamps are selected, the crystal is divided into n n n positions of interaction. The interactions in a given voxel are registered by time-pixel i with certain delays, denoted as first photon detection delays (FPDDs), which are determined by the scintillation pulse shape and the optical transport in the crystal. The FPDD probability distribution function (PDF) is estimated for each combination of voxel and time-pixel using the calibration events assigned to that voxel. First, the differences between the valid timestamps acquired by the time-pixel and the corresponding timestamps of the reference detector are calculated and then these values are used to estimate the FPDD PDFs with kernel-density estimation (KDE). This non-parametrical method reduces the influence of statistical fluctuations on the experimentally measured PDFs and, especially at the beginning of the PDFs, preserves more closely their shape compared to the previously proposed method, which employed a fit with an exponentially modified Gaussian . Finally, knowing the geometry of the calibration setup, the difference in the travel times of the paired annihilation quanta to the interaction points in the test and reference detectors can be determined and the zero of the FPDD distributions can be set at the moment of interaction inside the monolithic crystal. Thus, a group of N ts PDFs are obtained per crystal voxel, denoted as p t x y z i , , , , with t the delay relative to the time of interaction in the monolithic scintillator detector, x y z , , ( ) the coordinates of the voxel center, and i the index of the sensor time-pixel.
Once an unknown event is detected and its position of interaction is estimated, the event is assigned to the corresponding MLITE voxel and its timestamps are selected using the same validity conditions used to build the FPDD PDFs. Then, the time of interaction t int is estimated using the FPDD distributions of that voxel and a maximum likelihood algorithm. That is, the likelihood of having a set of valid timestamps t t , N 1, given an interaction time t int is defined as: The most likely time of interaction t int is obtained by finding the maximum of the likelihood function:^{

Position, energy, and time of interaction estimators for events with missing information
In part of the events not all light-pixels or time-pixels are acquired, e.g. a pixel could be in dead time due to the readout or recharge process resulting from an earlier trigger caused by a dark count or a previous scintillation event. If such events with missing information could be used in a practical way to still obtain accurate time, energy and position information, a significant improvement in detector sensitivity could be achieved. Moreover, the method could enable the use of detectors that have some malfunctioning or broken pixels. In fact, a complete light distribution generally is not essential if statistical methods are used to estimate the position of interaction, since they can be applied to only that part of the information that is available . However, the photon count on the missing pixels is still important for correctly estimating the energy deposited and/or when using analytical methods to estimate the position of interaction, such as the position pre-estimation methods and the DOI estimation algorithm described in sections 2.1.1 and 2.2. In such cases, missing pixel values can be estimated using analytical methods (e.g. interpolation/extrapolation techniques) based on calibration datasets , or statistical methods. Here, a new method based on average light distributions and the k-NN algorithm is introduced.
To calculate the average light distributions, a procedure similar to that described in section 2.3 to build LUT EN is followed. The monolithic crystal is subdivided into n n n x y z av-ld a v-ld av-ld × × equal voxels and each event of a large calibration dataset containing only complete light distributions is assigned to a voxel depending on the estimated position of interaction. All light distributions are subsequently normalized to the same total value and the average light distribution for each voxel is calculated.
When an incomplete light distribution is acquired, first a coarse energy discrimination is performed using different energy spectra built with events having different numbers of missing pixels, i.e. spectra built using only events with one missing pixel, two missing pixels, etc. Then the average light distribution that is most similar to the incomplete one is found using the k-NN algorithm. First, all pixels that are missing in the incomplete light distribution are set to zero also in the average light distribution dataset. Then, all distributions are normalized to the same total photon count and the one that is most similar (i.e. having the smallest Euclidean distance) to the incomplete distribution is selected. Finally, the values of the missing pixels of the incomplete light distribution are estimated as: with p inc-ld the pixel values of the incomplete light distribution before normalization, p av-ld the pixel values of the selected average light distribution, j N 1, , miss pix = … the indices of the missing pixels, and i N 1, , acq pix = … the indices of the acquired pixels. The estimated photon counts on the missing pixels are used to complement the measured light distribution, from which the total energy of the event is calculated. The complemented light distribution is also used for the position pre-estimation used to select reference events for the accelerated k-NN 1D algorithm. Next, the x,y-position of interaction is determined using the accelerated k-NN 1D method described in section 2.1.2. Whereas the complemented distribution is used for position pre-estimation, only the pixels that were acquired during the measurement are used for the definitive position estimation. That is, the pixels missing in the measured light distribution are also set to zero in the pre-selected reference light distributions and all distributions are normalized to the same total photon count. Finally, the complemented light distribution is used once more, this time to estimate the DOI using the method detailed in section 2.2.
If an event is missing part of the timestamps, the estimation of the time of interaction with the MLITE method is performed following the approach described in section 3, using only the available timestamps to define the likelihood function given in equation (2).

Monolithic scintillator detector and reference detectors
The monolithic scintillator detector tested in this work is based on a LYSO : Ce crystal (Crystal Photonics) with a base area of 32 mm × 32 mm, a thickness of 22 mm, and polished surfaces. The four lateral faces were covered with a specular reflector foil (Vikuiti ESR, 3M), whereas the top face was covered with Teflon tape. The 32 mm × 32 mm back surface was coupled to a DPC array developed by Philips Digital Photon Counting (PDPC), model DPC-3200-22-44 (figure 1), using a removable transparent silicone gel (Sylgard 527, Dow Corning).
The DPC array measures 32.6 mm × 32.6 mm and consists of 4 × 4 independent sensors (dies) at a regular pitch of 8 mm, each die containing 2 × 2 light-pixels whose dimensions are 3.2 mm × 3.8775 mm. Each pixel comprises 3200 microcells and is subdivided into 2 × 2 equal subpixels. Each microcell consists of a single photon avalanche photodiode (SPAD) and circuitry to actively quench and recharge the SPAD, read out its status, and selectively enable or disable it. This last feature makes it possible to disable the SPADs that show the highest dark count rates. In this work, the noisiest 5% of microcells on each die were switched off. More details about the sensor architecture can be found in Frach et al (2009Frach et al ( , 2010, Schulze (2014) and Schaart et al (2016).
The subpixels are used to set a statistical threshold for triggering the acquisition sequence of a die at the beginning of a scintillation event (Frach et al 2009, Tabacchini et al 2014. Trigger level MT =1 (PDPC notation, see Frach et al (2009)) was chosen in this work. With this setting, a die is triggered by the first SPAD firing, upon which a single timestamp is acquired and photon counting is continued during a user-defined validation time (set to 40 ns in this work). At the end of this interval, the die checks whether a higher-energy threshold criterion is met, so as to discard events triggered by dark counts. This validation threshold was set such that at least one pixel should have at least one photon counted on each subpixel (DPC validation threshold notation: '0 × 7F : AND'). If this condition is not met, the die undergoes a fast (~20 ns) recharge and reset cycle. Otherwise if it is met, the die continues counting for a userdefined integration time (set to 165 ns in this work) and then carries out a readout sequence (680 ns). At the end of this sequence it outputs the photon count of each of the four pixels as well as a single timestamp. Following the definitions given in section 3, a die thus corresponds to a time-pixel and contains four light-pixels.
The DPC array is equipped with neighbor logic, which makes it possible to force the readout of neighboring dies or the entire array after the validation of any single die. In this work, the option to force the readout of the entire array was activated.
Several reference detectors similar to the detector under test were assembled to acquire events in coincidence with it. Each detector was built using a different DPC sensor. Two detectors were based on two polished LSO : Ce crystals (Agile Engineering Inc., Knoxville, TN, USA) with dimensions of 32 mm × 32 mm × 20 mm and 16 mm × 16 mm × 20 mm. Three other detectors were assembled using Ca-codoped (0.2% in the melt) LSO : Ce crystals (Spurrier et al 2008) (produced at the Scintillation Materials Research Center, University of Tennessee, and provided by Agile Engineering Inc., Knoxville, TN, USA) with dimensions of 3 mm × 3 mm × 5 mm and a single polished 3 mm × 3 mm face that was coupled to a pixel of a DPC array. The small reference crystals were covered on all the faces with specular reflector foil.
The readout and data acquisition of all detectors was performed using field-programmable gate array based electronic boards and an accompanying computer program provided by PDPC (Schulze 2014).

Experimental setup
A slightly modified version of the experimental setup described in Borghi et al (2015) was used to perform the calibration and test measurements. Briefly, the setup is based on a paired collimator system that is placed in between the detector under test and a reference detector. It consists of a central tungsten cylinder containing a 22 Na point-source (Ø 0.5 mm, ~3.5 MBq, IDB Holland BV) and two sets of collimators used to obtain either a Ø 0.5 mm PB or a 0.5 mm wide FB of annihilation photons. The collimator system can be rotated by an angle of 90° so that the FB can be aligned perpendicularly to the x-or y-axis of the detector under test. This detector is mounted on two perpendicular linear stages driven by stepper motors (Physics Instruments, M-403.42S stages with C-663 controllers), which precisely move and position the detector in the plane perpendicular to the collimator axis.
The entire setup is placed inside a light-tight temperature chamber (Weiss WT 450/70). A Peltier element is mounted on the back of each detector to control its temperature with 0.1 °C precision. The working temperature of all detectors was set to −25 °C to reduce the dark count rate of the DPC arrays.

Data acquisition
Four different sets of events were acquired for calibrating and characterizing the position estimation response of the detector. These are referred to as the FB, perpendicular PB, side PB, and angular PB datasets. The 32 mm × 32 mm × 20 mm LSO : Ce detector was used as a reference detector to acquire the FB dataset while the 16 mm × 16 mm × 20 LSO : Ce detector was used for the other measurements. A count rate profile was performed prior to each measurement in order to find the center of the crystal with respect to the beam position. The use of the different datasets is detailed in section 6.
The FB dataset was obtained with the FB perpendicularly incident on the crystal front surface. First, the FB was aligned perpendicularly to the crystal x-axis and the x-subset was acquired collecting 12 800 full-energy events (see section 5.4 for data pre-processing) at a series of reference positions spaced 0.25 mm apart along the x-direction. Subsequently, the FB was rotated by 90° and the y-subset was acquired in an analogous way. The total acquisition time was about 6-7 h.
The perpendicular PB dataset was acquired with a perpendicularly incident PB at a grid of reference positions with a pitch of 0.25 mm covering the entire crystal front surface. At each point, 100 full-energy events were acquired. The total acquisition time was about 8 d.
The side PB dataset was collected by perpendicularly irradiating a crystal side surface with the PB at a grid of reference positions with a pitch of 1 mm in both the x-and z-(DOI-) directions, covering all possible DOIs and spanning from x = −15.5 mm to x = 0.5 mm. At each position ~4500 full-energy events were acquired.
Finally, the angular PB dataset was obtained by irradiating half of the crystal front surface at a 0.25 mm pitch reference grid with the PB at a 60° angle of incidence (figure 2). The PB was aligned parallel to the crystal y,z-plane and only the grid positions for which the entry and exit points of the line of irradiation intercepted the front and back face of the crystal were acquired. About 60 full-energy events were acquired per position.
For the calibration and characterization of the time response of the monolithic scintillator detector, the same preliminary measurements described in van Dam et al (2013) were performed. Briefly, before the crystal was mounted, the electronic time skews between the dies of the DPC array were determined with a pulsed laser (Hamamatsu PLP-04 laser, wavelength 633 nm, average pulse duration 50 ps). Furthermore, the CRT of each of the three possible combinations of 3 mm × 3 mm × 5 mm Ca-codoped LSO : Ce detectors was measured. This permitted us to calculate the single-detector time resolution for each detector. The best reference detector, with a single-detector time resolution of ~89 ps FWHM and an estimated CRT for two identical detectors (CRT ref ) ~125.9 ps, was chosen as a reference detector for the flood irradiations described below.
Three different datasets were then acquired for the timing calibration and characterization of the monolithic detector. To this end, the uncollimated 22 Na point source was placed in between this reference detector and the test detector. The distance between the source and the reference detector was kept at 12 mm, while the monolithic crystal was placed at 22 cm, 24 cm and 26 cm from the source so as to irradiate the whole monolithic crystal surface uniformly. About 4 million and 2 million full-energy coincidence events were acquired at 24 cm and the other distances, respectively. The resulting datasets will be referred to as the 22 cm flood irradiation (22 cm FI), 24 cm FI, and 26 cm FI datasets. A little less than an hour was required to acquire 1 million events during flood irradiation.

Data pre-processing
The first pre-processing step for all measurements was to select the events in which no pixel values and timestamps were missing (corresponding to about 65% of all events acquired with the detector settings and operating conditions used in this work) and for which the total energy deposited fell within the full-width-at-ten-maximum (FWTM) of the 511 keV photopeak. For timing measurements, the energy condition was also applied to the reference detector.
A correction to obtain a homogeneous response from the light sensor pixels was subsequently applied so as to obtain more accurate results from the analytical position estimation methods (i.e. the methods used for position pre-estimation and DOI estimation). The perpendicular PB and the FB datasets, which contain events evenly distributed over the crystal surface, were used to build LUTs for uniformity correction (LUT UC ). The assumption is that the average light distribution over all the events in both datasets should have the same value for all DPC pixels. Therefore, the LUT UC value for each pixel is defined as the ratio between the average value of all pixels and the value of that same pixel in the average light distribution. Since the light-pixels of the DPC array already have a highly uniform response, this correction was mostly performed to remove any effects due to the non-homogeneous quality of the optical coupling and/or possible deficiencies in the crystal wrapping. If necessary, this uniformity correction could be performed regularly, e.g. based on simple flood irradiations, to make the detector response more stable without having to repeat the entire k-NN calibration.

Results and discussion
6.1. Position estimation 6.1.1. Spatial resolution at perpendicular incidence. The perpendicular PB dataset (see section 5.3) was used to test the spatial resolution in the x-and y-directions. The positions of all events were estimated using the accelerated k-NN 1D method (see section 2.1) as well as the standard k-NN 1D method (Borghi et al 2015), which calculates the Euclidean distances between the light distributions for all the events in the reference datasets and not for only part of them, selected according to their pre-estimated position of interaction. The FB dataset was used as the reference dataset for both methods. The parameter k was set to 100 and 200 for the standard k-NN 1D and the accelerated k-NN 1D algorithms, respectively. In both cases the smoothing filter had a width of n 5 = . The FB dataset (the x-and y-subsets together) was also used to create the position pre-estimation LUTs, with n n 32  the misalignment correction procedure presented in Borghi et al (2015) was used to correct all results in this section. Table 1 shows the parameters characterizing the average spatial resolution obtained with both k-NN methods, obtained using all events in the PB dataset. The FWHM and FWTM parameters are derived from the 2D point-spread function (2D PSFs), defined as the normalized 2D histogram of the differences between the estimated x,y-positions of interaction and the known x,y positions of irradiation. The FWHM and FWTM are defined, respectively, as the full-width-at-half-and tenth-maximum of the cross-sections along the x-and y-directions of the 2D PSF across its maximum (figure 3, left).
Slightly better FWHM and FWTM values are obtained with the standard 1D k-NN method; however, since the 2D PSFs do not have a Gaussian shape this is only due to the difference in height of their central bins and does not represent a real difference in the positioning performance of the methods. For this reason we consider the other parameters reported in table 1 to be more comprehensive. The r 50% and r 90% values are derived from the normalized cumulative distribution functions (CDFs) of the x-, y-and total errors (figure 3, right), the total error being the length of the error vector. These parameters are defined as the values at which the CDF exceeds 0.5 (i.e. the median error) and 0.9, respectively. The mean absolute error (MAE) is defined as the mean absolute value of the (x-, y-or total) error. A comparison of the CDFs of the two k-NN methods in figure 3 (right) shows that their positioning accuracy is practically equivalent, in fact with slightly better r 50% , r 90% , and MAE values for the accelerated 1D k-NN method, as reported in table 1.
The   Figure 4 left-right respectively shows the average MAE value and bias magnitude (calculated as the average Euclidian norm of the bias vectors) in 1 mm × 1 mm regions on the detector front surface, each region containing 16 irradiation positions.
As expected, the spatial resolution is best at the center of the crystal and starts to degrade at about 4-5 mm distance from each edge, for the coordinate perpendicular to that edge (Li et al 2012, Seifert et al 2012a. This effect simultaneously occurs for both coordinates in the corners, where the MAE consequently is worse. Nevertheless, the MAE is <3 mm at practically all positions. The bias is negligible (<0.5 mm) except in ~3 mm wide regions near the crystal edges. It becomes substantial (>1.5 mm) only at about 1 mm from the edge, which is attributed to the truncation of the error distributions at the edges.
It has to be noted that all spatial resolution results also include the contributions due to the finite beam diameter, which in this work can be considered negligible, and the errors determined by events which underwent Compton interaction(s) inside the crystal. Compton  scattered events can travel a significant distance in thick scintillation crystals and therefore contribute to the tails in the error distributions (PSFs and CDFs) (Maas et al 2010).
For the accelerated 1D k-NN method, an average of ~8900 reference events were pre-selected to estimate the x-or y-coordinate of each unknown event, while ~1.6 million reference events were necessary with the standard method. A reduction by a factor of ~200 in the number of distance calculations, which is the most computing-intensive operation in the k-NN algorithm, has thus been achieved without any degradation of the spatial resolution. The computation time per unknown event using a non-optimized MATLAB implementation of the accelerated 1D k-NN method on a single core was ~5 ms, a significant improvement compared with the value of ~1 s needed with the standard method. This time could still be decreased substantially by replacing the MATLAB script by an optimized code.
6.1.2. DOI resolution. The whole FB dataset was used to create LUT DOI (see section 2.2), with n n 16 their x,y-position was estimated with the accelerated 1D k-NN method using the leave-one-out technique (Maas et al 2009), i.e. removing the event whose position is being estimated from the reference dataset used for k-NN. The side PB irradiation was used as a test dataset to determine the performance of the DOI estimation. Table 2 shows the results. The FWHM and FWTM values were calculated from the 1D distribution of the DOI errors, whereas r 50% and r 90% were obtained from the normalized error CDF. Figure 5 shows the DOI resolution and DOI bias as a function of the z-position of interaction. These plots were obtained by calculating the MAE and the mean error, respectively, for all events acquired at each DOI in the side PB dataset. As expected, the DOI resolution is best near the photosensor and degrades in the upper half of the crystal. The bias is less than 1.5 mm for almost all z-positions, becoming significant in the 3-4 mm near the front crystal surface and the 1-2 mm near the photosensor only, due to truncation of the error distributions. It is noted that these results were obtained with a homogeneous side irradiation; in the case of front irradiation, with more events interacting in the front half of the detector, the average DOI resolution would be slightly worse.
6.1.3. Spatial resolution at non-perpendicular incidence. The angular PB dataset was used to test the spatial resolution at an angle of incidence similar to that occurring in a clinical scanner for events originating at the edge of the radial field of view (FOV). The positioning error was determined for each event by estimating the position of interaction (x, y, z) and calculating its distance to the true line of response (LOR); see figure 2. The distance between the estimated x-position and the x-coordinate of the plane parallel to the y,z-plane that contains the LOR is denoted as the error in the x-direction. The distance between the estimated y,z-position and the Table 3. Spatial resolution for events incident at an angle of 60° on the crystal surface, obtained using the DOI estimated with the method presented in section 2.2 ('DOI') and a fixed DOI value of 7.5 mm ('no DOI'). LOR in the plane parallel to the y,z-plane that contains the LOR is called the y′-error. The resolution measures defined in section 6.1.1 were calculated on the basis of the x-and y′-errors.
To assess the importance of DOI estimation, the calculation was repeated after replacing z by a fixed DOI value of 7.5 mm, equal to the average DOI at 60° incidence. The results are summarized in table 3, while the 1D PSFs (i.e. the projections of the 2D PSFs on the x or y axis) are shown in figure 6. A small resolution degradation is noticeable for the x-coordinate. This could be explained by a small misalignment of the irradiation positions during the non-perpendicular irradiation, which is noticeable in figure 6 as a small offset between the error histograms. Also, the angular PB dataset on average has a smaller DOI than the perpendicular PB dataset and the x,y spatial resolution is worse in the front part of the crystal, which is further away from the DPC arrays.
As for the y′-error, the lower accuracy of the DOI estimation compared to the x,y estimation results in less accurate LOR positioning for non-perpendicular LORs. However, the MAE in the y′-direction is only ~17% higher compared to perpendicular irradiation, which is a large improvement over the ~75% degradation that occurs if the average DOI is used for LOR positioning.

Energy resolution
The FB dataset was used to create LUT EN (see section 2.3), with n n 16 x y en en = = and n 4 z en = .
The perpendicular PB dataset was subsequently used to test the influence of the position-dependent energy correction on the detector energy resolution. Only events belonging to the noncorrected 511 keV photopeak were used (see section 5.4). The energy resolutions before and after correction were determined with Gaussian fits. The center position of the 511 keV photopeak as a function of the x,y,z-position of interaction is shown in figure 7. At all DOI, it is possible to notice a region around (x, y) ~ (2,−10) in which the detector shows a reduced energy response, probably due to small defect(s) in the optical coupling. In the DOI layer between 16.5 mm and 22 mm, which is the closest to the DPC array, it is also possible to observe a different response between the regions above the dead areas of the photosensor ( y = −16,−8, 0, 8, 16) (figure 1), where the energy response is lower, and the regions centered above the DPC dies ( y = −12,−4, 4, 12), where the energy response is higher. This results in a difference in the photopeak position of up to ~200 cells. Despite this variation, the average energy resolution in the energy spectrum of the complete dataset without any correction still equals ~10.25% FWHM (figure 8). If the position-dependent energy correction is applied, an energy resolution of ~9.9% FWHM is achieved.

Time resolution
The timestamps of all events in datasets 22 cm FI, 24 cm FI, and 26 cm FI were corrected for the electronic time skews between the DPC dies (see section 5.3) and subsequently preselected using the method described in section 3. The optimized values for the time windows pre-computed on a temporal grid with a spacing equal to 1 DPC time-to-digital converter bin (10 ns/2 9 ≅ 19.5 ps) and stored in a LUT. Since the absolute time is irrelevant for determining the time resolution, the zero of the FPDD PDFs was left uncorrected for the irradiation geometry. The remaining half of the 24 cm FI dataset as well as the 22 cm FI and 26 cm FI datasets were used as test datasets. The time of interaction of each event was estimated using MLITE, with the values of the likelihood function t L t int ( ) | calculated on the same time grid as the FPDDs. The CRT for two identical monolithic scintillator detectors in coincidence was then calculated for each dataset as: where CRT exp is the experimentally measured CRT and CRT 2 ref 2 / is the squared value of the single-detector time resolution of the reference detector (see section 5.3). A deterministic timing method was also tested to assess the time resolution that can be achieved without any calibration procedure, apart from the measurement of the time skews between the DPC dies that remains necessary. Specifically, the average value of the first two timestamps registered by the monolithic detector was used as t int .
With MLITE, CRT exp is found to be 176 ps FWHM, resulting in a CRT of 214 ps FWHM for two monolithic scintillator detectors in coincidence. With the deterministic method, CRT exp equals 192 ps FWHM, resulting in a CRT of 241 ps FWHM. MLITE thus provides a ~10% improvement of the CRT. The deterministic method still provides a remarkable CRT for a crystal of these dimensions; it has to be noted, however, that this result could only be achieved with a strict selection of reliable timestamps (see section 3). The computation time needed with a MATLAB implementation of MLITE running on a single core was ~1 ms per event, which could be further decreased with an optimized code.
The , respectively, the value found for the LYSO : Ce monolithic detector is in excellent agreement with what could be expected from the previous measurements with Ca-codoped crystals. Figure 10 shows the three TOF-difference spectra derived from datasets 22 cm FI, 24 cm FI, and 26 cm FI. The CRTs are 175.5 ps FWHM, 176 ps FWHM, and 175 ps FWHM, respectively. The mean times of interaction at 24 cm and 26 cm, relative to the measurement at 22 cm, are 60.3 ps and 135.2 ps, corresponding to distances of 18.1 mm and 40.5 mm, respectively. The differences with the expected distances of 20 mm and 40 mm, respectively, can be explained by the reproducibility of the source position within the setup, which is estimated to be in the order of 1-2 mm.

Events with incomplete light distributions
To compare the results obtained with the same statistics, the datasets with complete light distributions that were used in sections 6.1-6.3 were artificially modified to determine the detector performance for events with incomplete light distributions. That is, four additional  versions were created for each dataset, in which the data of 1, 2, 3, or 4 randomly chosen dies (viz., four light-pixel values and one timestamp per die) were deleted. Random deletion is considered justified since the monolithic scintillator detector is operated with full neighbor logic, so dies are expected to be missing only because of dead times caused by dark counts, with essentially equal probability per die.
The FB dataset was used to calculate the average light distributions (section 4). The methods described in section 4 were used to process the resulting four versions of the perpendicular-PB and side-PB datasets. The analysis used in sections 6.1.1, 6.1.2 and 6.2 was then repeated to determine the spatial, DOI, and energy resolution as a function of the number of missing dies. Similarly, the analysis in section 6.3 was repeated on four versions of the flood irradiation datasets with different numbers of missing dies.
To put these results into perspective, a weighted average of the PSFs, CDFs, energy spectra, and TOF-difference spectra obtained with a variable number of missing dies was also calculated using as weights the fractions of events acquired with up to four missing dies during real measurements. The resulting distributions were used to estimate the performance that would be achieved for the detector settings and operating conditions used in this work if those distributions were accepted. The percentages of the events in which the data from 16, 15, 14, 13, and 12 dies were acquired are 68%, 18%, 6.5%, 2.5%, and 1%, respectively, i.e. ~96% of the total registered events.
6.4.1. Spatial resolution at perpendicular incidence. The parameters that characterize the spatial resolution in the x-and y-directions for perpendicularly incident events with up to four dies missing are reported in table 4. The corresponding error CDFs are shown in figure 11. A small degradation is noticeable for each additional missing die, but total degradations of only ~10% and ~15% are found for the MAE and r 50% , respectively, when four dies are missing in each event. For the acquisition conditions used in this work (AC), a degradation in the spatial resolution < 1.5 % would be obtained if events with up to four missing dies were accepted. 6.4.2. DOI resolution. The average DOI resolution as a function of the number of missing dies is reported in table 2. Figure 5 shows the DOI resolution and the DOI bias for incomplete light distributions as a function of DOI. The DOI resolution deteriorates by up to ~20% (MAE and r 50% ) at four missing dies and ~2% for the acquisition conditions used in this work, which is slightly more than the degradation observed in the x-and y-directions. This may be due to the fact that missing information has a twofold influence on the DOI estimation: it affects the x,y-position estimation as well as the SSPI value used to look up the DOI in LUT DOI .  It can be noticed that the estimated DOI of events with missing dies tends to be biased towards the crystal front surface, presumably because the algorithm that estimates the missing pixel values has more difficulty in recovering light distributions in which bright pixels are missing.
6.4.3. Energy resolution. The energy resolutions for incomplete datasets with and without estimated values of the missing pixels are reported in table 5. A degradation of ~55% is observed for the dataset with four missing dies without correction. This is reduced to <10% when using the estimated pixel values, resulting in an energy resolution of 10.6% FWHM, quite an acceptable value for LYSO : Ce. A degradation <1% would be obtained for real measurements accepting events with up to four missing dies. Figure 12 shows the energy spectra obtained with events having two and four missing dies, with and without complemented pixel values, together with the spectrum obtained with events having complete light distributions. It can be noticed that the correction algorithm accurately restores the position of the photopeak and reduces the tail on the left side of the peak.

Time resolution.
The CRT for events with missing information is shown in table 6. In the case of four missing timestamps, a deterioration of ~16% is observed with MLITE as well as with the average of the first two timestamps. In this case, MLITE still achieves an excellent CRT of 250 ps even though 25% of the timing information is missing. A degradation ~1.5 % would be obtained if events with up to four missing dies were accepted during the measurements.

Conclusions
New, time-efficient calibration procedures as well as position-, DOI-, energy-, and time-ofinteraction estimators for monolithic scintillator detectors were tested experimentally on a detector based on a 32 mm × 32 mm × 22 mm LYSO : Ce crystal and a DPC array. About 9 h were needed to fully calibrate the detector using a ~3.5 MBq 22 Na point-source. The possibility of avoiding PB irradiation, which was previously needed for time and DOI calibration and took about 8 d, reduces the calibration time by a factor of at least 20 (see section 5.3). Since the detector count rate during FB irradiation was lower than the maximum detector count rate, this time could still be reduced substantially if a stronger source were used. Moreover, considering that the new methods require FB and flood irradiations only, they open up the possibility of calibrating many detectors in parallel.
The accelerated position estimator based on 1D k-NN classification with pre-selection of the reference events reduces the calculation time by a factor of ~200 without noticeably reducing the spatial resolution. An average spatial resolution of 1.7 mm FWHM/1.6 mm MAE is achieved, which is better than or comparable to the best resolutions presented so far for monolithic scintillator detectors with a thickness of 15-20 mm (Li et al 2012, Borghi et al 2015. Additionally, a position-dependent analysis of the spatial resolution and bias showed that monolithic scintillator detectors can achieve quite homogeneous performance over almost all of their volume and that edge effects are reduced using crystals with a low aspect-ratio.
To our knowledge, this is the first time that the experimental DOI performance of a monolithic scintillator detector thicker than 10-15 mm is presented. A MAE better than 3 mm was obtained in the 22 mm LYSO crystal and it was demonstrated that this considerably improves the accuracy of LOR positioning at 60° incidence. This is expected to make the reconstructed image resolution of a typical clinical TOF-PET scanner much more homogeneous throughout its FOV.
A method to correct for the position-dependent energy response of monolithic scintillator detectors was shown to improve the energy resolution from 10.25% FWHM to 9.9% FWHM. This correction would not seem essential, as the uncorrected energy resolution is close to what is achievable with LYSO : Ce already. However, the method was shown to enhance detector performance and, since it can be easily implemented, it could be applied in systems made of many detectors to improve, e.g. the robustness to imperfections occurring during detector manufacturing.
A CRT of 214 ps FWHM was achieved with the MLITE method. This is in excellent agreement with previous results obtained with 20 mm thick monolithic crystals made from Ca-codoped LSO : Ce , taking into account the difference in the scintillation decay constants. This suggests that MLITE compensates for the optical photon transport similarly well, even though the present crystal is considerably wider than the previous ones.
A new, statistical method was introduced to estimate the photon count on missing photosensor pixels, which was tested on events in which the data from up to 4 out of 16 DPC dies (corresponding to 4, 8, 12, or 16 out of 64 light-pixels) was artificially removed. The results demonstrate that this method can restore sufficient information to profitably use the deterministic position pre-estimation algorithm of the accelerated 1D k-NN method while maintaining the robustness to missing information typical of statistical position estimation methods . The degradation of the DOI resolution is slightly more pronounced, although parallax errors can still be corrected adequately with up to four missing dies. Finally, the CRT obtained with MLITE showed limited deterioration for events with up to 4 out of 16 timestamps missing. A degradation of all the performance parameters <2 % would be observed if events with up to four missing dies were accepted during measurements with the detector settings and operating conditions used in this work, obtaining an improvement in the single-detector sensitivity of about 40% compared with the case in which only complete events are accepted.
In conclusion, this paper shows that monolithic scintillator detectors can be calibrated and operated in a practical way that should allow their use in systems consisting of many detectors. The spatial-, DOI-, energy-, and time-resolution achieved with a 32 mm × 32 mm × 22 mm LYSO : Ce crystal coupled to a DPC array shows an excellent match with the requirements of whole-body clinical TOF-PET imaging.