Compton-based prompt gamma imaging using ordered origin ensemble algorithm with resolution recovery in proton therapy

Prompt gamma ray (PG) imaging based on Compton camera (CC) is promising to realize in vivo verification during the proton therapy. However, the finite spatial and energy resolution of current CC, as well as the Doppler broaden effect, degrade the quality and resolution of PG images. In addition, due to the inherent geometrical complexity of Compton camera data, PG imaging can be time-consuming and difficult to reconstruct in real-time, while using standard techniques such as filtered back-projection or maximum likelihood-expectation maximization. In this paper, we propose three modifications of origin ensembles with resolution recovery (OE-RR) algorithm based on Markov chains to accelerate the convergence to equilibrium of OE-RR algorithm and improve the image quality. For evaluation, we performed a Monte Carlo simulation of a three-stage CZT Compton camera with resolution loss to detect the PG produced by a proton beam in a water phantom, and evaluate image quality of the gamma rays emitted during proton irradiation. The results show that our ordered OE-RR algorithm realized a good resolution recovery and accurate estimation of the position, including the peak and the distal falloff of the PG emission with remarkably faster reconstruction, thus demonstrating the feasibility of this new method in non-idealized PG-based proton range verification.

Currently, several methods have been proposed to solve the problem. The list-mode ordered subset expectation maximization (OSEM) with the shift-variant point spread functions (LM-OSEM-SV-PSFs) can recover the reconstructed degrade image by incorporating the SV-PSFs into the iteration process 9 . And this algorithm worked well for the reconstruction of proton-induced PGs 10 . However, the SV-PSF parameters estimation is time-consuming and have to change as the detection environment (including the detector geometry, field of view (FOV), source characteristics) changes, making the clinical real-time image reconstruction difficult to achieve. Besides, one previous study showed that the origin ensemble with resolution recovery (OE-RR) algorithm based on Markov chain 8 , which is an extended (stochastic origin ensembles) SOE algorithm 11 including resolution recovery (RR), has good performance in terms of image quality while clearly outperforming in reconstruction time. In their research, 511 keV multiple point-like sources with four groups CC were simulated and used to demonstrate the feasibility of this algorithm in resolution correction. However, the convergence to equilibrium of OE-RR was slower than the original SOE, and the store of the coefficients' matrix presented a challenge in the real-time reconstruction of large amounts of data such as PG. Meanwhile, as far as we know, no research has used the OE-RR algorithm to reconstruct the PGs induced by the proton beam.
In this study, PG induced by proton pencil beam were simulated and detected by a three-stage Compton camera with finite energy and spatial resolution. We assess the performance of PGs images reconstructed by OE-RR with Monte Carlo simulated data. We evaluate how well OE-RR could reproduce the peak and distal falloff of the proton beam, as well as how well it could improve the spatial resolution of reconstructions. In order to speed up the reconstruction by OE-RR and optimize the image quality, we made three modifications of OE-RR to improve its performance for PGI. Both the effects of the modifications were also evaluated. To evaluate the feasibility of range verification with non-idealized CC, a large water phantom (20 cm × 20 cm × 20 cm) similar to patients' size were chosen. Besides, the different reconstructions of characteristic PG were also evaluated.

Results
Origin ensembles based reconstruction. The PG emissions were produced by a 120 MeV proton pencil beam (2D Gaussian spatial profile σ = 3 mm, total number = 1.0 × 10 10 ), irradiating a large water phantom (20 cm length × 20 cm width × 20 cm depth) in our simulation. The 2D projection PG images were based on the volume of interest (VOI) 2000 × 2000 × 1 bins, the size of each bin being 0.1 mm × 0.1 mm × 0.1 mm. The center of the VOI in depth along the center of the phantom (center of the beam). They were 2D reconstruction since our reconstructions is to evaluate the difference between the reconstruction and the true value (Monte Carlo simulations) in the same spatial location (i.e. the same plane). The horizontal profiles were taken the central line of 2D projection images along the direction and perpendicular to the direction of proton beam, respectively. The horizontal profiles were created by taking the central line of the 2D image along the beam direction and transverse direction. The 2D projection images used for obtaining the horizontal profiles were followed by a 2D Gaussian post-smoothing filter with a small full-width half-maximum (FWHM) 2.0 mm to reduce the background noise and reproduce peaks better. The peak positions were obtained by recording the location of the point with a value of "1". The 80% and 50% falloff positions are obtained by linear interpolation with a maximum scale 0.1 mm. Three repeated independent simulations were implemented to reduce random error. Figure 1 shows the profiles of 2D projection 16 O PG images reconstructed by the SOE 11 and OE with event ordered modification 2 . The number of events used for reconstruction in Fig. 1 is about 70,000. Thus, 10 times ordered OE iterations are equal to 700,000 stochastic OE iterations with regard to the numbers of OE iterations. For a better comparison of OE algorithm, the interaction positions and deposited energies of events were exact, regardless of detection error. The narrower width at 80% positions and better reproduction of falloff and peak positions obtained by ordered OE in Fig. 1 and Table 1, show that event ordered used for OE iterations could provide a better reconstruction compared with event stochastic used. Besides, as shown in Table 1, the modification for event ordered can accelerate the reconstruction of OE algorithm due to its simpler iterative process without randomly selecting events.  16  As Table 1 shows, 100 times iterations for ordered OE algorithm with 10 4 events were sufficient to obtain the accurate estimation and high spatial resolution for PG imaging. Thus, the SOE without resolution recovery (abbreviated as SOE) and OE-RR in our study were ordered event used to accelerate the reconstruction process.    Figure 2 shows the horizontal profiles of the 2D projection 16 O PG images reconstructed by SOE, three kinds of OE-RR, as well as obtained by Monte Carlo simulation. The three kinds of OE-RR consist of OERR with the corrections of pre-calculation and deposited energies (OERR E,Pre ), OERR with the correction of the deposited energies (OERR E ) and OERR without the corrections of the pre-calculation and deposited energies (OERR). The method of pre-calculation correction for OERR E,Pre will be presented in section 4.2 (i.e. pre-calculation for initial states of OE-RR), and the method of deposited energies correction for OERR E,Pre and OERR E will be presented in section 4.3. The detection of PG used a three stage Cadmium Zinc Telluride (CZT) detector with spatial resolution of 1 mm in landscape and 1 mm in the depth, energy resolution of 1.64% for 137 Cs (662 keV). The total coincidence events used for reconstruction was about 137000, 123000, 70000 and 75000 for 14 N, 12 C, 15 O and 16 O de-excitations, respectively. Figure 3 compares the horizontal profiles of the 2D projection different characteristic PG images obtained by origin ensembles (OE) algorithms and Monte Carlo simulation.
As shown in Figs 2 and 3, it can be observed that the FWHM or the width at 80% obtained by both three kinds of OE-RR is less than that of SOE. Besides, compared with SOE, the peak positions and fall-off distal obtained by OE-RR have a better agreement with the true value (i.e. Monte Carlo simulation). Thus, the OE-RR reconstructions provide better spatial resolution while the SOE reconstructions are distorted due to the finite spatial and energy resolution of CC. The improvement by OE-RR is shown in Tables 2 and 3, specifically.  Table 3. The width at 80% positions(mm) obtained by the Monte Carlo simulation (true), SOE, OE-RR with the corrections of pre-calculation and deposited energies (OERR E,Pre ), OERR with the correction of the deposited energies (OERR E ) and OERR without the corrections of the pre-calculation and deposited energies (OERR). Values in the table represent the mean (mm) ± standard deviation (mm) of three repeated independent simulations. Moreover, even though the reconstructions obtained by three OERRs were similar, the peak positions reconstructed by OERR E,Pre were generally closer to that obtained by Monte Carlo simulation such as the 14 N + 12 C + 15 O + 16 O PGs within 2 mm error in Table 2. In addition, as shown in Table 3, the width at 80% positions along the beam direction obtained by the OERR E,Pre are less than the others for 12 C, 14 Table 2 shows, for the large phantom simulation, both the reconstructions by OE-RR algorithms could provide accurate estimation for the 80% falloff positions of PGs distribution within 2 mm. Besides, the reconstruction by OERR E,Pre could predict the peak positions of 12 C, 15 O, 16 O and 14 N + 12 C + 15 O + 16 O PGs obtained by Monte Carlo simulation with an accuracy of less than 1.9 mm. Moreover, OERR E,Pre based reconstructions reproduced the peak position and the falloff of 16 O with the accuracy of less than 1.65 mm. These results demonstrate the feasibility of OERR E,Pre based reconstruction for range verification with large phantom.
The results of OERR E and OERR in Table 3 show that the modification of deposited energies worked well for high energy PG like 16 O (6.13 MeV), but for lower energy such as 14 N (2.31 MeV), the modification has little effects. Besides, by comparing the results of OERR E,Pre and OERR E in Tables 2 and 3, we find that the modification of pre-calculation of probability densities in VOI performed better for 12 C, 15 O and 16 O since their true PG distribution were more concentrated. That is the OERR E,Pre led more convergent results compared with OERR E . However, for 14 N PG, whose distribution is like a line. The peak position obtained by OERR E,Pre is worse than that obtained by OERR E . Considering that 14 N PG distribution is in worse agreement with Bragg peak compared with 12 C, 15 O and 16 O, OERR E,Pre is more suitable for PGs reconstruction.     Reconstructed PGs image. The water phantom (10 cm length × 10 cm width × 20 cm depth) were chosen for further PG measurement and evaluation of reconstructions. A 120 MeV proton pencil beam (2D Gaussian spatial profile σ = 5 mm, total number = 1 × 10 9 ) was irradiated to the phantom to simulate the small area treatment in proton therapy. The three-stage CZT detectors with spatial resolution of 2 mm in landscape and 1 mm in the depth direction, energy resolution of 1.64% for 137 Cs (662 keV) were used. The number of coincidence triple-events used for reconstruction was about 184 000, and the number of coincidence double-events was about 64 000. The total coincidence events used for reconstruction was about 88 000, 74 000, 41 000 and 45 000 for 14 N, 12 Fig. 10, respectively. The results show that the distal edges of the OE-RR reconstructions were found in good agreement with that obtained by Monte Carlo (true value). Besides, the less FWHM of OE-RR reconstruction provided a better spatial resolution compared with SOE reconstruction. However, for PG derived from 14 N de-excitation, OE-RR could not improve the agreement between the Monte Carlo PG and reconstructed PG.
Range verification. Table 4 shows the depths of the peak position, 80% and 50% falloff positions for 14 N + 12 C + 15 O + 16 O, 14 N, 12 C, 15 O and 16 O PG emissions obtained by Monte Carlo simulation, SOE without RR and OE-RR. Compared with the reconstruction using SOE, OE-RR-based reconstruction provided better estimates of both the peak position, 80% and 50% falloff positions. Compared with the peak position, the 80% and 50% falloff positions can be more accurately estimated for OE-based reconstruction. Our results also show that the reconstruction using OE-RR can be used to predict the 80% and 50% falloff positions of 14 N + 12 C + 15 O + 16 O, 12 C, 15 O and 16 O PG emissions with an accuracy of less than 1.6 mm. For PG emission derived from 14 N, the difference between the peak position and 50% falloff positions of reconstructed by OE-RR and that of Monte Carlo simulation was large (>6 mm). However, the OE-RR reconstruction can predict the 80% falloff positions of it to within 1.4 mm. The mean errors of the peak, 80% and 50% falloff positions(mm) of PG emissions between the Monte Carlo simulation (true) and that calculated by SOE and OE-RR are shown in Table 5.  To further evaluate the performance of the OE-RR algorithm, we also studied the reconstructions for different count levels of proton. As shown in Table 6, protons from 10 8 to 10 10 were simulated. for 14 N + 12 C + 15 O + 16 O, 12 C, 15 O and 16 O PG emissions, it's observed that the reconstructions were similar when we changed the count levels of proton. For 14 N PG emission, as the count level increases, the estimated peak positions obtained by OE-RR are closer to the true values. The results in Table 6 also indicate that 12 C, 15 O and 16 O PG emissions reconstructed by OE-RR were more suited for in proton beam range monitoring due to their higher accuracy of reconstruction and less dependence on the numbers of incident protons. Since the falloff positions at 80% and 50% predicted by using different count levels of protons have an accuracy of less than 1.5 mm, it is proved that the reconstruction of OE-RR algorithm can be used to accurately predicted the 80% and 50% falloff positions for PG emissions induced by different count levels of protons.
Reconstruction time. For this study, a 64-bit Linux computer, with a 2.50 GHz Intel i5-7200U CPU, was used to run the OE algorithm written in C++. For 100 times ordered OE iterations for 248 000 events (i.e. 24.8 million OE iterations), the reconstruction time for OE-RR was 28 s while the reconstruction time was 10 s for SOE algorithm. Both of them were utilized one thread with a single core.

Discussion
In this paper, we studied the origin ensembles based reconstruction for PG imaging induced by proton pencil beam and detected by a three-stage CC in proton therapy from simulation. We made three modifications for OE-RR and investigated the improvement of each modification. Finally, we used the optimized OE-RR algorithm for PG reconstruction and evaluated its performance. Like previous studies 4, 9 , the results of PG imaging show that the spatial resolution degradation caused by finite spatial and energy resolution of detection system reduces the accuracy of range verification based on reconstruction. The results also show that the optimized OE-RR algorithm realized better resolution recovery and more accurate estimation in range verification compared with SOE algorithm, as well as demonstrate the feasibility of this novel method in non-idealized PG-based proton range verification. For 14 N + 12 C + 15 O + 16 O PG emissions, the OE-RR based reconstructions provided a better spatial resolution and more accurate estimation for the peak and falloff positions (i.e. 4.1 mm to 1.9 mm in Table 4) compared with the reconstructions obtained by SOE algorithm. Besides, the OE-RR based reconstructions can be used to predict the 14 N + 12 C + 15 O + 16 O PGs within 0.2 mm at 80% falloff positions and 0.5 mm at 50% falloff positions, respectively. Thus, it's promising to monitor the beam range during proton therapy due to the accuracy of within 2 mm 4,12 .
As can be seen from the results, the PG produced by different kinds of atomic de-excitations (e.g. 14 N, 12 C, 15 O, 16 O) showed different correlations with the proton depth-dose curve, corresponding a previous study 13 . Thus, we also evaluated the reconstruction images of each characteristic, respectively. Our results show the OE-RR reconstruction can be used to predict the 80% and 50% falloff positions of 12 C, 15 O and 16 O PG emissions within 1.2 mm or less, as well as predicted the peak position with an accuracy of less than 1.6 mm. Although the reconstructed 14 N PG had a poor agreement with that obtained by Monte Carlo simulation at the peak position, the falloff at 80% and 50% positions were accurately predicted within 1.4 mm. According to one recent study 13 , the PGs derived from 12 C and 15 O de-excitation (i.e. 4.44 MeV and 5.25 MeV) have the similar distribution to the proton depth-dose curve, while the PGs derived from 16 O de-excitation (i.e. 6.13 MeV) have the closest falloff correlation with the dose deposition curve in proton therapy. Therefore, the reconstruction of PG using OE-RR is able to accurately monitor in vivo dose deposition in proton therapy.
For the reconstruction of multi-energy PGs (e.g. 14 N + 12 C + 15 O + 16 O), the reconstruction process of the algorithm is the same as that of a single energy. However, the reconstructed images of multi-energy had a lower agreement with the Monte Carlo distribution, compared to that using single characteristic peak energy window (e.g. 12 C, 15 O or 16 O). Since the 14 N PG distribution had a lower agreement with dose deposition of proton beam compared with that of 12 C and 15 O 13 , the reconstruction contained 14 N PG degraded the image quality and the accuracy of the image. But the good agreement of 80% and 50% falloff positions (within 1.4 mm) demonstrates the feasibility of this optimized OE-RR for multi-energy events, as well as the multiple interactions in detection system(e.g. triple events).
The reconstruction time for OE-RR is more than that of SOE algorithm due to the correction operation for deposited energy and intersected positions in per iteration (e.g. N 1 × 11 times correction for N 1 triple events, and N 2 × 8 times correction for N 2 double events, respectively). However, the speed of OE-based reconstruction was more than five times faster than PSF-based algorithm 10 , and the quality of reconstruction has good agreement shown above. Moreover, the reconstruction of OE-RR could further speed up by using parallel computing based on graphics processing unit (GPU), enabling the OE-RR algorithm to realize real-time image of PG.
Our results also indicated the requirement of the proposed method in the non-idealized system. For the large phantom used in our simulation, which is closer to the clinical application, the Compton camera (CC) is suggested to have the spatial resolution of less than 1 mm (i.e. the largest size of pixels is less than 2 mm × 2 mm), as well as the high energies resolution (i.e. semiconductor detectors like CZT or high purity germanium is more suitable due to their measurement error of deposited energies is within 1% for PGs). For the small phantom which can be detected with less distance between the PGs and detectors (i.e. nearly at 15 cm), the spatial resolution equals 2 mm of CC is sufficient for the accurate estimation in range verification. However, in both cases, the energy resolution of CC is suggested to within 1% for accurate reconstructions. Since the CZT detectors can work at room temperature without coolant, the size of the detection device can be less in volume and easier to move during treatment compared with that made by scintillation detector. The most complicated part in CC imaging is the processing of electronics and the records of coincident events. But we believe the problem can be solved well in real-world according the experiments implemented by Polf, J. C et al. 4 . As for the imaging geometry or setup, since the reconstruction time of OE-RR is dependent on the number of detected events but there is almost no relation between the geometrical complexity of the imaging space (i.e. volume of interest) and the reconstruction time, the proposed method is promising to realize the real-time 3D   reconstructions regardless of the choice about the size or bins. Besides, since the computing burden of OE reconstructions is only dependent on the number of events and generation of random numbers, the requirement of the CPU or storage capacity of imaging equipment is easy to meet. The difficulties for different treatment sites are generally dependent on the distance between the CC (i.e. the first scattering detector) and the treatment site for the proposed method. The detection efficiency decreases as the distance increases, while the imaging quality also decreases as the distance increase due to the spatial degraded performance of cone intersection in CC imaging. Thus, the distance within 25 cm in real-world monitoring should be better. For the monitoring at longer distance, higher spatial (less than 1 mm) and energy resolution (less than 1%) of detectors may be required.

