Multiple scatter correction for single plane Compton camera imaging in nuclear decommissioning

This paper focuses on the application of Compton cameras in nuclear decommissioning, aiming to develop efficient and affordable radiological characterization methods. Compton cameras offer advantages over traditional Anger cameras, including increased sensitivity and a larger field of view. We showcase the use of monolithic scintillation detectors in a Compton camera setup to compute two‐dimensional reconstructions of the source location of incident gamma rays. The presence of multiple scattering effects in the large, spatially non‐resolved detectors can lead to spatial artifacts in the reconstructed source activity. To address this, the paper proposes an extended forward model that accounts for double scattering interactions within a single detector element. Experimental and simulation results demonstrate the effectiveness of the proposed approach.

Emission tomography based on Compton cameras is a promising imaging approach capable of registering and reconstructing the incident angle of gamma rays.Due to their design without collimation medium, the sensitivity of Compton cameras is greatly increased compared to Anger cameras and allows measurements with a large field of view.Example usages of Compton cameras include nuclear decommissioning [1][2][3], homeland security [4,5], astronomical imaging [6] or medical imaging applications like hadron therapy or single photon emission computed tomography [7][8][9][10][11][12].Due to these inherently different applications, various sizes and geometries of Compton cameras have been considered, the typical parameters that are considered for optimization in the construction of setups are usually sensitivity (overall or in specific energy ranges), spatial resolution and energy resolution [8,13].
The radiological activity in the gamma ray energy regions relevant for nuclear decommissioning can be measured with high accuracy and even in the case of low activity or large background radiation.The setup we consider consists of several monolithic scintillation detectors (plastic and CeBr 3 crystal detectors) in which the photons interact with the detector material, leading to measurements.The emission activity is reconstructed from measurements of photons that interact in at least two detectors in the Compton camera.
The measurements and validating simulations show that a Compton camera with large detectors is capable of reaching the necessary sensitivity in the nuclear decommissioning environments, but is also prone to perturbing effects due to photon interactions in the camera that are not represented by the forward model.In particular, we found that the number of measured coincidences in which the photon scattered multiple times in the scatterer is large.The induced model mismatch leads to spatial artifacts in the reconstruction of the source activity.Some works have already considered models in which coincidences of three or more interactions are explicitly modeled, but usually for photon interactions between more than two detectors, not multiple non-distinguishable interactions in the same detector element [8,14].
We propose an extension to the usual modeling of the system matrix in Compton cameras.Our model accounts for the probability of double scattering interactions in a single detector pixel, where subsequent scattering positions can not be resolved and multiple interactions not be distinguished from single scattering from the data alone.The mathematical model is described in Section 2, the measurement setup in Section 3 and numerical results presented in Section 4.

MATHEMATICAL MODEL
The mathematical modeling of the data acquisition of a Compton camera consists in the formulation of the laws underlying the kinematics of gamma ray photons in the detector material of the camera.The corresponding probabilities have been formulated, for example, in the works [15][16][17][18].In the following, we recall this model and extend it to suit our measurement setup.

