Improved image quality using monolithic scintillator detectors with dual-sided readout in a whole-body TOF-PET ring: a simulation study

We have recently built and characterized the performance of a monolithic scintillator detector based on a 32 mm  ×  32 mm  ×  22 mm LYSO:Ce crystal read out by digital silicon photomultiplier (dSiPM) arrays coupled to the crystal front and back surfaces in a dual-sided readout (DSR) configuration. The detector spatial resolution appeared to be markedly better than that of a detector consisting of the same crystal with conventional back-sided readout (BSR). Here, we aim to evaluate the influence of this difference in the detector spatial response on the quality of reconstructed images, so as to quantify the potential benefit of the DSR approach for high-resolution, whole-body time-of-flight (TOF) positron emission tomography (PET) applications. We perform Monte Carlo simulations of clinical PET systems based on BSR and DSR detectors, using the results of our detector characterization experiments to model the detector spatial responses. We subsequently quantify the improvement in image quality obtained with DSR compared to BSR, using clinically relevant metrics such as the contrast recovery coefficient (CRC) and the area under the localized receiver operating characteristic curve (ALROC). Finally, we compare the results with simulated rings of pixelated detectors with DOI capability. Our results show that the DSR detector produces significantly higher CRC and increased ALROC values than the BSR detector. The comparison with pixelated systems indicates that one would need to choose a crystal size of 3.2 mm with three DOI layers to match the performance of the BSR detector, while a pixel size of 1.3 mm with three DOI layers would be required to get on par with the DSR detector.


