Compton imaging reconstruction methods: a comparative performance study of direct back-projection, SOE, a new Bayesian algorithm and a new Compton inversion method applied to real data with Caliste

Abstract— Compton imaging is one of the main methods to localize radioactive hotspots, which emit high-energy gamma-ray photons, above 200 keV. Most of the Compton imaging systems are composed by at least two detection layers or one 3D position sensitive detector. In this study, we demonstrate the application of a new miniature pixelated single plane detector to Compton imaging. In this configuration, we do not have the information on interaction depth but we successfully test its ability to perform Compton localization by means of comparing different Compton reconstruction algorithms applied to real data measured with our single plane detection system.

� the energy of the scattered photon, � the energy deposited by Compton interaction and the scattering angle, the energy and total relativistic momentum conservations lead to the formula: Where � � = 511 keV is the energy mass of an electron at rest.
Regarding the reconstruction method, Compton kinematics give only the angle between the direction of the original photon and the scattered photon, which gives a cone with an angle of possible origins in 3D.
Usual Compton imaging systems measure the Compton scattering interaction, which deposits � , followed by the photoelectric interaction of the scattered photon. Most of the systems use two detection layers, one scatter layer and one absorption layer, which fully absorbs the scattered photon [5]. Other systems perform Compton imaging with three or more layers, allowing reconstruction even if the final photon is not fully absorbed, assuming that the photon is scattered at least two times [6]. Finally, 3D-position sensitive systems, which are able to localize the interactions in the detector in 3D, are also used as Compton camera [7].
In this paper, we demonstrate the feasibility of Compton imaging with a single plane pixelated miniature detector, namely Caliste [8], by exploiting its advantages, especially its spectral and spatial resolution through Compton event selection. We study different Compton reconstruction algorithms to achieve the best localization performance with our detector: Direct Back-Projection (DBP) of the cones obtained by Compton kinematics, Stochastic Origin Ensemble with Resolution Recovery [9] and we introduce a new Bayesian approach and a new inversion of the Compton DBP.
Other algorithms for Compton reconstruction exist, such as Filtered Back-Projection or List-Mode Maximum Likelihood Expectation Maximization (LM-MLEM) [10]. LM-MLEM has been unsuccessfully tested for our particular detection configuration.
The paper is structured as follows: Section II focuses on the detection system properties. Section III is dedicated to the Compton events selection. Section IV describes the different algorithms we have developed and implemented for reconstruction and Section V presents the results we get from real data acquisitions.

II. CALISTE: A MINIATURE CDTE SPECTRO-IMAGER
Caliste is a hybrid CdTe detector module, initially developed for high-energy astronomy in space, declined in different versions depending on its applications. The version used for this study is Caliste-HD [11], [12].
Caliste-HD is a pixelated imaging spectrometer made of a 16×16 pixels array with 625-μm pitch patterned on a 1 cm 2 monolithic CdTe crystal, 1 mm thick. The CdTe crystal is placed on top of a 3D electronics package where eight ultra-low noise read-out ASICs IDeF-X HD [13] are embedded. In such an architecture, each pixel is a miniature and independent spectrometer channel. The package size is 10×10×16.5 mm 3 . Caliste-HD performance has an excellent spectral response with 670 eV FWHM at 60 keV and 4.1 keV at 662 keV for single events in the summed spectrum of all pixels, with an energy range from 2 keV to 1 MeV for single events. An operating temperature of -10°C and an electric field of 300 V/mm are appropriate to reach these performances. Figure 1 shows a picture of Caliste-HD. We use these characteristics to perform Compton imaging with this detector integrated in a very compact system, namely WIX-HD, with a total weight of 1 kg. WIX-HD includes one Caliste-HD module installed into a small vacuum tight enclosure connected to a data acquisition electronics system. The cooling is realized by miniature thermos-electric Peltier coolers. This system is slightly different to the usual Compton imaging gamma-ray cameras. It is indeed composed of one planar layer of detection without interaction depth measurement, which causes disadvantages in terms of angular resolution because of the uncertainty on the interaction depth. If the two deposited energies are compatible with Compton scattering, we are not able to distinguish the Compton scattering interaction and the photoelectric absorption at first glance. Moreover, the planar configuration of our detection system implies that the detector plane is a symmetric plane in the reconstructed image, so that the field of view is not 360°×180° but 180°×180° instead.
Like every Compton imaging system, the localization resolution is essentially limited by the Doppler broadening effect [14].
However, we use the other advantages of the system, its spectroscopic performances and its small size, associated to specific algorithms, in order to perform Compton reconstruction.

