Among valid events acquired within 511 keV energy window (e.g. 435–585 keV), ICS events can be distinguished from PE by using energy and position information of multiple interactions obtained from photosensors. It is important to note that the energy window range may vary depending on the energy resolution of the detector elements in the PET system. Among multiple interaction information, determining the first interaction position of an ICS event is the most important part. Various methodologies have been suggested for recovering ICS events into the first interaction position. In this section, different methods for ICS events recovery will be discussed. Additionally, we assume that object-scattered event is excluded by the energy window. Although distinguishing random events within the coincidence window is challenging, they can potentially be rejected by considering factors such as the distance between detector blocks.

4.1 Energy comparison and Compton kinematics4.1.1 Energy comparison

The simplest way to estimate the first interaction position of ICS events is to compare the energy absorbed by each crystal. Based on the acquired energy signal, the first interaction can be assigned to the crystal element where (1) maximum energy deposited, (2) second maximum energy deposited, and 3) CoG position from weighted-energy signals as shown in Fig. 6.

Fig. 6

ICS recovery scheme based on energy comparison. Here, we assume a case when an incident gamma-ray with energy E0 undergoes two consecutive Compton scatterings at crystals C2 and C4, and finally absorbed at crystal C5, and energy is deposited at three crystal positions with energies E1, E2, and E3, respectively. C = (x, y) mm represents the position of each crystal center in the x and y directions

Comanor et al. [27] compared methods for recovering the first interaction position of ICS events based on energy information. They investigated and compared the maximum energy signal, second maximum signal, and joint algorithm for a 60-cm-diameter PET system with 2.5 cm axial field-of-view (FOV) consisted of 3 × 3 × 30 mm3 BGO crystals 1-to-1 coupled to 8 × 8 photosensor arrays. The joint algorithm designates the crystal as the first interaction based on the energy difference between two crystals where ICS occurred. When energy difference is larger than a specific threshold maximum signal is designated, and when energy difference is smaller than a specific threshold second maximum signal is designated. As a result, the second maximum signal and joint algorithm showed misidentification ratio of 12% which was better than that of maximum signal algorithm with 22% misidentification ratio.

Shao et al. [7] compared the accuracy of the first interaction positioning methods using maximum energy, maximum Z designating the first interaction position with the highest z-coordinated among the interaction positions, minimum Z, and weighted-energy scheme (Fig. 6). The maximum Z scheme exhibited highest ICS positioning accuracy for both 2 × 2 × 10 mm3 BGO and LSO, especially effective when gamma ray incident angle varied from 0 to 30 degrees. The minimum Z and maximum energy showed similar results, while the weighted energy scheme showed the worst performance.

Surti and Karp [28] investigated the positioning accuracy and analyzed the reconstructed image quality of using maximum energy deposition and second maximum energy, and compared them with weighted energy centroid positioning. The 1-to-1 detector coupling configuration was used while varying LSO crystal thickness from 1 to 3 cm and the crystal cross-section from 1 × 1 mm2 to 4 × 4 mm2. The positioning error, the distance from the true first interaction position to the recovered interaction position, was highest for the maximum energy in all cases, followed by weighted energy and second maximum energy. The second maximum energy deposition scheme demonstrated improved detector intrinsic spatial resolution and showed lower mean positioning error for crystal thicknesses of 1, 2, and 3 cm with all cross-section ranges. The system simulation was also conducted with a detector consisted of 4 × 4 × 20 mm3 LSO crystal and a system with an 85 cm ring diameter and a 70 cm axial FOV. The second maximum energy deposition scheme demonstrated improved system spatial resolution in both axial and transaxial directions and the contrast recovery coefficient values were enhanced by 30–40% compared to the maximum energy deposition scheme. (Fig. 7).

Fig. 7

Comparison of energy comparison methods. (a) Reconstructed spatial resolution (equivalent width: total count in profile/peak count) of a point source in air plotted as a function of radial position for a 70 cm long whole-body scanner using 4 × 4 × 20 mm.3 LSO crystals. (b) Reconstructed images of the central transverse slices for the simulated lesion phantoms with 1 cm (top row) and 0.5 cm (bottom row) diameter lesions. Moving left to right in each row are images from Anger, 1st max crystal, and 2nd max crystal positioning algorithms. The water-filled cylinder was 35 cm in diameter and 70 cm long, and sphere uptake was 3:1 with respect to background. (From [28]; with permission)