Introduction
The spatial resolution of current, clinical time-of-flight (TOF) positron emission tomography (PET) scanners is largely determined by the physical size of the individual crystals, which typically are 4-5 mm wide and 20-30 mm thick. Recent Monte Carlo studies have shown that further improvement of the spatial resolution of clinical TOF-PET systems is desirable, as it can improve image quality if a sufficient number of events are acquired (Surti et al 2013, Thoen et al 2013. In particular, it was demonstrated that spatial resolution affects small lesion detectability and lesion uptake estimation, two important tasks in oncology.
The main challenge when designing high-resolution detectors for clinical PET is the ability to improve the spatial resolution without sacrificing sensitivity, time resolution, and energy resolution. In pixelated detectors, spatial resolution improvements can be achieved by reducing the size of the crystal elements. For instance, a detector module with a 2 mm crystal pitch and 22 mm thickness has been presented for potential use in clinical systems (Marcinkowski et al 2014), while a complete high-resolution whole-body PET system was recently built, using detector blocks with 2.35 × 2.35 mm 2 LYSO crystals (Wai-Hoi et al 2015). Decreasing the pixel size even further is usually considered for small-animal scanners only (Yamamoto et al 2013, Nishikido et al 2015, Ko et al 2016.
In general, the light collection process becomes inefficient for crystals with a very highaspect ratio (i.e. the ratio between crystal thickness and width), so maintaining good energy and time resolution can be challenging (Gundacker et al 2014, ter Weele et al 2015. An interesting alternative to pixelated detectors is the monolithic scintillator detector, in which a continuous scintillator (typically a cuboidal crystal) is coupled to an array of photosensors. The position of interaction of the detected gamma rays is estimated from the light intensity distributions measured on the photosensor elements, typically using a statistical positioning algorithm. Many positioning algorithms have been proposed in literature to estimate both the transversal coordinates (Bruyndonckx et al 2004, Miyaoka et al 2010, van Dam et al 2011a, Marcinkowski et al 2016 and the depth-of-interaction (DOI) (Lerche et al 2005, Ling et al 2007, van Dam et al 2011b. Monolithic scintillator detectors have already been shown to achieve an excellent combination of spatial resolution, coincidence resolving time (CRT), and detection efficiency (Seifert et al 2012, van Dam et al 2013. Until recently, the need for sophisticated calibration procedures has been considered the main limitation in the practical use of monolithic scintillator detectors. However, we have recently presented practical methods for the efficient calibration of such detectors . Using these methods, we presented a complete characterization of a monolithic scintillator detector based on a 32 mm × 32 mm × 22 mm LYSO:Ce crystal and a digital silicon photomultiplier (dSiPM) array coupled to the back surface of the crystal, i.e. the surface facing away from the scanner axis (Borghi et al 2016b).
With the aim of further improving the detector performance for the next generation of TOF-PET scanners and to push the technology closer to its physical limits, a monolithic scintillator detector with dSiPM arrays coupled to both the crystal front and back surface in a so-called dual-sided readout (DSR) configuration was subsequently developed and fully characterized (Borghi et al 2016a). The idea of using a single crystal read out by two photodetectors, on both sides, was originally proposed in the context of conventional pixelated detectors, mainly to provide DOI capability (Moses and Derenzo 1994) and/or to improve timing performance . In monolithic scintillator detectors, however, the DSR configuration has the additional advantage that the information provided by the second photodetector array significantly improves the spatial resolution in the transverse coordinates (van der Laan et al 2006(van der Laan et al , 2012. In this work we aim to quantify the influence of the improvement in detector spatial resolution obtained with the DSR approach on the quality of reconstructed images acquired with a whole-body clinical TOF-PET scanner based on monolithic scintillator detectors. To this end, we perform Monte Carlo simulations of full rings of BSR and DSR monolithic scintillator detectors, modelling the detector response according to the results of characterization experiments performed on prototypes of the two types of detector at our labs (Borghi et al 2016a(Borghi et al , 2016b. Using clinically relevant metrics for the evaluation of image quality, such as the contrast recovery coefficient (CRC) and the area under the localized receiver operating characteristic curve (ALROC), we quantify the improvement offered by the DSR configuration in comparison to the more conventional BSR geometry. Finally, we simulate rings of pixelated detectors with three-layer DOI capability and a variety of pixel widths, so as to determine which decrease in pixel size produces an improvement in image quality that is equivalent to the improvement provided by the additional photosensor of the DSR monolithic scintillator detector.

Experimental characterization of monolithic scintillator detector response
The experimental characterization of the two types of monolithic scintillator detector considered in this work has been reported in detail elsewhere (Borghi et al 2016a(Borghi et al , 2016b. In brief, the detectors are based on dSiPM arrays (Philips Digital Photon Counting, DPC-3200-22), which have dimensions of 32 mm × 32 mm and consist of 8 × 8 dSiPM pixels at a pitch of 4 mm (Frach et al 2009, Schaart et al 2016. Either one or two dSiPM arrays were optically coupled to 32 mm × 32 mm × 22 mm commercially available LYSO:Ce crystals (Crystal Photonics), in BSR or DSR configuration, respectively.
The crystals were irradiated with a perpendicular incident pencil beam of 511 keV annihilation photons from a 22 Na point source, using a tungsten collimator with an aperture of 0.5 mm. The 32 mm × 32 mm surface of each crystal was irradiated at a grid of 128 × 128 reference positions with a pitch of 0.25 mm. A detailed description of the irradiation setup used for this measurement is reported in Borghi et al (2015). At each irradiation point, 100 events in the energy photo-peak were registered, for a total of approximately 1.6 million events. For each gamma event, we recorded the corresponding light distribution as measured by the dSiPM array(s), i.e. 128 pixel photon counts for the DSR and 64 for the BSR detector.
Using our k-NN classification methods (Borghi et al 2016a(Borghi et al , 2016b, the position of the gamma interaction point of each event in the dataset was estimated from the measured light distribution. This resulted in a list of 'true' irradiation positions and the corresponding 'experimental detector response'. A similar list composed of approximately 1.7 million events was produced from a side irradiation experiment, using the same irradiation setup as described in Borghi et al (2015). Here, the DOI was estimated using the clustering method presented in Borghi et al (2016aBorghi et al ( , 2016b, in which the squared sum of the pixel intensities was used as a measure to discriminate between different DOI levels. The resulting data were used to model the spatial response of the BSR and DSR monolithic scintillator detectors as explained in the following section.

Monte carlo simulations
2.2.1. Simulated system geometry. Monte Carlo simulations of the transport of annihilation photons were performed using GATE (geant4 application for tomographic emission) (Jan et al 2004). The whole-body PET system considered for this study has an inner diameter of 65 cm and an axial field-of-view (AFOV) of 19.2 cm, see figure 1. The choice of such a relatively compact geometry is motivated by the possible integration of the PET ring into an MRI system for simultaneous PET-MR imaging. The PET scanner is composed of 32 sectors, each of which contains 2 × 6 (transverse × axial direction) detectors. Each single PET detector was modelled as an LYSO block with transverse dimensions of 32 mm × 32 mm and a thickness of 22 mm, to match the size of our prototype detectors. No gaps were assumed between crystals belonging to the same sector, so that each sector has an area of 6.4 cm × 19.2 cm. Two types of detector were considered in this study: (i) monolithic scintillator detectors (DSR and BSR), and (ii) pixelated detectors with DOI capability (different sizes of the pixels were investigated).
As the primary focus of this work is on spatial resolution, the modelled detectors only differ in their spatial response, as explained in the next section. All other parameters were fixed, so as to study solely the influence of differences in the detector spatial response on the reconstructed image quality. Specifically, an energy resolution of 11% FWHM and a CRT of 200 ps FWHM were assumed for both detectors. The value of the CRT used lies in between the average values (147 ps FWHM for DSR and 214 ps FWHM for BSR) measured with the two prototype monolithic scintillator detectors (Borghi et al 2016a(Borghi et al , 2016b. The modelled energy resolution thus represents a conservative value that can be achieved with both types of detector.

Detector response.
Each gamma ray was tracked and all the interactions (photoelectric effect, Compton scattering, Rayleigh scattering, electron ionization) of all particles within the crystal were recorded as discrete interactions (so-called hits). Such a collection of hits (called event) was then used as the basis for modelling the detector spatial response. The output of the model is a single point in space, which represents the detector estimate of the position of interaction of the gamma photon.
To model the spatial response of a monolithic scintillator detector we first selected, for each event, the coordinates of the first interaction of the gamma ray within the crystal (which we interpret as the 'exact' position of interaction). We then blurred this position of first interaction by randomly sampling from the experimental detector response corresponding to a true irradiation position equal or closest to the position of the first hit. It is emphasized that the effects of Compton scattering within the detector are automatically taken into account by this procedure. Position blurring was carried out separately for the transverse coordinates and for the depth of interaction, using the data from the front-and side-irradiation experiments, respectively. It must be noted that, while the experimental detector responses were collected with a perpendicularly incident beam, the annihilation photons in a PET acquisition enter the crystal at different angles of incidences. It has been shown, however, that there is very little degradation in the detector response of monolithic detectors for angles of incidence of up to 60 degrees (Borghi et al 2016b).
To model the spatial response of pixelated detectors, we considered each 32 × 32 × 22 mm 3 block as subdivided into N x × N y × N z smaller scintillator crystals (pixels) of equal size, where N z represents the number of DOI layers and N x × N y are the numbers of crystals in the transverse dimensions. For every event, we calculated the coordinates of the energy-weighted centroid of all hits. The detector response was then defined as the center of the pixel which contains the energy centroid.
We used N z = 3 (three DOI layers), which allowed us to reduce the parallax error and obtain a wide range of performance by only varying the transverse size of the pixels. Different values of N x = N y were investigated; more specifically, we considered systems with the following values for N x = N y : 32, 24, 16, 11, 10 and 8, corresponding to pixelated arrays with pixel sizes (size = 32 mm / N x ) ranging from 1 mm to 4 mm. It must be noted that, although we do not have any specific detector design in mind, the proposed model for the pixelated detector is representative for detector blocks based on Anger logic with light-sharing read-out which make use of less-than-one photosensor per crystal element (Gross-Weege et al 2016).

Simulated source
The simulated source used in this work was a 27 cm diameter, 19 cm long, uniform cylinder containing a total of 16 hot spheres. All spheres had a diameter of 0.5 cm and were distributed at two radial positions, 5 cm and 10 cm, and two axial positions, AFOV/4 and -AFOV/4 (figure 2). The activity uptake was 7:1 for the contrast recovery study. For the lesion detectability study, the activity uptake had to be lowered to 4:1 to create a more challenging detectability task capable of pointing out differences in the performance of the investigated systems.
Since this work required a large number of simulation runs, some approximations were necessary to reduce the computational burden, as explained in the following. No ion or positron source was modelled (implying no blurring effects due to positron range); instead, each decay was assumed to generate two simultaneous 511 keV gamma rays emitted at an angle close to 180° (i.e. photon acollinearity was modelled). As the scatter fraction and absorption are the same for all simulated systems, no attenuating medium was modelled. Unless stated otherwise, the number of true coincidences acquired in a single scan was 25 million, which is at the high end of the range of current clinical 18 F-FDG imaging.

Data analysis
2.4.1. Image reconstruction. Images were reconstructed using a list mode TOF ML-EM code (Vandenberghe et al 2000). Siddon ray tracing (Siddon 1985) was used to calculate the elements of the system matrix. No resolution modelling was used for image reconstruction. Random coincidences were removed from the list mode file (random fractions are expected to be the same for all systems). Thus, only true coincidences were used for image reconstruction. The voxel size of the reconstructed images was 2 × 2 × 2 mm 3 .

Contrast recovery study.
In order to assess the ability of each investigated system to estimate small-lesion uptake, we used the CRC metric. For each sphere i in the lesion phantom, the CRC was defined as: where C H,i was measured as the mean image voxel value in a spherical region-of-interest (ROI) drawn over the hot sphere under consideration; similarly, C B,i was determined as the mean voxel value in a spherical ROI drawn at a location axially symmetric to that of the hot sphere; 'uptake' is the input lesion contrast value (equal to 7 in this work). As a figure-of-merit (FOM) for comparing different systems we report the average CRC, where the average is performed over all 16 spheres and 10 different image realizations.
In order to separately study the influence of the various physical effects that affect the CRC, we performed a number of additional simulations, starting from an 'ideal' acquisition system and progressively enabling all the main effects that degrade the spatial resolution and, therefore, the CRC. Specifically, we studied the CRC curves for the following systems: (i) the ideal system with ideal source, in which the source emits perfectly collinear gamma rays and the detectors are able to identify the exact lines-of-response (LORs) (i.e. the detectors can perfectly localize the position of the first gamma interaction); (ii) a system with ideal detectors as under (i), but with a source model that includes photon acollinearity; (iii) a system with a somewhat more realistic detector model, in which the detector spatial response is equal to the energy-weighted centroid of all the energy depositions in the crystal. Although there is no particular physical meaning associated with the energy centroid, we used it as a means to quantify the spatial blurring due to the fact that current clinical detectors cannot distinguish between multiple energy depositions in the crystal. Finally (iv), we quantified the additional loss of contrast due to the partial volume effect. Details of this calculation are reported in the appendix.

Lesion detectability.
In addition to the CRC, we also investigated lesion detectability for the lesion phantom using a mathematical observer based on a generalized scan-statistic model (Popescu andLewitt 2006, Surti andKarp 2011). Local contrast (c) at a given point in the image space was measured as the ratio of the mean counts in a circular ROI centered at that point (diameter 5 mm) and the mean counts in an annular ROI around the circular region (inner diameter 6 mm, outer diameter 20 mm). Local contrast was measured over 160 lesion copies (16 spheres in the lesion phantom × 10 image realizations).
Kernel density estimation (KDE) was used to estimate the lesion contrast probability density function f(c) from the 160 samples. Similarly, the noise nodule distribution p(c) was estimated by measuring the local contrast for each background voxel in the same image plane containing the lesions (excluding the voxels containing the lesion). Only contrast values > c 1.0 were considered for the analysis.
Finally, using the generalized scan-statistic model described in Popescu and Lewitt (2006), we were able to calculate the max-scan distribution g(c), representing the probability density function of the local contrast of the most suspicious noise nodule in a background image (i.e. image with no lesion present).
The ALROC was calculated from first principles and used as a metric for lesion detection and localization. ALROC values were calculated for each system investigated and as a function of the number of EM iterations. Since we found that the maximum ALROC was reached already after the first iteration for all systems, we will only show results for this choice of the iteration number. The error bars on the ALROC values were estimated as two standard deviations over 100 bootstrap replicates.

Experimental detector characterization
To illustrate the differences in spatial response between the BSR and DSR monolithic scintillator detectors, some of the experimental results from Borghi et al (2016aBorghi et al ( , 2016b, measured using the methods summarized in section 2.1, are reproduced in figure 3. The left-hand plot shows a cross-section along one of the transverse coordinate axes (viz., the x-axis) of the 2D probability density function of the positioning error on the entry point of perpendicular 511-keV gamma rays (i.e. the point spread function, or PSF), for the DSR (red line) and the BSR (black line) detector. It is noted that the volume under the 2D PSFs equal unity and that the cross-sections shown in figure 3 (left) have not been renormalized. Clearly, a substantially more peaked distribution is obtained with DSR, with a full-width-half maximum (FWHM) of 1.1 mm and a full-width-tenth maximum (FWTM) of 2.5 mm, compared to the BSR detector (FWHM of 1.7 mm and FWTM of 5.0 mm). Similar values were found for the y-coordinate.
The side irradiation experiment (right-hand plot in figure 3) resulted in a DOI resolution of 2.4 mm FWHM and 5.6 mm FWTM for the DSR detector, compared to 3.7 mm FWHM and 11.1 mm FWTM measured with the BSR detector. In the following we will show how the better spatial response of the DSR detector translates into improved image quality when these detectors are used in a TOF-PET system. Figure 4 shows the CRC curves for the systems based on DSR (red line) and BSR (black line) monolithic scintillator detectors. For comparison, the CRC values for systems based on pixelated detectors with three-layer DOI capability are also shown. Decreasing the pixel size from 4 mm to 1 mm results in higher CRC values, consistent with (Surti et al 2013, Thoen et al 2013. For the purpose of this work, we emphasize the results obtained for pixel sizes of 1.3 mm (blue line) and 3.2 mm (green line) with three-layer DOI classification, which we found to best match the performance of the DSR and BSR detectors, respectively. The comparison between the two types of detectors shows that, with regard to the improvement in CRC values, the impact of the additional photosensor of the DSR monolithic scintillator detector is 'analogous' to a reduction in pixel size by about 2 mm in a pixelated detector.

Contrast recovery
It must be noted that, in absolute terms, the CRC values are well below 0.5 for all systems investigated. An analysis of the contributions of the various effects responsible for the loss of contrast can be made by looking at figure 5. The purple line represents the performance of the 'ideal' system with an ideal source, defined in section 2.4.2. In this case there is no positioning error in the LORs, hence the partial volume effect is the only cause of signal loss. As reported in the appendix, the partial volume effect is expected to produce a decrease in contrast of about 0.2, given the lesion size and the voxel size used in this study. This is consistent with the CRC value of 0.8 observed for the 'ideal' system.
When photon acollinearity is enabled in the Monte Carlo simulation, significantly lower CRC values are obtained, as shown by the green line in figure 5. To put this in perspective, the spatial blurring effect due to photon acollinearity for our scanner geometry is ~1.5 mm; note that this work assumes a smaller scanner diameter (65 cm) than typical for commercial whole-body scanner designs (85-90 cm), resulting in a commensurately smaller influence of acollinearity.
Finally, a significant drop in the CRC values is observed (blue line) when the positioning error due to multiple interactions within the crystal is taken into account (by using the energy centroid for positioning). For reference, figure 5 also shows the CRC curve for the DSR monolithic scintillator detector (red line). If we look at the detector contribution only, it appears that the performance of the DSR monolithic detector is essentially limited by the physics of gamma ray interaction within the scintillator. It must be noted, however, that we are only considering a single-hit estimator in this work. In fact, leaving computational issues aside, the multiple signals measured by the photosensor arrays of the monolithic scintillator detector (light intensities and time stamps) may allow for more sophisticated models (Hunter et al 2013 that, in principle, can be included within the image reconstruction algorithm (Barrett et al 1997), so as to potentially improve the image quality.  In order to study the magnitude of the main physical factors affecting contrast recovery, we start by modelling an 'ideal', system in which all physical causes of spatial blurring have been eliminated (purple line) and progressively switch on the various physical effects that degrade the spatial resolution and, therefore, the CRC. See text for details. V Tabacchini et al Phys. Med. Biol. 62 (2017) 2018

Lesion detectability
For illustrative purposes, figure 6 shows some sample images of the lesion phantom, from a total of 10 replicates used to determine the corresponding ALROC for each investigated system. Figure 7 shows the results of our comparative study. ALROC values of 0.76 ± 0.04 and 0.52 ± 0.04 were found for the DSR and BSR monolithic scintillator detectors, respectively. Similar values of 0.69 ± 0.04 and 0.48 ± 0.04 were obtained for the corresponding 1.3 mm and 3.2 mm pixelated systems with three-layer DOI capability.
The observed difference in ALROC values for the BSR and DSR detectors is mainly due to the increased lesion local contrast of the DSR detector as shown in figure 8 where the lesion contrast (f, red dotted line) and the max-scan noise nodule contrast (g, black solid line) are plotted for (a) the BSR and (b) the DSR detector and also for the pixelated detector with (c) 3.2 mm and (d) 1.3 mm pixel size. The noise contrast distribution is roughly independent of the system since it is expected to depend mainly on the time resolution and the number of counts in a single scan, which, for the results in figure 8, were fixed parameters.  Our findings suggest that the DSR detector has the potential to achieve the same lesion detectability performance as the BSR detector with fewer counts. In order to investigate this possibility quantitatively, we determined the ALROC at different statistics (20, 15 and 10 million true coincidences). The results are shown in figure 9. It is indeed observed that the DSR detector performs similar to a BSR detector at a lower number of counts. In other words, with regard to small-lesion detectability, the improved spatial response of the DSR (leading  to higher lesion contrast) over the BSR detector can be seen as equivalent to an increase in sensitivity (of up to 20%).
As a final remark, we would like to emphasize that the focus of this work was on the spatial response of the different detectors investigated. In particular, the results presented in this work did not take into account other parameters that affect image noise, such as TOF resolution, even though experimental work has shown an improved CRT of the DSR detector compared to the BSR detector (Borghi et al 2016a(Borghi et al , 2016b. Other practical aspects such as the cost and complexity of each proposed design were also left aside.

Conclusions
This work presents Monte Carlo simulations of clinical rings of monolithic scintillator detectors that utilize experimental data for modelling the detector response. The results show that the better spatial response of a detector with DSR, compared to a crystal with BSR, translates into a significant improvement in image quality for high-resolution studies. Specifically, when imaging small lesions (i.e. 0.5 cm diameter), the system based on DSR monolithic scintillator detectors produces, on average, images with significantly higher CRC values and increased lesion detectability as measured by the ALROC. The improvement in the ALROC was observed at different levels of count statistics, ranging from 10 million to 25 million true coincidences. It was found that the DSR detector can achieve the same lesion detectability performance as the BSR detector with about ~20% less counts.
For a comparison with more conventional PET detector technology, we also simulated rings of pixelated detectors with DOI capability. It was found that, to produce the same improvement in image quality as the one provided by the additional photosensor of the DSR detector, one would need to reduce the size of the individual crystals in such a pixelated system from 3.2 mm with three DOI layers (to match the performance of the BSR detector) to 1.3 mm with three DOI layers (to match the DSR detector).
As a final remark, it must be noted that these results are specific to the imaging task and metrics considered in this study (including the object imaged and the number of true events acquired). Since we are not targeting this study to a specific application, but rather trying to draw broader conclusions, the metrics used are meant to be representative of imaging performance and closely tied to the physical performance of the instrument. Furthermore, when looking at the absolute values obtained for the performance parameters, it must be kept in mind that object scatter and attenuation were neglected, both of which affect signal-to-noise ratio and detectability measures. where V v represents the volume of the voxel. The average value of ( ) r h over the hot sphere S is: The parameter C H can be used to compute the CRC value as defined in equation (1). It must be noted that C H depends on several parameters, such as the size of the sphere, the size of the voxel, the contrast, and the position of the sphere relative to the lattice points. In our study we used C = 7, the hot spheres had a 0.5 mm diameter, and the voxel size was 2 × 2 × 2 mm 3 . We solved equation (A.5) using Monte Carlo integration, for two different cases (see figure A1): (i) the sphere is centered on a lattice point (best case) and (ii) the center of the sphere is located on the center of the voxel (worst case). We found that the CRC varies from 0.8 to 0.6. For our particular choice of the position of the spheres in the simulated source (figure 2) we expect on average a CRC close to 0.8. Figure A1. Voxelized representation of a 5-mm-diameter hot sphere in a uniform background with a 7:1 contrast. The voxel size is 2 × 2 × 2 mm 3 . Two different alignments of the sphere with the voxel grid are shown, which result in different amounts of spill-out. Although the figure illustrates this effect in 2D, it is actually a 3D effect.