III. EVENTS SELECTION
The first step of Compton reconstruction is the selection of the relevant events to perform the reconstruction. In a measurement with Caliste-HD, we first select multiple events, which are events, registered in time coincidence into the detector. The time between the Compton scattering interaction and the second interaction is indeed small compared to the time resolution of the detector of 20 μs.
The second rule is that the two events do not occur in neighbouring pixels, so that there is no possible confusion with charge sharing caused by a single interaction between the two pixels [15]. The events are also accepted if they are recorded in two separated clusters of pixels. Figure 2 shows some different and typical multiple events detection patterns.
The position in 2D associated with the events is considered to be the center of the pixel, when only single pixel hits are detected such as illustrated as type a on Figure 2. Conversely, the interaction position is considered to be the centroid of the cluster weighted by the energy depositions in the cluster's pixels, see Figure 2, types b and c. For instance, in the case of events type c, the centroids of each clusters are calculated separately to determine the most probable interaction coordinates in the detector array.   Once these events are discriminated, we use energy conservation to select events with one Compton scattering interaction with an energy deposition � and a photoelectric absorption of the scattered photon � .
In our case, we assume having a prior knowledge of � as a result of a spectral source identification performed on a singleevent spectrum such as shown in Figure 3. This information can be obtained with a Bayesian Convolutional Neural Network [16], especially in case of small counting statistics. Figure 4 shows a spectrum of summed multiple events, separated by at least one pixel, for a 137 Cs source measurement. The 662 keV line corresponds to reconstructed double events with a Compton scattering in the detector followed by a photoelectric absorption of the scattered photon, so that the total energy deposited corresponds to the full energy of the original photon.
The cases of single 662 keV Compton events whose scattered photons escape from the detector form the Compton edge at 478 keV as illustrated in Figure 3. However, as expected in our simulations, 662 keV double Compton events whose second scattered photon escapes from the detector tend to shift the Compton edge towards 500 keV as shown in Figure 5. In that case, the maximum reachable energy deposition is ��� = 555 keV as expressed in eq. 3.
Whatever the Compton interaction is, single or double, the backscatter peak remains at 184 keV for 662 keV photons emitted by a 137 Cs source (eq. 4).
(eq. 4) Figure 5: Energy selection with a 137 Cs source measurement on double events with Caliste-HD. The shape of the line includes a narrow peak populated with double events type a (45% of the events), a larger bump populated with double events type b (43%) and an even broader bump populated with double events type c (12%) As shown in Figure 5, the line shape at 662 keV is composed of three event populations: the narrow peak structure is populated with events of type a (45%) summed up with two other bumpy populations from events of type b and c. The energy of the main line emission of the source is 662 keV. Because of the slightly non-linear detector response [17] and charge loss due to split events, we use a selection window, between 647 keV and 670 keV, which is considerably larger than the detector resolution for single events. The selection is made in the main peak of the spectrum so that we keep only the best events in energy to reduce the influence of the energy uncertainty on the computation of the scattering angle .
When these events are selected, we have to determine the order of the two measured interactions, i.e. which of the interactions is most likely corresponding to the Compton interaction. When one of the energies recorded is below the minimal energy of a Compton scattered photon � , this interaction is associated to the Compton scattering and the other one to the photoelectric absorption.
In this case, we distinguish the two interactions. However, the Doppler broadening effect may result in an energy of the scattered photon lower than � , which induces an error in the identification of the interaction order.
When the two deposited energies are between � and � − � , we are not able to distinguish the two interactions with this rule. This situation happens when the energy of the original photon is above � � � � � = 255.5 keV. One solution is to compute the two possibilities accepting that one is necessarily wrong, which causes noise in the reconstructed image, making the treatment more challenging. Alternatively, considering that the probability of photoelectric interaction is lower when the energy of the photon is higher, we attribute the photoelectric absorption to the lower deposited energy.
Writing � and � as the two measured energy deposits, the selection rule is: , � is attributed to the photoelectric absorption;  otherwise, � is attributed to the photoelectric absorption.
Eventually, the formula for Compton kinematics (eq. 1) can be used with different parameters, considering the energy conservation (eq. 2). We can use � + � instead of � , � − � instead of � and � − � instead of � . Analyzing the sensitivity of in (eq. 1) on the energy measurements � and � according to the different forms of the expression, we conclude that using � − � instead of � leads to the minimization of the maximal error on . Thus, we use the formula:

IV. RECONSTRUCTION ALGORITHMS
After Compton event selection, we perform Compton imaging using different algorithms. In this study, we have implemented four algorithms that we explain in this section and we evaluate their performance with our detector system on real data in Section V.

A. Direct Back-Projection (DBP)
Also called Simple Back-Projection, this algorithm is the direct application of the Compton kinematics. We note , the position of the Compton interaction, the position of the photoelectric absorption and the angle computed according to equation 5. We determine the equation of the cone of axis [ ), of apex and angle .
We compute the intersection of the cone with a sphere of radius , where is the expected distance from the detector to the source. In most cases the distance is unknown and r is set to a sufficiently high value so that we consider an infinite radius length comparatively to the detector size. This intersection is called the projection of the cone.
This algorithm is applied to each individual event and the cones are added iteratively, so that we get a back-projection map in the end. For a point source, this algorithm localizes the source at the intersections of the cones, with an uncertainty caused by the spatial and spectral resolutions of the detector.
However, when two point-sources or more with the same emission energy � are present in the field of view, the cones created by one source turns to be a noise source for the other sources, wherever they are on the reconstructed image. The problem is naturally more important when the sources are close to each other, as we will see in the results in Section V. Consequently, more advanced algorithms are required to be able to localize and separate multiple point sources emitting at the same energy.

B. Stochastic Origin Ensemble with Resolution Recovery (SOE-RR)
Stochastic Origin Ensemble [18] is a Monte-Carlo Markov Chain algorithm applied to Compton imaging. It associates each event to a pixel of the reconstructed image. For each iteration, each event might move to another pixel according to a transition probability.
For initialization of this iterative method, one particular azimuthal direction on the back-projected cone is randomly selected. The selection can be uniform along the projection of the cone or it is also possible to use a prior distribution, based on a DBP calculation for instance.
A density map ( ) is then computed on the reconstructed image, counting the number of events attributed to each pixel .
At each iteration, for each event , another position is randomly selected on the back-projected cone. � being the new selected position, � the old position, the transition probability is: In eq. 6, if the detection efficiency of each pixel is known, we can multiply this term by the ratio of the efficiency of the new and the old pixel.
After a sufficient number of iterations, depending on the counting statistics, the sources positions, and the number of sources in the field of view, the algorithm converges to a stationary state, which gives the most probable positions of the sources.
The SOE version with Resolution Recovery [9] constructs the back-projection cone randomly according to the spatial and spectral resolutions of the detector. In this work, we consider only the spatial resolution: for each interaction, we pick up a position randomly around the estimated position of the event, uniformly in a size of one detector pixel, and we randomly select the interaction depth, uniformly in the detector volume. The measured energy is kept as measured, without any modification.
Finally, all the events are processed at once, by parallelizing the computation in order to accelerate the calculation as many iteration are needed to converge, typically 20 000).

C. Bayesian algorithm
We propose a new application of a Bayesian algorithm [19] to Compton reconstruction. The idea of this algorithm is to draw each cone by accumulating the information thanks to the previous events and by applying Bayes' theorem.
On a pixelated reconstructed image, we attribute a counter ( ) to each pixel , firstly initialized to 0. We also attribute a prior probability ( ) of the origin of the photons, which could be either uniformly initialized or initialized with a prior knowledge depending on the geometry of the detector. For instance, in the case of Caliste and its direct environment, we use a uniform prior in the field of view of the detector.
For each event , we compute the following posterior probability for each pixel to be the origin of this measurement thanks to Bayes' theorem: Where ( ) is the current prior probability of observing an event originated from the pixel .
( | ) is the probability to observe the event among all of the observable Compton events originated from the pixel . The rigorous computation of this probability requires a heavy simulation of the detector response in terms of Compton events for each position of the reconstructed image. With this approach, a histogram of all recorded events is made, considering the positions of the first and the second interaction and the energy deposited in the second interaction. The bin is defined according to the detector configuration. For instance, a bin of one pixel for the two interaction positions and 1 keV in energy can be used.
However, this simulation is very long to compute to get a sufficient number of examples. For � = 662 keV, � − � = 478 and considering that neighbouring interactions are excluded, the number of bins required is: 256 × (256 − 9) × 478 = 30 224 896 This is the order of magnitude of the number of events to simulate in order to get a fair approximation of ( | ). Moreover, this simulation must be performed for every pixel , i.e., 2025 times for a pixelization of 4°×4° in a field of view of 180°×180°.
Consequently, we chose to approximate ( | ) by considering the function ( | ), which is chosen equal to 1 in the pixels belonging to the cone defined by and equal to 0 for the other pixels. This approach is similar to the approximation made in some implementations of LM-MLEM [10].
In eq. 7, ( ) is the probability to observe the event . Its computation is not necessary by using the rule that ∑ ( | ) � =
We update the prior probability with: This update is made after accumulating some events to avoid that the first events put a probability of 0 to some pixels. We empirically chose to update ( ) after processing 10 % of the total number of events.
After processing all events, ( ) gives a proportion map of the presence of radioactive sources.
Reference [19] gives more details on possible variations of this Bayesian based algorithm.

D. Compton inversion
The last algorithm we propose is a Compton inversion algorithm of the image initially computed by a Direct Back-Projection method.
A passing matrix is obtained by simulating the presence of a source in each pixel of the source space image, and applying the DBP to these simulated data, so that we derive, at each possible position in space, the expected DBP pixelated with the same resolution as the source space image.
The line of contains the simulated DBP, in flattened data, of a source placed in the pixel of the source space image.
If we note , the DBP of a real acquisition and , the expected intensities of the source in each pixel of the source space image, the equation to be solved is written as: is a matrix with high dimensions, of a typical size of 8100 2 for a field of view of 180°×180° with a discretization of 2°×2°. Its inversion is not feasible; therefore, we may approach an optimal solution of this inverse problem minimizing instead: The minimization process is constrained to positive elements of . A gradient descent algorithm computes this minimization, since this is a convex problem.
Once the optimal solution is found, we apply a TVregularization [20], which is used in Compton imaging with other algorithms [21], in order to denoise the source space image by minimizing ‖∇ ‖ � .

A. Test configuration
We take real data with Caliste-HD to apply and test the algorithms. A 1.7 MBq 137 Cs source is placed in a distance of 30 cm from the detector at different longitudes and latitudes. The corresponding dose rate is ~1.4 μSv/h at the detector level.
Two separate measurements are acquired: one measurement with a source directly placed on axis, at longitude = 0° and latitude = 0° for a total exposure time of 35 h. The corresponding count rate at Caliste level is 50 counts/s/cm 2 . A second measurement with the source moved to the right at longitude = 10° and = 0° for a total exposure time of 57 h. Count rate is the same. The coordinates system is explained on Figure 6. From these data, we can build a third augmented data set combining randomly all events from the two independent measurements into a single list.
The count rate being low, (50 counts/s at 30 cm), the chance of spurious double events from the two sources -that we cannot record combining the two event lists -is considered to be negligible (<< 1 count/s). This way, we obtain a representative data set that we would have recorded with two 137 Cs sources simultaneously installed in the field of view the detector. Choosing the number of events in the two lists, allows to virtually modulate the sources relative intensities. It also allows to test our methods with an arbitrary number of sources in the field of view even having only one source available in the lab.
In the results section, we study the Compton reconstruction processed on the measurement with the source at coordinates (10°,0°) and the "augmented data set" with two sources at coordinates (0°,0°) and (10°,0°). Other configurations have been successfully tested with similar performance of the algorithms.
The algorithms are applied with a field of view of 180°×180°. To help the reader to see the results, the images of Figure 7A-D are zoomed in the region of interest, between -30° and 30° for the longitude and between -15° and 15° for the latitude. Based on experience, we apply a discretization of 2°×2° for DBP and Compton inversion, 4°×4° for the Bayesian algorithm and 5°×5° for SOE-RR. We apply a spline interpolation for the representation on the reconstructed image. SOE-RR is performed with 20 000 iterations.

B. Results
The total number of detected events with the acquisition at (10°,0°) is 5 423 829. Applying the energy selection according to the method described in section III, 8119 Compton events are extracted, corresponding to 0.15 % of the total number of recorded events.
With the two sources mix augmented data set, the total number of events is 8 614 247 and 13 076 Compton events are extracted, also corresponding to 0.15 % of the events.   Computation time is evaluated with the following hardware setup: Intel® Xeon® Silver 4114, CPU @ 2,20 GHz. The algorithms are implemented in Python using NumPy library.

C. Comments
The Direct Back-Projection algorithm shows the worst performances in terms of angular resolution. However, it shows the fastest computation time, 12 s to process 4500 events. SOE-RR shows a way better localization performance, especially with the power to separate two sources 10° apart. However, the main limitation is related to the computation time, nine to ten times slower than DBP. Moreover, as the number of iterations to compute is not pre-determined, the computation time depends on the number of iterations. Nevertheless, this algorithm considers the uncertainties in the measurement thanks to the Resolution Recovery version, which is important to get a more robust result.
The Bayesian algorithm shows intermediate performances in terms of localization. On Figure 7C right, the left source is not localized at the exact position: this is a discretization effect in the 4°×4° pixelization reconstructed image. In fact, the probability of presence is higher in the neighbouring pixel of the real position of the source. The two sources are however well separated, with a quick computation time similar to the Direct Back-Projection algorithm. Actually, the number of operations is the same as in DBP, since the list of events is scanned only once.
We conclude here that SOE-RR is the most promising algorithm to solve our Compton reconstruction problem with a single layer Compton imager as long as computation time is not an issue. Our Bayesian algorithm gives a fast and relatively accurate result, possibly used in other methods as a better prior than the Direct Back-Projection reconstruction.
Compton inversion algorithm is less convincing than SOE-RR and the Bayesian algorithm for the detection of point sources. However, a preliminary study based on simulations shows that this algorithm will reach better performances in the case of extended sources.

VI. CONCLUSION AND OUTLOOKS
Compton imaging can be performed with a single plane miniature pixelated detector with fine pitch and high energy resolution. It requires a particular attention on the selection of the Compton events, given the configuration and the characteristics of the Caliste-HD detector, followed by a reconstruction algorithm. We tested four algorithms on real data and demonstrated that SOE-RR and a new Bayesian algorithm, are able to separate two sources of 137 Cs separated by only10° which a classical Direct Back Projection cannot do.
In forthcoming studies, we envision to further develop the advanced SOE-RR algorithms, noticeably including the energy resolution and iterative event selection. As computation time might be a severe limitation for use in real time applications, we will cascade several algorithms using one as a prior for the next. An evaluation of the sensitivity of the system must be carried out, since we use here the whole statistics of acquired photons leading to more than 8000 Compton events, but it is important to evaluate the necessary number of measured Compton events to perform the reconstruction algorithms. In addition, the case of extended sources will be studied, with real data acquisition, in order to evaluate the capability of our system and algorithms to image extended sources.