Energy comparison approach provides simple way to recover ICS events to the first interaction position, however, this approach requires 1-to-1 coupling of scintillation crystal to photosenstor pixel to obtain individual energy deposition at each crystal within the detector array. Overall, when comparing energy information for determining the first interaction position, the second-highest energy deposition method seems to offer better ICS recovery. This implies that a significant portion of Compton scattering occurs as forward scattering within the PET system, resulting in lower energy deposition at the initial interaction position (which also matches with Fig. 3b). However, it is not always accurate to assume that the second maximum energy deposition corresponds to the initial interaction position; the maximum energy deposited position can be the first interaction position when backscattering occurs. Hence, incorporating Compton kinematics in energy and position information of each crystal could help to improve positioning accuracy.

4.1.2 Compton kinematics-based method

Rafecas et al. [29] presented six recovery methods for ICS events in the MADPET-II, an LSO-APD-PET prototype PET system. They introduced Compton kinematics and Klein-Nishina equations to identify the crystal position where first interaction occurred. Six different methods were compared. (1) The averaging over LoRs (Av) scheme involves assigning a weight of 1/2 to each possible LoR and summing them. (2) The maximum energy scheme, as described in the previous section, draws LoRs based on the crystal with the highest signal. (3) The Compton kinematics scheme estimates the correct path by utilizing energy information and interaction position information through the Compton interaction formula (Eq. 1) to determine the scattering angle. (4) The modified Compton kinematics scheme considers errors in calculating the scattering angle. (5) The Klein-Nishina scheme assigns the trajectory with a higher probability between two trajectories using the Klein-Nishina formula. (6) The hybrid method is another modified Compton kinematics scheme, where the maximum energy scheme is used instead of Compton kinematics scheme under certain circumstances. Simulation study was conducted and reconstructed image analysis was performed.

The Compton kinematics scheme estimates the scattering trajectory by utilizing energy and interaction position information through Compton interaction formula (Eq. 1). Here, several assumptions are made: E0 equals 511 keV, the incident gamma-ray scatters only once and deposits the scattered radiation energy E2 in the second crystal, and the scattered radiation does not interact with non-sensitive detector materials, implying E0 = E1 + E2 = 511 keV. We can calculate the maximum energy transfer as E1 = 340.7 keV from Eq. 1, thus providing constraints for E1 and E2:

$$<E}_\le 340.7\text170.3\text\le _<511\text$$

(3)

Also, the scattering angle (θ) can be calculated in two different ways as described below, where \(\overrightarrow_}\) and \(\overrightarrow_}\) are the distance vectors corresponding to the primary and secondary trajectories, as indicated in Fig. 6. Ideally, the scattering angle calculated using measured energy information according to Eq. 4 and that calculated geometrically using the trajectories (Eq. 5) should yield the same value, but in real situations, they differ due to energy and position uncertainties.

$$\mathit_=2- \frac_}_-_}$$

(4)

$$\mathit_= \frac_}\bullet \overrightarrow_}}_}\right|\left|\overrightarrow_}\right|}$$

(5)

The following explains how the first interaction position is finally determined in the Compton kinematics scheme, based on the measured energy signals from the photosensor, E1 and E2.

(1)

If E1 < 170.3 keV, E2 will be larger than 340.7 keV. In this case, the position with E1 is identified as the first interaction position.

(2)

If 170.3 keV < E1 and E2 < 340.7 keV, choose the one with the smaller \(_\) value between the two possible scenarios (E1 or E2 corresponds to the first interaction). The \(_\) represents the absolute difference between Eqs. 4 and 5, as calculated in the following equation:

$$_=\left|cos_-cos_\right|$$

(6)

The modified Compton kinematics scheme uses a modified \(_\) value that incorporates uncertainties of energy and position information (\(_^\)). Detailed error calculations can be found in this paper.

$$_=\frac_-cos_\right|}_^}$$

(7)

The Klein–Nishina scheme is based on Klein–Nishina formula described in Eq. 2. For two possible cases, probability of Compton scattering is calculated based on Klein-Nishina formula and LoR is chosen for the one that has higher probability. This approach allows for two possible implementations, depending on the method used to calculate \(cos\theta\): either by energy considerations (KN–E, Eq. 4) or by geometrical analysis (KN–G, Eq. 5).

Finally, the hybrid scheme is a modified Compton kinematics scheme by applying the maximum energy selection criterion under certain circumstances. For triple coincidences events, this method classifies two possible LoRs of ICS event into module differences (Mij), axial slice differences (Sij), and layer sums (Lij) (Fig. 8a). Depending on the (Mij, Sij, Lij) combinations, the identification fraction was pre-calculated using the Compton kinematics. For the combinations where the identification fraction was lower, maximum energy scheme was used instead of Compton kinematics.

Fig. 8

Comparison of ICS recovery methods for a dual layer, high resolution PET scanner. (a) Axial and frontal view of MADPET-II. The possible LoRs classified in (M, S, L) combinations with module differences (Mij), slice differences (Sij), and layer sums (Lij). Right figure shows the labelling of the channels in MADPET-II system. (b) Identification fraction (IF) for Maximum Energy (ME, ▼), Compton kinematics (Co, ●), and Hybrid scheme (HC, □) as a function of module difference M. For each case, the columns represent different values of L, and the rows represent values of S

The correct identification fraction of the first interaction position was highest with the value of 59.3% for the hybrid scheme, followed by 56.6% for Compton kinematics, 56.4% for modified Compton kinematics, 50.6% for maximum energy, and 49.9% for Klein-Nishina scheme as shown in Fig. 8b. The SNR of the reconstructed hot-rod PET images with the hybrid scheme was better than with other methods, except for the Av scheme. Contrast was also highest with the hybrid method.

Abbaszadeh et al. [26] also uses Compton kinematics to recover intra- and inter-crystal scatter in a CZT-based PET system. In 'cross-strip' CZT detectors, positioning is determined at the intersection of the responding anode strip and cathode strip. Consequently, in cases of intra- and inter-crystal scatter, more than one anode and cathode may respond, leading to the occurrence of four or more intersections. In such cases, energy-based methods are not applicable due to detector characteristics. Therefore, Compton kinematics are utilized to determine possible angles of trajectories. When applying Compton kinematics-based method with a low energy threshold of 10 keV, sensitivity improvements were observed with 1.43 times higher value and contrast-to-noise ratio (CNR) improvement from 5.81 to 12.53.

4.2 Statistical approaches

Another approach to determining the first interaction position in ICS events involves employing optimization algorithms based on statistical or mathematical principles. Given that ICS events follows Compton kinematics and Klein–Nishina cross-section, Bayesian and Maximum Likelihood (ML) methods leverage probabilistic inference and likelihood maximization to estimate the first interaction position. These approaches use probabilistic models, thereby refining accuracy by integrating prior knowledge and updating probabilities that is based on Compton kinematics and Klein-Nishina formula.

When a 511 keV gamma-ray interacts at different positions after Compton scattering \(}=[_,_,\dots ]\) with energy \(}=[_,_,\dots ]\), observed detector response signal is represented by \(}=[_,\dots , _]\). We assume that each detector response follows Gaussian distribution with mean (\(\mu (},})\)) and the standard deviation (\(\sigma\)). When ICS event occurs at crystal position \(}=[_,_]\) with energy deposition \(}=[_,_]\), the probability that we observe detector response \(}=\left[_, \dots _\right]\) from a 8 × 8 photosensor array can be described as \(P\left(},}|}\right)\) and the likelihood function is given in \(\mathcal(},})\). For an estimate for maximizing \(\mathcal(},})\) is given as Eq. 10 is known as ML estimate. Moreover, by incorporating a prior information or probability (e.g. in this case prior information will be Compton kinematics or Klein-Nishina cross-section), we can estimate parameters in a Bayesian framework. This is called Maximum a Posteriori (MAP) estimation as given in Eq. 11. Weighting parameter \(\beta\) can be used to adjust relative importance of each component.

$$}\left(E,x|S\right)=\prod_}=1}^}}\frac}\sqrt}}}\mathbf[-\frac}}^}}}_}}-}\left(E,x\right)\right)}^]$$

(8)

$$\mathcal\left(E,x|S\right)=\left(}\left(E,x|S\right)\right)$$

(9)

$$}}}\left(E,x\right)=\underset},}}}\mathcal(E,x|S)$$

(10)

$$}}}\left(E,x\right)=\underset},}}}\mathcal^}}}}_}}}}}}^}}$$

(11)

Pratx and Levin [30] presented a Bayesian method for reconstructing sequences of interactions in detectors to accurately identify the first interaction in multiple interaction events in their cross-strip CZT system. The MAP estimation based on Bayesian framework was used to optimize the sequence consistency while considering the trajectory total cross-section. The study compared the performance of MAP with other positioning methods–energy-weighted positioning and minimum distance positioning using simulations of CZT PET systems. MAP showed the first interaction positioning accuracy of 85.2% and significantly reduced the mispositioning of events compared to minimum distance method, leading to improved image quality, higher CNR, and better spatial resolution and thereby making it a promising technique for high-resolution PET imaging.