Materials and Method
Simulation and Detection of the PGs. The PGs emission and detection system were simulated by using Gate version 8.0 and Geant4 version 10.03 with QGSP_INCLXX physics list. The physics list QGSP_INCLXX, an experimental list, uses the Liege model to describe the inelastic interaction of protons and neutrons and it is suitable for simulating proton-induced interaction 14 . The standard physics option 3 was used to model the physical processes 15 . In the simulation, a 120 MeV proton pencil beam with a two-dimensional Gaussian spatial profile (σ = 3 mm~5 mm) was used, irradiating a large box water phantom (20 cm length × 20 cm width × 20 cm Depth) or a small water phantom (10 cm length × 10 cm width × 20 cm Depth). From 10 8 to 10 10 protons were simulated 16 .
To detect PGs, we use a three-stage CC, consisting of three separate detection stages, and each is composed of several pixelated Cadmium Zinc Telluride (CZT) crystals. For this study, CZT was chosen for its high interaction cross section for gamma rays in our range of interest (up to ~7 MeV), as well as its high spatial and energy resolution, and the ability to record the depth information of interaction position by the pulse height and shape analysis of the anode signal 4 . Moreover, the three-stage structure has higher detection efficiency compared to two-stage for PG, because it can record both double-events (one scatter in the first detector and then are absorbed by the second detector) and triple-events (two times scatter in two different detectors and third interact in the other detector). For the detection the PGs from the small phantom, the first two stages contain a 25 × 25 array of pixelated CZT crystals (2 mm × 2 mm × 15 mm each) for a total CZT area of 5 cm × 5 cm, and the third stage contains a 50 × 50 array with the same pixelated CZT crystals for a total CZT area of 10 cm × 10 cm (As shown in Fig. 11). For the detection the PGs from the large phantom with the same total area of CZT, the pixelated CZT crystals were both set as 1 mm × 1 mm × 15 mm each due to the higher spatial resolution requirement of detectors for the farther detection distance. Adjacent stages ware separated from 5 cm. It was assumed that both detectors had a resolution of 1 mm in landscape and 1 mm in the depth direction. The CZT detectors' energy resolution, δ CZT (FWHM) is obtained by 17 : where E is the deposited energy, σ E is the standard deviation of spectral peak. For this study, both diploid events and triple events were recorded. The a and b parameters used in our simulation are −0.2892 and 0.4332, respectively. Both the parameters Referred to CZT pixel array detectors produced by IMDETEK Corporation Ltd, whose energy resolution is 5.13% for 241 Am (59.5 keV) and 1.64% for 137 Cs (662 keV) 17 . We performed an energy cut to only detect PGs with a total energy (determined using the Compton energy formulas 2,16 in the range above 1 MeV to reject 511 keV annihilation gammas and most of background radiation. The simulation results of energy deposited of proton and PG are shown below. As shown in Fig. 12, the relative intensity of PG decreases rapidly at the end of proton range and presents a close correspondence with the distal fall-off of BP, proving that the simulation results are reasonable. As shown in Fig. 13, four main PG energy spectral peaks are: 2.31 MeV ( 14 N de-excitation), 4.44 MeV ( 12 C de-excitation), 5.25 MeV ( 15 O de-excitation) and 6.13 MeV ( 16 O de-excitation) 18 . Since PG derived from 14 N, 12 C, 15 O and 16 O de-excitation exhibits different correlations with proton energy deposited in depth 13 , we evaluated the reconstruction by them respectively. The total energy windows of coincidence events were set as within ±0.2 MeV of the four known PG energy spectral peaks, to alleviate the affect due to the incomplete absorption and background radiation. For our designed Compton camera, the efficiency was 4.1 × 10 −4 (E total > 1 MeV) events per incident proton.

