Calibration of Gamma Ray Impacts in Monolithic-Based Detectors Using Voronoi Diagrams

Molecular imaging systems, such as positron emission tomography (PET), use detectors providing energy and a 3-D interaction position of a gamma ray within a scintillation block. Monolithic crystals are becoming an alternative to crystal arrays in PET. However, calibration processes are required to correct for nonuniformities, mainly produced by the truncation of the scintillation light distribution at the edges. We propose a calibration method based on the Voronoi diagrams. We have used <inline-formula> <tex-math notation="LaTeX">$50 \times 50 \times 15$ </tex-math></inline-formula> mm<sup>3</sup> LYSO blocks coupled to a <inline-formula> <tex-math notation="LaTeX">$12\times 12$ </tex-math></inline-formula> SiPMs array. We have first studied two different interpolation algorithms: 1) weighted average method (WAM) and 2) natural neighbor (NN). We have compared them with an existing calibration based on 1-D monomials. Here, the crystal was laterally black painted and a retroreflector (RR) layer added to the entrance face. The NN exhibited the best results in terms of <italic>XY</italic> impact position, depth of Interaction, and energy, allowing us to calibrate the whole scintillation volume. Later, the NN interpolation has been tested against different crystal surface treatments, allowing always to correct edge effects. Best energy resolutions were observed when using the reflective layers (12%–14%). However, better linearity was observed with the treatments using black paint. In particular, we obtained the best overall performance when lateral black paint is combined with the RR.