Lage et al. [31] presented an inter-detector scattering recovery method for recovering triple coincidences in PET without the need of additional energy resolution requirements. The authors distributed triple events (T) among their possible LoRs based on the relative proportions of double coincidences (D1-2 and D1-3) in those LoRs, representing a ML solution (Fig. 9). LoR information is presented in histogram format and described in count values; D and T components are also described as counts in LoR histogram format. This study adapted a real preclinical Argus PET/CT scanner to acquire and process triple coincidences and developed a normalization procedure specific to triple events. By including triple coincidences using their method, the peak noise-equivalent count (NEC) rates of the scanner increased by 26.6 and 32% for mouse- and rat-sized objects, respectively, leading to improved image quality with better signal-to-noise ratio while preserving spatial resolution and contrast. This ML-based inter-detector scattering recovery method also outperformed other approaches based on maximum energy, Compton kinematics, and machine learning in a simulation study [32] with the scanner geometry of SimPET, a MRI-compatible PET insert [33]. Furthermore, this approach has been extended to ICS and its impact on image quality of brain PET scanners has been evaluated [34].

Fig. 9

An ML approach without energy measurement. Example PET scanner with two possible LoRs for ICS events. The orange color crystal indicates the first interaction position and the green color crystal indicates the scattered crystal position

Gross-Weege et al. [35] introduced a ML-based positioning algorithm for preclinical PET scanners. The ML algorithm compared expected light distributions with measured light distributions using Probability Density Functions (PDFs) generated from single-gamma-interaction models based on measured data. This approach led to a sensitivity gain of up to 19% over the CoG algorithm, without significantly impacting energy resolution or image quality. The study also demonstrated that the ML algorithm was less affected by missing channel information and could improve energy resolution and image quality by rejecting events that do not comply with the single-gamma-interaction model.

Lee et al. [36] treated the detector observation as a linear problem and offered two identification methods, pseudoinverse matrix calculation and convex constrained optimization and compared them with the maximum energy scheme. The study included simulation and experimental evaluations. The simulation study showed that the proposed convex optimization method yielded robust energy estimation and high ICS identification rates of 0.98 for the one-to-one coupling and individual readout with 8 × 8 3 × 3 × 20 mm3 crystal array and 0.76 for light-sharing detectors with 12 × 12 2.5 × 2.5 × 20 mm3 crystal, respectively. The experimental study showed a resolution improvement after recovering the identified ICS events into the first interaction position. The average intrinsic spatial resolutions improved by 1.5 times higher after applying ICS recovery using convex optimization method (Fig. 10). This proposed convex optimization method was applicable in typical pixelated light-sharing designs, even with a multiplexing readout scheme.

Fig. 10

The impact of ICS recovery on detector performance. The intrinsic resolution profiles of four crystals in the one-to-one coupling detector with (a) the individual readout scheme and (b) the RC sum readout scheme. (From [36], with permission)

4.3 Artificial intelligence

Artificial intelligence (AI) technology has gained significant attention and emerged as a predominant force across biomedical imaging and interpretation research [37,38,39,40,41,42]. Its applications extend beyond PET image processing, as AI technologies are actively employed to enhance the performance of PET detectors and readout technologies [3, 43, 44]. Hence many researcher groups have begun incorporating efforts to recover ICS events with various approaches using the AI technology.

Michaud et al. [45] pioneered the use of an AI-based method employing artificial neural network (ANN) to recover ICS events in PET systems. They utilized ANN on triple coincidence events to identify the correct LoR. The input of ANN consisted of position and energy information of triple coincidence data, which underwent preprocessing steps, such as normalization, energy sorting, and removal of symmetric events. A feed-forward multi-layer ANN structure was employed and trained using Monte Carlo simulation data of LabPET I [46] to acquire the best LoR recovery performance. The proposed ANN algorithm achieved a LoR recovery rate of 75 and a 55% sensitivity increase by including triple coincidences within 360–660 keV energy window. Despite some resolution degradation resulting from the inclusion of triple coincidence events rather than their rejection, the method's versatility across various conditions and pixelated detector systems renders it promising for PET imaging.