Optimized OE-RR algorithm. The SOE algorithm is a Monte Carlo Markov chain method which uses the
Metropolis-Hastings algorithm 19 . Rather than back-project the entire cone defined by detected gammas in the reconstruction, the SOE algorithm reproduces the image by considering only a single representative point on the conical surface. The representative points are initially chosen by randomly selecting a point on the conical surface and within the phantom volume. Each iteration of the algorithm attempts to improve the reconstruction by exchanging the current representative points with points where the probability of a gamma originating is higher. The reconstruction time using SOE approach is dependent on the number of detected events but there is almost no relation between the geometrical complexity of the region of response and the reconstruction time. Thus, SOE has a wide horizon of applications in emission tomography. The SOE algorithm has been used for PET, SPECT and CC image reconstruction 11,20 . For Compton-based PG images, SOE algorithm performed well with ideal CC 2 , but it couldn't reproduce the PG distribution with a real-world CC because of the finite spatial and energy resolution 4 .
In the ideal case, the half conical surface corresponding to each CC event used in the OE reconstruction is defined by the exact deposited energies and the exact locations of interactions in the CC detectors. In fact, the conical surface, which is created using the measured detectors' outputs, may not include the true event origin due to the limitations of the detection systems mentioned above, creating the deterioration of the spatial resolution in the reconstructed image.
In order to correct the finite energy and spatial resolutions, rather than using the output data of the measured detector to create the cones and implement the SOE iteration, OE-RR algorithm obtained the representative points on the "guess" cones determined by randomly corrected data corresponding to the measured data, which samples from the distributions of the positions of interactions in projection elements and deposited energies. After many iterations similar to SOE using these imaginary conical surfaces, the distribution of the sources will eventually be reconstructed.
In this study, we propose three modifications of OE-RR algorithm to improve its performance for PG image, and compare performance with the algorithm as described by Andreyev et al. 11 .
The first difference is that we referred to the modification of SOE in Mackin et al. 2 , stepped through the detected gammas one after another and calculated the acceptance probabilities A(Y s → Y s+1 ) to implement ensemble transitions, where Y s and Y s+1 are two subsequent, whereas Andreyev et al. 8,11 select a detected gamma at random. Stepping sequentially through the list of detected gammas and ensuring that each gamma is tested once for each iteration improves the agreement between the MC true origin distribution and the SOE reconstructed origin distribution 2 .
The second difference is that we used random sampling method in sequence to obtain the initial event density distribution in the volume of interest (VOI) based on corrected values, whereas Andreyev et al. 8,11 used the measured value. The sampling was uniform on the conical surface (which was in VOI at the same time) defined by obtaining samples from the distribution of the locations of interaction in projection elements and deposited energies. We evenly sampled k (k ≥ 1) points in VOI for each event. Since the true value of the interacted position and deposited energy is different from the measured value, using the corrected value (or called "guess" value) to generate the initial distribution is closer to the equilibration region, making the algorithm arriving faster to the equilibrium. Moreover, the sampling process can be implemented in parallel by dividing the event into several subsets. Thus, we can obtain a reasonable initial event density distribution in VOI with little time, while the convergence to equilibrium of the algorithm is remarkably faster.
The third difference is the correction for finite energy resolution. Since the PGs' characteristic is different from the radioisotopes used in medical imaging, as well as the different detection methods, the correction method used in this study was different from that mentioned in Andreyev et al. 8 , aiming at improving the performance of OE-RR in PG image.
Correction for finite energy resolution. In order to be able to "guess" the true deposited energy based on the measured energy, the posterior distribution of the true deposited energy is used. Different from the energies of photons emitted by the radioisotopes used in medical imaging, the exact energies of proton-induced PG are always unknown because of the prompt gammas' secondary interaction in phantom. Besides, the detected events used for reconstruction contains triple events. The energy correction is implemented on the individual interactions stemming from a single source photon. The uncertainty of deposited energies in different stage were assumed to be independent. Thus, the corrected deposited energies used in our study is given as 3 for the three stages, respectively. The real scattering angle is β  . The parameters measured by the camera are the locations L 1 , L 2 , L 3 , and the energies E 1 , E 1 and E 3 , from which the Compton angle β is calculated. The point S where the photon is emitted lies on the cone having the apex in  L 1 , the axis collinear with the scattered ray direction and the half-opening angle β  .
where E i is the energy that is actually measured, i is the number of stage, σ E i is the statistical standard deviation of E i , N is normal distribution defined by E i and σ E i . For the CZT simulated in our study, σ E i was given as where η is the coefficient of proportionality between the sigma and the deposited energy is known 8 (in our case 1.5%). In fact, the likelihoods | E E p( ) i i has Gaussian shapes centered around E i and standard deviations σ E i . where the measurement E i is assumed to be derived from the Gaussian distribution with true deposited energy E i and standard deviation σ E i . We assume that σ σ ≈ E E i i , then we can obtain a reasonable "guess" of ′ E i by using (2), since the true deposited energy E i could be derived form σ N E ( , ) i E i in the case that E i is the true deposited energy. The corrected deposited energies defined by (2) is used in each Markov step of the OE algorithm to determine the cone angle.
Correction for finite spatial resolution. For the correction of spatial resolution, we modify the detection data based on the size of the detection unit in x-y plane and the measure error in z (depth) direction as shown in Fig. 14. The data measured by the detector is = L x y z ( , , ) where the x y ( , ) i i is the central coordinate of the interacted detection unit (projection element) in x-y plane. For the case that the scale of the detection unit is much smaller than the detection distance, it can be considered that the actual intersection location is uniformly distributed in the detection unit. Assuming that the spatial scale of the i-th detection unit is D i , the correction value ′ = ′ ′ ′ L x y z ( , , ) i i i i is given by, where U is the uniform distribution, δ iz is the measure error in the z direction for i-th stage.
When an event is considered by the OE algorithm, a new position of interaction within the detector unit is selected according to the distribution defined by (4) ~ (6). Unlike the correction for finite energy resolution where two random number was generated to determine a sample from the normal distribution, the implementation of the correction for spatial resolution requires a generation of six random numbers for double events, night random numbers for triple events, respectively. Three random numbers are needed to simulate the detection uncertainty in the x, y, and z directions of the detector elements for each of the detectors.
Reconstruction using OE-RR. The modifications and corrections described above were combined in the optimized OE-RR algorithm to allow for the full modeling of the Compton camera resolution due to limited detectors' resolutions. Taking the detected triple event as an example and assuming the information for each detected event that we measured is (L 1 , L 2 , L 3 , E 1 , E 2 , E 3 ), in the OE-RR reconstruction the following steps were executed: 1. For each detected event, the identifier was randomly defined when input and store the list-mode data. We used the random sampling method in sequence to obtain the initial event density distribution in the VOI based on the corrected values ′ ′ ′ ′ ′ ′ L L L E E E ( , , , Scientific RepoRts | (2019) 9:1133 | https://doi.org/10.1038/s41598-018-37623-2 moves, but on average the number of event origins in each voxel remains the same. Finally, the weighted number of events in each voxel determined the probability distribution of sources in VOI.