Compton transform
The measurements in a Compton camera are coincident energy depositions in two different detectors in the camera.To trigger such a coincidence, a photon of initial energy  0 has to scatter on an electron in the scintillation detector material, where it deposits a part Δ of its energy.This energy difference Δ is the value measured in the scattering detector.The scattered photon, now with remaining energy  1 =  0 − Δ and a new trajectory direction must subsequently be absorbed in a photoelectric effect in a second detector.The two energy values  1 , Δ are registered in the camera and matched to a coincidence by defining a coincidence trigger time window.
The directional information about incident photons is inferred from the two energy measurements via the Compton formula where  is the scattering angle of the photon and    2 the energy of an electron at rest.Assuming a monochromatic radiation source, meaning  0 is fixed, implies that (2.1) is a bijective relation between  and  1 , abbreviated as cos  =  1 ( 1 ,  0 ). 0 being fixed, the expected number of photons that scatter at  and are absorbed at  with energy  1 is given by the following integral transform of the emission intensity () at any source location : where the integration locus becomes a cone (, ,  1 ) = { ∶ cos  =  1 ( 1 ,  0 ) ℎ  =  − ∠( − ,  − )} with a surface element d() due to the relation (2.1).Herein, (, , ) is the probability that a scattering event at  and an absorption event at  take place.Note that if  ∈ (, ,  1 ), we can write (, , ) d =  1 (, ) 2 (, ,  1 ) d(), with  1 modeling the probability that a photon emitted at  reaches  and is scattered there, and  2 including the remaining terms that depend only on  1 , but not on  (differential cross section, attenuation between  and  and the probability to be absorbed at ).The integral operator in (2.2) defines the forward model of the Compton camera image reconstruction problem.
In contrast to many other Compton camera setups, we do not include uncertainty about the interaction positions in our forward model but rather consider entirely spatially unresolved detectors.That significantly reduces the cost of the camera and the electronics, but has to be accounted for in the forward model.The measurement does not consist of interaction positions ,  and corresponding energy values  1 , but rather only of the respective detector index and energy values.The expected number of coincidences involving scatter detector   and absorption detector   with absorption energy  1 is hence the integral of (2.2) over the whole detector volumes Figure 1 shows the Compton camera setup considered here with seven spatially unresolved scintillation detectors.In this setup, the measurement (2.3) is an integral of the spatially resolved data  1  over two cylinder volumes.The integral operator (2.2) and variants of it have been analyzed under the names Compton transform or conical Radon transform before.If the detectors are spatially resolved and complete data are available, the associated inverse problem of inverting the Compton transform, that is, solving  1  =  for  is overdetermined and explicit inversion using, for example, harmonic analysis or backprojection formulas can be considered, see ref. [19] for an overview.
In practice, analytical inversion techniques are rarely used due to the low coincidence statistics, background noise and errors in spatial and energy measurements.Since the data follows a Poisson statistic, the standard technique of solving the Compton camera image reconstruction problem is the maximum likelihood expectation maximization algorithm [15].While in spatially resolved detectors where the number of measurement bins greatly exceeds the number of measured coincidences, list-mode variants of maximum likelihood expectation maximization (MLEM) are typically used, in the underdetermined problem (2.3) the standard bin-mode MLEM is computationally more efficient.
One main complication in the Compton camera image reconstruction task is the presence of large amounts of noise in the data.Partly, the presence of this noise is uncontrollable since, for example, natural background radiation can trigger coincidences in the camera.However, a number of possible scenarios can trigger coincidences in the camera that are not described by the forward model (2.2) but still originate from the source that we want to reconstruct.The following (non-exhaustive) list of other possible events all imply coincidental energy depositions.
• photons that scatter more than once in   before being registered in   , • photons that interact more than once in   , the last of these interactions not necessarily being an absorption, • photons that interact outside the camera or in the housing or electronics of the camera before interacting in the camera, • two different incident photons in two different detectors at almost the same time, misleading the camera to register a coincidence.
Some events from these categories can be removed from the data by appropriate preprocessing steps, for example, by ensuring complete absorption, that is, taking only coincidences whose total deposited energy is close to a source emission energy (if the source nuclide is known).In Compton cameras with fairly large scintillation detectors, events from the first two categories (more than one interaction in either   or   ) are particularly prevalent.Theoretically, a photon could scatter an arbitrary number of times in   before reaching   .The forward model (2.2) can be understood as the first order Born approximation of scattering in   of the actual incident photon flux on   .
Figure 2 shows histograms of coincidences from the Compton camera considered here, generated by Monte Carlo simulations.It can be seen that the simplest coincidences with one scattering and one absorption modeled by (2.2) are very rare, multiple scattering events in either detector dominate among the measurements.
It can be argued that multiple interactions in   are not too problematic since the distribution of Compton angles  in the scatterer is not changed too much.The events with multiple scattering events in   , however, are not represented well by the forward model.In the case of several consecutive interactions of a photon in the same detector, the detector measures the total deposited energy, that is, the sum Δ of energies from several interactions.When trying to reconstruct from these events with the forward model (2.2), the scattering angles inferred from (2.1) are wrong.The histograms hence indicate that a second order Born approximation could reduce the model mismatch significantly and improve reconstruction quality.

Adjusted model for multiple scattering -a second order born approximation
We aim to extend the forward model (2.3) Herein, C1 models the differential cross section of the first scattering event, attenuation in   between  1 and  2 and the probability to scatter again at  2 .Since we do not assume spatial resolution, the expected number of (twice scattered) coincidences with absorption energy  2 is given by Note that evaluation of this forward model is of increased computational complexity compared to the the first order Born approximation, since we have to integrate over three additional dimensions.

MATERIAL AND CAMERA SETUP
The Compton camera we consider consists of a total of seven monolithic detectors.The scattering and absorption detectors are placed in one plane, the setup is hence a single plane Compton camera.The detectors are of cylindrical shape, with all symmetry axes aligned parallel to each other.Each cylinder's height and diameter are three'' both.The six outer detectors are scattering detectors made of PVT plastic.The material has a favorably large ratio of scattering attenuation coefficient to absorption coefficient in the relevant energy ranges, so that scattering is predominant in these detectors.The central detector is a CeBr 3 crystal scintillation detector from Hellma materials which acts as the scatterer.The expected energy resolution of the crystal is determined by calibration measurements to be approx.4% FWHM (full width half maximum) at a photon energy of 662 keV.
Figure 1 shows the geometry of the Compton camera and a depiction of the whole setup from the front of the camera.The energy depositions are registered by a photomultiplier and a cathode attached to one of the cylinder faces.The coincidence measurements take place between the six detector pairs (each of the six outer plastic detectors builds a pair with the CeBr 3 detector).
In order to validate the camera's measurements and differentiate between different types of interactions, surrogate measurements are computed by Monte Carlo simulations.These computations are carried out by recreating a digital version of the same camera setup in the physical Monte Carlo simulation package FLUKA.