Recently, Wu et al. [47] explored the use of a feed-forward deep neural network (DNN) architecture for ICS recovery in 3D position sensitive detectors [48], which provide 3D position information: x, y, and z (depth-of-interaction, DOI). The authors utilized simulation data from a brain PET system with a 25 cm diameter and 5.4 cm length. This brain PET system consisted of detector modules arranged with 6 × 16 × 16 SiPM pixels, each 1-to-1 side-coupled with 3 × 3 × 3 mm3 LGSO crystals. The simulation assumed 12% energy resolution, an energy window of 400–600 keV, and a minimum detectable energy for any interaction of 10 keV. Here, a fully connected DNN was used with input comprising interaction position (x, y, z) and deposited energy information, capable of accommodating up to 10 scattering interactions per event. The first interaction position was used in training as label information. A softmax output layer produced a probability distribution for first interaction positions up to 10 possible positions; the position with the highest probability was then chosen. This DNN outperformed traditional algorithms such as maximum energy and lowest energy methods, achieving positioning accuracy ranging from 0.75 to 0.68 depending on the number of interactions per annihilation photons. Moreover, the DNN model showed better performance in image quality metrics such as contrast recovery (CR), CNR, and Modulation Transfer Function (MTF) values. The study suggests that DNN approach can enhance image quality and quantitation figures of merit in PET systems.

A CNN (convolutional neural network)-based approach has been also proposed to recover ICS in conventional pixelated PET detectors with various crystal sizes [49, 50]. In these studies, the authors have developed two distinct network architectures: ICS-eNet, which estimates energies in individual crystals, and ICS-cNet, which directly determines the position of the first interaction crystal (Fig. 11). Simulations were conducted to generate input data with 8 × 8 to 21 × 21 LSO crystal arrays with length of 20 mm, along with the same photosensor array of 8 × 8 covering area of 25.8 × 25.8 mm2. In the simulation study, the CNN approach showed positioning accuracy of up to 90% for 1-to-1 coupled detector design, while the accuracy decreased as the degree of light sharing increased [49]. The authors also presented experimental validation using the proposed CNN network on the system level with 17 cm diameter virtual ring PET system. The CNN network achieved improvements of 11–46%, 33–50% in volume resolution in reconstructed image, and 47–64%, for 8 × 8, 12 × 12, and 21 × 21 arrays after ICS recovery, respectively [50].

Fig. 11

Structures of ICS-eNet and ICS-cNet. (From [49]; with permission)

Petersen et al. [51] employed deep learning approach for their depth-encoding Prism-PET detector, achieving promising results for a light sharing detector. In this study, a simulated Prism-PET detector [52]—comprising a segmented light guide that enhances spatial and DOI resolution in a multi-crystal, single-ended readout— with 16 × 16 crystal array comprised of 1.4 × 1.4 × 20 mm3 LYSO coupled to 8 × 8 photosensor array was assumed. A Bayesian estimation method utilized a scatter absorption model as prior information and a detector response model as likelihood. For the machine learning approaches, they employed an ANN and a hybrid U-Net architecture (U-Net + autoencoder), which uses 8 × 8 sensor signal distribution as input. The output of the ANN was the energy and the DOI position (z) and the output of U-Net was the first interaction position and the DOI position (z), respectively. Due to the high degree of light sharing in the detector design, the Bayesian method yielded an estimation error of 20.5 keV in energy and 3.1 mm in DOI, while ANN achieved 26.2 keV in energy and 2.9 mm in DOI. By utilizing the hybrid U-Net architecture (Recovery-Net), 83% accuracy in crystal identification and 3.0 mm DOI estimation error was achieved, and image quality was significantly improved (Fig. 12). Furthermore, the results of several investigations [32, 47, 51, 53], in which the detectors with DOI information were utilized, suggest that incorporating DOI information can enhance the accuracy of ICS event position estimation.

Fig. 12

ICS recovery for a DOI detectors. (a) The count profiles were generated using three positioning approaches: centroid-based (no recovery), Recovery-Net, and an ideal detector (representing the upper limit of what's achievable). The displayed profile corresponds to a single pair of exactly opposed crystals, each positioned at the center of a Prism-PET detector module. (b) Line profiles were taken through 1.0 and 1.5 mm spots in the reconstructed image of an ultra-micro hot spot phantom. (c) The reconstructed images of the ultra-micro hot spot phantom were obtained under centroid-based, Recovery-Net, and ideal detector positioning schemes. (From [51]; with permission)