I. INTRODUCTION
R ADIATION detectors are extensively used in the field of nuclear and atomic physics, characterizing particles interacting with them. This requires precise determination of their deposited energy and 3-D impact coordinates. These quantities are accurately estimated employing calibration procedures addressing nonuniformity responses of the detectors [1]- [4]. In particular, γ -ray detectors are of special interest in both high energy and medical physics. They are key components of molecular imaging systems, such as gamma cameras, single-photon emission computed tomography (SPECT), or positron emission tomography (PET) scanners [5].
In the particular case of PET, providing accurate reconstructed process requires: 1) precise determination of the energy, XY planar coordinates as well as depth of interaction (DOI) of the γ -ray within the scintillation crystal; 2) timing calibration when this information is included in the reconstruction process; and 3) a correction of nonuniformities, as a result of the different detector components or manufacturing processes between different blocks [6]- [8].
Radiation detectors for PET are, in most of the cases, based on pixelated or monolithic scintillation crystals coupled to high-density photodetectors [9]. The advantages of each scintillator configuration have been extensively described in the literature (see [10] and references therein). On the one hand, in pixelated-based detectors, the estimation of the 2-D photon interaction position is basically carried out by identifying the pixel that provides the maximum signal value. Photon DOI estimation typically requires additional hardware [11], [12] or the use of algorithms, such as maximum likelihood (ML) [13]. Notice that DOI information is especially important for small aperture scanners configurations such as in small animal imaging or organ-dedicated systems [14] since it allows one to correct for the parallax error [15]. On the other hand, in monolithic-based detectors, the generated distribution of scintillation photons covers many photosensors and the position of the photodetector element with the maximum signal does not always correspond to the estimated centroid of the scintillation light distribution (LD). However, it is possible to characterize the LD profiles allowing to directly estimate the 3-D coordinates of the γ -ray interaction [16]. As a drawback, there might be scintillation light reflection from the inner faces of the scintillator, as well as truncation of the LD as a result of the finite detector size, producing a mispositioning of the γ -ray impact. This effect, known as edge effect or bias, is typically characterized by a shift in the impact position determination toward the crystal center that becomes stronger at the edges of the block.
The most accurate approach for the estimation of the impact position in monolithic blocks would be to readout every photosensor element. Statistical methods, such as ML algorithms [17] or k-nearest neighbor (k-NN) methods [18], [19], This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ have been proposed. These methods make it possible to determine the photon impact position by comparison of the LD shape at different interaction positions with a set of calibration data stored in look-up tables. Recently, it has been demonstrated the possibility to apply artificial neural network (ANN) [20] or gradient tree boosting (GTB) algorithm [21], [22] for this purpose. However, those methods typically require large acquisition times, and hardware with accurate positioning capabilities. As an alternative, the analytical methods that model the relation between the source position and the measured photodetector pixel signals, using for example weighted least-squares (WLSs), can be employed without prior calibration dataset [23], [24]. These methods could also be applied when using reduction readout schemes which reduce the complexity and cost when compared to reading out every single-photosensor element [25].
Electronic configurations based on the networks of passive components have been proposed to reduce the number of readout channels. The most traditional is the so-called anger logic [26], which returns only four output signals allowing to implement the center of gravity (CoG) algorithm to estimate the XY interaction position. A modification of the resistivenetwork, providing the collected charge (scintillation photons) for all rows and columns of the photosensor arrays, permits to better measure the centroid of the LD and makes it possible to additionally estimate the photon DOI [27]. However, the edge effect is still present in the peripheral region of the detector when statistical and analytical methods are used, especially when some type of reflection treatment is applied to the crystal faces. Modifications of the CoG, such as the so-called raise-to-power (RTP) algorithm [28], allows one to partially mitigate that effect, but the need of calibration procedures is still required.
The edge effect of the estimated XY positions (anger logic, CoG, RTP, etc.), can be further reduced by means of bilinear [29], 2-D polynomial [30], or 1-D polynomial [31] interpolation methods using calibration masks. These approaches have been applied in several systems based on monolithic crystals [6], [7], [31]. However, while the interpolation between mask points is well achieved, extrapolation beyond the outmost point sources (mask holes) may lead to some artifacts and, typically, events at the edges are rejected reducing system sensitivity.
Improvements on the calibration procedure for monolithicbased PET detectors would enhance their performance. This article focuses on the development of a calibration method to accurately determine 3-D position and energy of γ -ray impacts for the whole monolithic crystal volume and, thus, also increasing system sensitivity. The proposed method is based on the mathematical structures called Voronoi diagrams [32]. Those diagrams are currently being used for automatic target volume definition during treatment planning in radiotherapy [33], image correction of deformable motion in PET image reconstruction [34] or, in pixelated-based PET detectors, for energy calibration procedures [35] and crystal identification in a multiple-layer configurations [36]. Moreover, Voronoi diagrams are also used  in sports [37], [38], chemistry [39], [40], astronomy [41], medicine [42], and image processing [43].

A. Detector Block
Two thick LYSO monolithic scintillators with dimensions of 50×50×15 mm 3 have been used (see Fig. 1). All crystal faces are polished and one of the 50 × 50 mm 2 surfaces (exit face) coupled to a photosensor array by means of optical grease (BC-630, Saint Gobain) to reduce light transmission losses. We have tested six different treatments on the same crystals, as depicted in Fig. 1 and described in Table I For each study, the same crystal and SiPM array was used. Each SiPM array (ARRAYC-30035-144P-PCB, SensL/OnSemi) was composed by 12 × 12 photosensors and covering an approximate area of 50 × 50 mm 2 . Each individual photosensor has an active area of 3 × 3 mm 2 and the pitch is 4.2 mm in both directions. Each array is typically operated at a bias voltage of 31 V, 6.5 V over the breakdown voltage. Since SiPMs are sensitive to temperature variations, the detectors are kept at stable temperatures of approximately 15 • C-18 • C, also reducing dark noise contributions.
The detector readout provides information for each row and column of the SiPM array allowing to characterize the scintillation LD profiles. This readout provides a reduction from 144 SiPMs output signals to only 24. The 24 signals are digitized using custom ADC boards (12-bit precision), using an integration window of about 250 ns. All acquisitions were performed in coincidence mode using an identical reference detector, and a coincidence window of 5 ns. The digitized and synchronized signals are sent to a workstation where they are processed (see description below). Fig. 2 shows a schematic of the steps followed during the data calibration. First, using the signals provided by the analog readout, the 3-D photon impact coordinates, and energy of each coincidence event are estimated applying the RTP algorithm. XY coordinates are calculated by raising the 12 digitized signals for each projection to the power of two [44], [45], before CoG calculation. The DOI coordinate, Z, is estimated for each event as the average for rows and columns (r, c) of the ratio of the sum of all 12 signals (photon energy, E) to its maximum value (E/I max ) r,c [46]. This process is labeled as RTP2 in Fig. 2(a).

B. Calibration Method Based on Voronoi Diagrams
1) Detector Calibration Procedure: Reference data sets, named calibration maps, have been acquired for each crystal treatment using an array of 11 × 11 22 Na sources that covers an area of 46 × 46 mm 2 . The first row of sources is located at 2 mm from the crystal border and the pitch is 4.6 mm. The sources were mechanically collimated using a tungsten mask of 24-mm thickness and 1.2-mm drilled holes [see Fig. 2 An array with a smaller pitch would have implied overlapping of the sources at the edges of the flood maps. During data processing, each detector area is binned in 600×600 pixels. Once the calibration map has been acquired, the detector calibration is done as follows.
1) 3-D photon impact coordinates (and energy) result in flood maps as the one shown in Fig. 3(a). All 121 calibration sources were correctly identified. A software collimation in the range of 1 • (depending on this article) was applied between both detectors helping to better determine the distribution centroids. That means that only events whose line of response (LOR) is contained within this specific angle are considered [47].  Fig. 3(a). The detector calibration method is based on Voronoi diagrams. These diagrams, also known as Dirichlet tessellation, are defined as the partitioning of the plane into various convex polygons T i , named Voronoi cells, each of them containing one generating point, named Voronoi point. An arbitrary point lies within a specified Voronoi cell if, and only if, the distance from this point to the Voronoi point of its associated polygon is smaller than all other distances between this point and the remaining Voronoi points [32]. This can be mathematically expressed as where is the Voronoi diagram generated by X, with X the set of the m Voronoi points, and d(x, xi) is the Euclidean distance on R 2 . In our case, since we are using a calibration mask of 11 × 11 sources (m = 121), the crystal surface is divided into 11 × 11 Voronoi cells [see Fig. 3(b)]. The yellow line delimits the "intuitive interior" of the 121 Voronoi points, named convex hull. In this article, the Voronoi diagrams have been computed using a specific MATLAB function (voronoi) [48].
Once the Voronoi cells are determined, five Voronoi factors for each Voronoi cell, one for X, one for Y, one for energy and two for Z, are independently calculated. For the XY and energy Here, (x, y, E) known source and (x, y, E) floodmap source correspond to known and measured XY coordinates and energy of each calibration source. Regarding the energy correction, a Gaussian fit to the energy spectra of each Voronoi cell was applied providing the photopeak value in ADC units, (E) floodmap source . We have used the central cell, (E) known source , to convert the ADC units to 511 keV. For the DOI calibration, we have obtained the limits of the E/I max histograms of each Voronoi cell, named a and b, by using the analytical expression for the DOI distribution [27] To calibrate the measured E/I max units to mm we have considered two Voronoi factors for each Voronoi cell (VoronoiFactor (Z1,Z2) ) that correspond to the two parameters of a linear fit considering the limits a − σ int and b + σ int equal to 0 and 15 mm (crystal thickness), respectively.
2) Interpolation Methods: Once the Voronoi factors are calculated, each recorded event (x, y) is added to its corresponding Voronoi cell. However, they are calibrated by the interpolation of the Voronoi factors of the closest cells, as depicted in Figs. 4 and 5. Thus, the chosen interpolation method plays an important role. They cannot be indiscriminately used and, therefore, it is important to understand their principles and limitations. Moreover, the events near the edge should not be discarded, so the selected method must allow one to also extrapolate the events that are beyond the convex hull, represented again in Fig. 4 by the yellow line.
Two different interpolation methods have been studied.  Weighted average method (WAM) of the Voronoi factors. Each calibration factor of one recorded event (f j ) corresponding to each position coordinate and to the energy, is calculated as follows: y) is the distance between the recorded data point and the Voronoi points and k is a chosen exponent. For each recorded event inside the convex hull, three closest Voronoi points have been selected (n = 3; at left, right, up, or down) [see Fig. 5(a)]. The calibration factors of the recorded events that are located outside the convex hull are determined directly from the Voronoi factors of their associated Voronoi cell, without any interpolation from the closest cells. We have tested different k values (not shown here), suggesting the best results for the squared power (k = 2) and therefore, for simplicity, we are only providing these results in this article.
Natural neighbor (NN) method of the Voronoi factors. In this article, this interpolation has been made with a specific MATLAB function (scatteredInterpolant) [48], [49]. The recorded event (x, y) is added as a new Voronoi point and, therefore, generates a new cell that intersects with certain Voronoi cells of the former set of 121 ones, named NNs [see Fig. 5(b)]. The interpolation is a weighted average area of the NNs of the recorded data, where f j is calculated as follows: where y) is the Voronoi cell of the recorded event and T i (x, y) is the NNs cells. Moreover, the calibration factors of the events located outside the convex hull are calculated using an extrapolation method based on a least squares approximation of the gradient at the boundary of the convex hull [48].
In the WAM, weights are distance-based and, therefore, the weights assigned to each Voronoi cell diminishes as the distance from the recorded data to the Voronoi point increases.
In the NN interpolation, weights are area-based and, therefore, depend on the area of the NN cells which are inside of the new Voronoi cell. Here, larger areas result in the larger influence of the corresponding Voronoi factor on the interpolation value.
The calibrated data for each recorded are calculated multiplying the estimated 3-D coordinates and energy event (x med , y med , z med , and E med ) by its appropriate calibration factor (f j ), as follows:

C. Experimental Validation
To validate the calibration data process carried out using the 11 × 11 22 Na sources array (4.6-mm pitch), a second set of data has been acquired using another array but with 9 × 9 22 Na sources (5-mm pitch). The studied cases are as follows.
1) Comparison of the Voronoi interpolation methods and 1-D monomial approach. The 1-D monomial interpolation [31] also makes use of the 11 × 11 sources array to perform the calibration. We have used the crystals with RR treatment since they have been successfully tested before [44]. The 9 × 9 22 Na sources array was placed in front of the crystal under evaluation. The position of these sources is different from those used during the calibration (11 × 11). Having sources outside the convex hull (less than 2 mm to the crystal edge) would have been challenging to resolve. The reference detector was located at 416 mm and just software collimation of 0.6 • (total aperture) was applied. 2) Crystal Surface Treatments: The aim here was to test the calibration using the NN interpolation method for  the six different crystal treatments because as it will be shown in the results section, it exhibited the best performance in the case 1). In order to reduce possible random coincidences, the 9 × 9 22 Na sources were additionally collimated using a tungsten mask with 1.2 mm in diameter drilled holes in a 24-mm thick block. The collimator was attached to the entrance face of the detector under calibration. The reference detector was placed at 110 mm and a software collimation of 1.1 • (total aperture) was applied after the calibration. Concerning the XY calibration, we have evaluated the central row and column of the 9 × 9 22 Na sources array [see horizontal and vertical yellow bands in Fig. 6(a) left], before and after calibration, for all cases. The known mechanical source position defines the true position of the sources. The centroid of each source in the flood maps is calculated using a multi-Gaussian fit. Several parameters have been used to evaluate the performance of the XY calibration. 1) Bias: Difference between measured and true source position. 2) Linearity: Mean distance between sources. We have also calculated the confidence intervals (CIs). We provide the 65% CI, this is a range where one can be 65% certain it contains the mean value. Finally, the energy calibration has been evaluated with the sources located along one diagonal in the flood maps (see also diagonal yellow band in Fig. 6) both before and after calibration. The energy resolution is calculated as E(FWHM)/E centroid using a Gaussian plus a linear distribution.

A. Comparison of Voronoi Interpolation and 1-D Monomial
The noncalibrated flood map of the 9×9 22 Na sources array acquired with the RR crystal is shown in Fig. 6(a) left. A slight shift of the sources along the Y-axis is observed due to a misalignment of the array positioning system. We have calibrated the 3-D impact positions and energy using the interpolation methods described above (named WAM and NN), and the 1-D monomial approach, also shown in Fig. 6(a). The profiles of the central row of sources, together with a multi-Gaussian fit (red line), are shown in Fig. 6(b). Fig. 7 shows the bias, distances between sources, spatial resolution (FWHM), and the R bias coefficient, for the four studied cases. The panels on the left correspond to the central row (X-axis) of the flood maps and the panels on the right to the central column (Y-axis). The average standard deviation (X-and Y-axes) of the bias is 1.15 mm for the noncalibrated data. This value is reduced to 0.34, 0.33, and 0.35 mm for the calibrated data using the WAM, NN, and 1-D approaches, respectively.
The linearity and CIs are listed in Table II for both X-and Y-axes. More accurate values are obtained with the Voronoi approaches, closer to the actual value of 5 mm (see also Fig. 7). The average linearity for the X-and Y-axes is 5.3±0.2, 5.1±0.2, 5.1±0.1, and 5.2±0.2 mm, for the noncalibrated data and WAM, NN, and 1-D cases, respectively.
The spatial resolution FWHM as a function of the mechanical source position is also depicted in Fig. 7. FWHM values of the noncalibrated sources closest to the crystal edge diminish for both x-and y-axes due to the image compression. FWHM values worsen for the sources closest to the edge in the x-axis when the WAM is applied. However, this effect is reduced using the NN or 1-D approaches. The average FWHM for the  The DOI profiles for the noncalibrated (E/I units) and calibrated data (mm units), for the three studied ROI, are shown in Fig. 8. The DOI histograms for the noncalibrated data showed different shapes depending on the analyzed detector area [see Fig. 8(a) top panel]. After DOI calibration, the histograms resembled the expected one. That means a larger number of events at the crystal entrance (15 mm side) and lower closer to the photosensor side (0 mm). This happens for the three ROI and calibration approaches. Fig. 8(b) shows the determined upper limits, named b, for the three ROI. Fig. 9 depicts the photopeak centroid position, normalized to the central one, for the nine sources across the diagonal. The noncalibrated data exhibited lower photopeak values of about 13±1% at the crystal edges. After energy calibration, all interpolation methods returned differences as small as 5±1% at the edges, and below 2±1% in the central region. The energy resolution for the whole scintillation volume resulted in 17.4±0.2%, 13.3±0.1%, 13.0±0.1%, and 13.2±0.1%, for the noncalibrated, WAM, NN, and 1-D approaches, respectively.

B. Crystal Surface Treatments
We have obtained noncalibrated data of the 9 × 9 22 Na sources array for the six different crystal treatments. The NN interpolation method was used to calibrate the 3-D impact positions and energy for all crystal treatments. Fig. 10 left shows the flood maps of the noncalibrated data and on the right side the calibrated flood maps. The profiles for the central column (y-axis) of each flood map are also shown at the sides.    Fig. 12(a) shows the average standard deviation for the x-and y-axes of the bias for each treatment before and after calibration. The largest value (1.5±0.2 mm) obtained for the White treatment is reduced to 0.3±0.1 mm after calibration, as also found for other treatments. Fig. 12(b) depicts the average linearity. Some excess of the linearity in the noncalibrated data is observed for the treatments including lateral walls black painted. This is the opposite for both specular (ESR) or diffused reflection (White) treatments. The average linearity for all calibrated cases agrees well with the actual value of 5 mm.
The average spatial resolution FWHM is plotted in Fig. 12(c). Regarding the calibrated spatial resolution FWHM, the best value is obtained when using the RR treatment (about 1.6±0.1 mm on average). This value is slightly better than that obtained in Section III-A due to the additional mechanical collimation (see Section IV for more details). The worst spatial resolution is observed when reflective materials are used (ESR and White) approaching 2.2±0.1 mm. The average R bias coefficient for each treatment is shown in Fig. 12(d). The calibrated data exhibit results nearing zero. In particular, the RR treatment results in 0.2±0.1 mm 2 , whereas the ESR treatment shows the highest value of 0.9±0.1 mm 2 .
The energy resolution for the whole scintillation volume before and after calibration is plotted in Fig. 13(a). In contrast to the spatial resolution, ESR and White treatments exhibited the best energy resolution (12%-14%). A deterioration to 16% and 17% is observed when the lateral black paint is used in combination with other reflective materials, and to about 22% if the entire block is black, most likely due to a poorer collection of scintillation photons. Fig. 13(b) shows the photopeak variation for each treatment as a function of the γ -ray impact position for a diagonal of sources across the detector surface. The largest variation of the photopeak positions is observed at the detector edges reaching 24%, 22%, 7%, 11%, 17%, and 29% (±1%) for the Black, RR, ESR, White, B+ESR, and B+W treatments, respectively. After calibration these values are significantly reduced to only: 5%, 2%, 1%, 2%, and 2% (±1%), respectively. Fig. 14 depicts the DOI histograms corresponding to each treatment, before and after calibration. The noncalibrated data is plotted in E/I units and all other data using the NN interpolation in millimeters. The treatments with reflective material show a shift of the DOI histograms to higher values, due to the increase of collected energy affecting the estimator (E/I). After calibration, there is a good agreement in between all distributions. Indeed, Fig. 15 shows the determined lower and upper limits, a and b, for the three ROI. The upper limits exhibit a position dependence. After calibration, most of the treatments showed almost no variation of these parameters independently of the studied region. This indicates the possibility to also calibrated the photon DOI with this methodology.
The larger variation in the noncalibrated data was observed for the ESR treatment, which is nevertheless reduced after data calibration.

IV. DISCUSSION
In this article, we have introduced a calibration method based on the Voronoi diagrams to accurately determine the energy and 3-D impact positions in γ -ray detectors based on monolithic scintillation crystals.
In a first set of experiments, we have evaluated the WAM, NN, and 1-D calibration methods using a 9 × 9 22 Na sources array acquired with the RR crystal (lateral black painted and an RR layer at the entrance). On the one hand, the noncalibrated image shows the inward compression at the edges of the crystal due to the truncation of the LD. This results in an overestimation of the spatial resolution FWHM values of the sources closest to the crystal edge, as expected. Interestingly, the linearity of this noncalibrated data did not show a large mispositioning of the sources at the edges. We assumed this happened because the sources at the edges are at 5 mm from the crystal edge, where this effect is less pronounced given the photosensor density (12 × 12 SiPMs), surface treatment, readout granularity (12 + 12 signals), crystal aspect ratio (15-mm thickness and 50-mm size) and RTP2 calculation of the estimated impact position. On the other hand, when the 1-D approach is applied, we observed that it is hard to resolve data that occurs beyond the outermost calibration sources. In particular, the upper row of the 1-D approach is almost vanished, due to the slight shift of the 9 × 9 array with respect to the detector center (see Fig. 6).
The interpolation methods based on the Voronoi diagrams made it possible to calibrate the whole scintillation volume, thus increasing the detector sensitivity. When the WAM was applied, the FWHM values of the sources closest to the X edge worsened. This was caused because these sources are located between two Voronoi cells (see Fig. 4) and, thus, this method has some position challenges, producing an elongation of the sources. On the y-axis, the sources closest to the edge are not located between two Voronoi cells (see Fig. 4) due to the aforementioned 9 × 9 sources array shift and, therefore, this elongation was not observed. However, when the NN interpolation method was applied, the elongation was not shown in any axis; providing a correct interpolation in both directions.
The different shapes of the noncalibrated DOI distributions for the three ROIs are produced by a stronger scintillation light truncation toward the crystal edges and corners. After calibration, DOI distributions for all regions exhibited a similar behavior. The use of an RR layer, when compared to a totally black painted crystal, is characterized by an excess of the slope of the exponential attenuation curve [50]. Regarding the photopeak energy dependency with the impact position, a larger photopeak variation at the crystal corners was found in the noncalibrated data, especially for configurations using black painted walls, as expected. This is because in monolithic crystals with black painted walls, there is a scintillation light collection dependence with the γ -ray impact position due to the scintillation light truncation at the edges and, thus, certain  losses of light transferred and collected at the photosensors. These effects were compensated after calibration. Summarizing, the accuracy on the impact determination (FWHM, bias, and distance) using the NN interpolation performs better than the other tested approaches, especially recovering impacts at the edges of the crystal without deteriorating the FWHM, as it occurs for the WAM interpolation. Independently of the used interpolation method, an improvement on the energy and DOI performance is always observed. It is worth mentioning that the calibration of data described in this article could be implemented in reconfigurable devices such as field programmable gate arrays (FPGAs), by means of lookup tables for each Voronoi factor and detector unit.
After the validation of the Voronoi interpolation methods, in particular when using the NN approach, an additional study was carried out for different crystal treatments. Notice that in this article we used an additional mechanical collimation as compared with the previous study in order to better resolve the sources at the edges of the crystals. In the noncalibrated images is discernible how the compression effect increases when the White treatment is used. Thus, larger bias values of these data are obtained at the edges. However, for all treatments, this bias is significantly reduced when the data is calibrated using the NN interpolation. Moreover, the linearity for all the treatments agrees well with the actual value. We found the best FWHM for the RR treatment. A deterioration of the calibrated spatial resolution FWHM for the ESR and White treatments is observed, most likely due to the fact that the LD is no longer preserved, and the stronger edge effect. In contrast to the spatial resolution behavior, ESR and White treatments exhibited the best energy resolution (12%-14%) due to the larger collection of optical photons. This light collection increase also causes a shift of the measured DOI histograms to higher values. However, after DOI calibration we observed a general good qualitative agreement of all profiles with the expected gamma ray attenuation distribution. For all cases we made it also possible to calibrate the a and b parameters for the entire block.

V. CONCLUSION
We have described and validated a calibration method for monolithic crystals of large dimensions based on analytic calculations. We have combined this with high-density photosensors arrays and readout electronics. Such readout makes use of a reduction scheme of signals resulting in the number of rows plus columns of the photosensor arrays. This has been shown to be a good sampling to return accurate energy and 3-D impact coordinates of the γ -rays in monolithic blocks. The aim of this article has been to calibrate these measured energy and impact coordinates to the expected values. This has been done using the so-called Voronoi diagrams and interpolation methods. We have studied a few interpolation methods, but in this article, we have shown the two that achieved better performance, named WAM (k = 2) and NN.
We have first evaluated the two interpolation methods and compared them with an additional approach based on 1-D monomial that we have successfully used in former detector designs and systems. The tested interpolations are used not only for the XY impact position but also for the energy and now for the photon DOI. If we compared to the 1-D method, both procedures based on the Voronoi diagrams allowed one to calibrate the whole scintillation volume, without rejecting events at the detector edges, increasing the system sensitivity. Nevertheless, we observed some better calibration performance when using the NN method.
Finally, we made use of the NN interpolation method to evaluate the response of the described methodology to six different crystal treatments. Although black paint is widely used to preserve the LD expecting more accurate 3-D impact determination, other crystal treatments that enhance the light extraction using diffuse or specular reflectors are also possible. We have demonstrated that it is possible to use the NN interpolation with a variety of crystal treatments. Therefore, the selection of the crystal treatment would be a tradeoff of the aimed system geometry and application.