NUMERICAL RESULTS
In order to validate the extended model taking into account first and second scattering in the scatter detector, we compare the analytical models (2.3) and (2.6) with the spectra created by simulations and real measurements.
For the computation of the forward model, we have to discretize the activity distribution () as well as the data space.For the former, we consider the activity to be constant inside each pixel of a grid imposed on a two-dimensional plane at a fixed distance of 1 m to the camera.The plane corresponds to a wall to which point sources were attached in the measurements with the real camera.We consider 137-Cs point sources which emit monochromatic gamma rays with photon energy 662 keV.Since the monolithic detectors are spatially non-resolved and  0 = 662 keV is fixed, the data space only consists of the absorption energy  (corresponding to  1 in the single scattering and  2 in the double scattering case) for every of the six detector pairs.For each pair, the energy measurement is binned in 48 intervals that scale with the energy resolution of the detector, the data hence consists of counts in 48 × 6 = 288 bins, while the reconstructed image is chosen to be of size 64 × 64.For each combination of source pixel, detector pair and energy bin we compute the forward model (2.3) for single scattering and (2.6) for double scattering.The integrals over detector volumes   and   are computed by effectively discretizing the detector volume and summing up (2.2) or (2.5), respectively, for a grid of equidistantly spaced interaction points.
Since the image reconstruction problem is underdetermined here and the data still contains noise and background radiation as described in Section 2, the reconstruction procedure needs to be regularized.Stopping the MLEM iteration early provides a form of regularization.We further choose to slightly vary the MLEM algorithm and regularize the reconstruction by imposing prior information on the image space.In order to prevent large oscillations in the reconstructed image since we know that large parts of the true source distribution are zero, we choose total variation (TV) regularization.We use the iterative algorithmproposed in ref. [20] under the name EM-TV to solve TV regularized imaging inverse problems with Poisson distributed data.
The reconstructions of a point source show the effect of the unrepresented multiple scattering effects in the data.Point sources in the center of the wall are reconstructed well with the simpler model (2.3) since the response in all six detector pairs is symmetric.Once point sources are placed outside the wall's center (see Figure 3), reconstructions using the single scattering model are typically off by several centimeters or decimeters, depending on the distance between camera and wall and the incident photon energy.With the adjusted model (2.6), this effect can be corrected for since the estimated spectrum is closer to the true dynamics in the detectors.

A C K N O W L E D G M E N T S
The authors acknowledge support from the German Federal Ministry of Education and Research BMBF (project QGRIS, 15S59431 B).
Open access funding enabled and organized by Projekt DEAL.

F I G U R E 1
In (A), the relative positions of the seven cylindrical detectors in the camera are shown.The center detector in blue is a CeBr 3 crystal detector, the six outer green detectors are polyvinyltoluene (PVT) plastic scintillators.Image (B) shows the whole camera.Section 3 further discusses this setup.

F I G U R E 2
The figure shows coincidence data from a simulated 137-Cs source with perfect energy resolution and without background in the setup described in Section 3. (A) shows the coincidences grouped by number of scattering events in the scatterer, (B) the coincidences grouped by number of interactions in the absorber.In panel (C), we show the energy deposited in the absorber detector.With model (2.2), that is, only coincidence data from a single Compton scattering in the scatterer, the minimum possible energy to be deposited in the absorber at 184 keV corresponds to a scattering angle of 180 • , indicated by the red line.Due to multiple scattering in the scatterer, many coincidences with absorber energies around and below the red lines are registered.

F I G U R E 3
The figure shows the same simulated spectrum as Figure2(C) for a 137-Cs point source and one of the six detector pairs.We now compare the simulated data with the predictions of the forward models (2.3) in (A) and (2.6) in (B) for this point source.The model that takes into account double scattering is able to capture the contribution of lower energy photons in the absorber that is not present in the single scattered data.In panels (C) and (D) we show reconstructions of a 137-Cs point source, now with data from the camera.Panel (C) is computed using the forward model (2.3), panel (D) with (2.6) instead.The white dot indicates the location of the true point source.
Scattering at  1 ,  2 ∈   with intermediate energy  1 and scattering angles  1 ,  2 , respectively, and photoelectric absorption with energy  2 at  ∈   .Since the energy measured by the detector   is the sum of both energy depositions during the scattering events, we do not have access to the intermediate energy  1 or the individual angles  1 ,  2 .Applying the Compton formula (2.1) twice and summing up gives the relation cos  1 + cos  2 =  1 ( 1 ,  0 ) +  1 ( 2 ,  1 ) = 2 −    2 scattering angles and the measured energies, meaning we only have access to the value cos  1 + cos  2 .The implied geometrical relation between the energies and interaction positions now involves two cones.The source point obeys  ∈ ( 1 ,  2 ,  1 ( 1 ,  0 )) with the intermediate energy  1 , while the first interaction point obeys  1 ∈ ( 2 , ,  1 ( 2 ,  1 )).The expected number of photons measured with absorption energy  2 and whose last two interaction positions are  2 ,  is hence the integral over the two cone surfaces and all possible intermediate energies  1 ∈ ( 2 ,  0 ): by explicitly modeling the interactions in which photons scatter twice in the scatterer   .The terms which describe the underlying kinematics are similar, but we now have to consider three interactions: