Modeling and Reconstruction Strategy for Compton Scattering Tomography with Scintillation Crystals

: The recent development of energy-resolved scintillation crystals opens the way to build novel imaging concepts based on the variable energy. Among them, Compton scattering tomography (CST) is one of the most ambitious concepts. Akin to Computerized Tomography (CT), it consists in probing the attenuation map of an object of interest using external ionizing sources but strives to exploit the scattered radiation as an imaging agent. For medical applications, the scattered radiation represents 70 to 80% when the energy of the source is larger than 100 keV and results from the Compton effect. This phenomenon stands for the collision of a photon with an electron and rules the change of course and loss of energy undergone by the photon. In this article, we propose a modeling for the scattered radiation assuming polychromatic sources such as 60Co and scintillation crystals such as LBC:Ce. Further, we design a general strategy for reconstructing the electron density of the target specimen. Our results are illustrated for toy objects.


Introduction
Invented and theorized by the Nobel medal award-winners G. Houndsfield and A.M. Cormack between 1963 and 1979, the concept of Computerized Tomography (CT) has become an essential way to investigate the inside of a human body or of any type of medium. Its principle relies on the phenomenon of attenuation, characterized by the map µ E , suffered by a photon-beam which travels through matter. This phenomenon is ruled by the Beer-Lambert law with which quantifies the loss of intensity I(x, θ, E) of the photon beam by the exponential of the total attenuation along the travel path between x and x + Tθ. For instance, in a fan-beam CT-scan, an X-ray tube illuminates a target with a given intensity. A set of detectors located outside the target, typically on an annulus, will collect the incoming photon flux and characterize the attenuation map of the medium according to the Beer-Lambert law Equation (1). Afterwards, the general (inverse) problem consists in reconstructing the attenuation map µ E from the measured beam intensities for different angular views. We refer to [1,2] for a general review on reconstruction algorithms for imaging techniques. Since the advent of CT, many imaging concepts have emerged and the need in imaging has grown. One can mention Single Photon Emission CT (SPECT), Positron Emission Tomography (PET) and also Cone-Beam CT for the imaging systems based on an ionizing source.
Due to technological limitations, the factor energy has first been ignored as a potential parameter for an imaging system. However, the development of crystals and detector fabrication technologies enabling the collection of photons and separating them by range of energies has opened the way to many more imaging architectures capable of enhancing the image quality, optimizing the acquisition process or compensating for some limitations (such as limited angle issues), see [3][4][5][6][7][8][9][10]. Furthermore, the use of high-energy X-rays has increased over the years and shows interesting properties for industrial applications, see [11][12][13]. In order to provide a reliable information for an imaging system based on Compton  In this communication, we will consider the example of LBC:Ce scintillators, which provides a better energy resolution than CeBr3 scintillators and estimates their impact on our approach. We refer to [14,15] for an exhaustive study on CeBr3. LBC:Ce scintillation crystals can now be produced with 5 cm diameter, 5cm thickness and achieve about 2.8% FWHM at 662 keV [16,17]. Furthermore, one can note that LBC:Ce is self-radioactive but only in the range from 1550 keV to 2250 keV [18] which is over the measurement range (355 keV, 1332 keV) considered here. Table 1 delivers the characteristics of LBC:Ce and CeBr3 crystals which motivates the use of LBC:Ce crystals for our application, see [17]. Regarding the density of the crystals, we should properly distinguish between crystals and scanning objects. Indeed, both object and crystal will scatter the radiation to some extent. Since this study focuses on the scattering induced by the object, the crystals are assumed to perfectly absorb the scattered radiation. The physical interactions between photons and matter can be distinguished into: Thomson-Rayleigh scattering, photoelectric absorption, Compton scattering and pair production. In the classic range of applications of the X-rays or γ-rays, i.e., >60 keV, the photoelectric absorption and the Compton scattering are the dominant phenomena, see [19]. While the photoelectric absorption plays an important role in the attenuation of the photon beam, a measured photon either suffers no interaction (primary radiation) or is scattered (scattered radiation); this is the reason why the Compton scattering is more natural to exploit than the photoelectric effect as an imaging agent.
The Compton effect stands for the collision of a photon with an electron. The photon transfers a part of its energy E 0 to the electron which experiences a recoil while the photon is scattered by an (scattering) angle ω with the axis of propagation, see Figure 1. The energy of the photon after scattering is expressed by the Compton formula [20], where mc 2 = 511 keV represents the energy of an electron at rest. Usually energy resolved cameras are collimated and combined with polychromatic sources, leading to multi-channel CT-data and corresponding estimations of the attenuation map. Here, we intend to focus instead on the phenomenon of Compton scattering which enables a modeling of the energetic data in terms of electron density and for corresponding emission energies, see [21][22][23][24][25][26][27][28][29][30][31][32][33][34][35][36]. The purpose of this communication is to validate the possibility to use scintillation crystals and their energetic sensitivity in a fan-beam CT scan from a bichromatic ionizing source (here the Cobalt-60). To achieve this, the modeling and the reconstruction strategy are presented in Section 2. Section 3 presents the architecture and properties of the designed scanner with synthetic data as well as reconstructions. The validation of our approach with the properties of the LBC:Ce scintillators motivates the potential of the imaging system in particular with better resolved crystals or semiconductor detectors such as CZT (CdZnTe) [37].

Modeling and Processing of the Compton Scattered Data
We first assume the source to be monochromatic, i.e., it emits photons with same energy E 0 , for the sake of clarity. For sufficiently large E 0 , larger than 80 keV in medical applications for instance, the Compton effect represents a substantial part of the radiation as more than 70% of the emitted radiation is scattered within the whole body. As depicted by Figure 2, the Compton effect ruled by Equation (2) brings an interesting diversity into the measured spectrum as a monochromatic source leads to a polychromatic radiation measured by scintillation crystals. The development of scintillation crystals such as LBC:Ce and the improvement of their energy resolution delivers a new dimension to explore.
Neglecting Thomson-Rayleigh scattering and the pair production, we focus on the Compton scattering to interpret our data and decompose the spectrum Spec(E, s, d) measured at a detector d with energy E as follows In this equation, g 0 represents the primary radiation which crossed the object without being subject to the Compton effect. It corresponds to the signal measured in conventional CT, Equation (1). The functions g n (E, s, d) correspond to the photons that were measured at d with incoming energy E after n scattering events. The schematic curve with marks • depicted in Figure 2 (right), often called Compton scattered spectrum, is typical from what one observes when measuring the spectrum of a monochromatic source, see for instance [38]. We refer to [39] for a detailed review on the technology of detection.
first-order second-order third-order

E0
Compton effect (Right) a schematic spectrum of the detector (•) is illustrated via its decomposition in terms of type of scattered data and depicts the response to a monochromatic source due to the Compton effect.
In the following, we study how the scattered radiation g n , n = 1, 2, . . . behaves and explain how to deal with it. To simplify the notations, we denote E n+1 = E (E n , ω n+1 ) with E n the scattered energy after n scattering events and ω n+1 the (n + 1)th scattering angle.

The Scattered Flux
We first focus on the first-order scattering as depicted in Figures 1 and 3 for the notations. The travel of a scattered photon beam follows the same scheme: • The photon beam is emitted by the point source s in a differential solid angle dΩ c with an energy E 0 ; • Photons may be absorbed or scattered along the path following the Beer-Lambert law (Equation (1)) leading to the attenuation of the beam by the weight A E 0 (s, x); • A part of the beam may collide with an electron at position x 1 which belongs to a differential volume dx; • Due to Compton scattering, the photon is scattered by an angle ω 1 within a solid angle dΩ c which follows the Klein-Nishina probability dσ • Again the scattered beam is attenuated by a weight A E 1 (x, d) due to absorption or scattering of higher order; • The attenuated and scattered beam finally reaches the point detector at d.
This scheme leads to model the variation of the number of photons g 1 scattered at x and detected at d with energy E 1 by The second variable in the intensity I allows the source to be anisotropic. This formula describes the evolution of the first scattered radiation which is detected at a given energy and at a given detector position. Regarding the first-order scattering, the Compton Formula (2) relates the energy to the scattering angle and thus delivers a unique geometry for each energy measured. Indeed, all scattering points x responsible for a detected scattered photon at energy E 1 belong to T(ω 1 , s, d) defines in 3D the lemon part of a spindle torus (ω 1 < π/2) and the apple part of a spindle torus (ω 1 > π/2). In 2D, this reduces to two circular-arcs as depicted in Figure 4a. Figure 4b shows the scanning of a specimen by paired circular-arcs for different energy levels or equivalently different scattering angles ω 1 . An integral modeling of g 1 can then be obtained by integrating Equation (4) over the whole domain of the specimen and leads to a generalized Radon transform along T(ω 1 , s, d), noted C 1 (n e ), We refer to [41,42] for more details.
As explained above, the spectrum cannot be reduced to the first-order scattering and multiple scattering must be taken into account. However, the question of its modeling arises. Physically, the same reasoning as for the first-order scattering can be iterated. The second-order scattering is illustrated in Figure 3. The computation of the successive scattering events leads thus to with x 0 = s and x n+1 = d. The components of the spectral data are then the functions g n obtained by integrating Equation (6) over the support of the object. The difficulty here is the interpretation of this variation as an operator C n (n e ) due to the complicated relations between successive scattering angles and the measured energy. This work was done for n = 2 in [41] for the 3D case and in [42] for the 2D case. The model is similar to C 1 but with a trickier integration support and a quadratic dependency on the electron density.

Extension to Polychromatic γ-ray Sources
Monochromatic sources are of course difficult to produce and extremely expensive. While X-ray tubes produce a wide spectrum with few characteristic peaks, γ-ray sources emit essentially photons with energies at their characteristic peaks. In both cases, the spectrum of an ionizing source can be decomposed into with T f lat the smooth and wide part of the spectrum and T In this work, we consider the widely used Cobalt-60 which produces two characteristic peaks at E using the integral modeling for g n .

Modeling the Impact of the Scintillation Crystals
For 60Co sources, LBC:Ce crystals demonstrates a 2.8% FWHM at 662 keV, 2.1% FWHM at 1.173 MeV and 2% FWHM at 1.332 MeV. This resolution as well as its diameter 5cm has a huge impact on our modeling and on the quality of the scattered data. Since an exact estimation of the point spread function in our model is extremely difficult due to the complexity of scattering inside the crystals, it can be relevant to simplify this model to a convolution, noted * , with a Gaussian normal distribution with standard deviation satisfying for LBC:Ce crystals In addition, the size, shape and potential capsule of the crystal, noted here Cr(d), can have an important impact on the resolution and the measurement, see Figure 5. Indeed, our models C n assume the crystal to be a point while Cr(d) is in practice a cylinder. Assuming that each element of a crystal behaves as a detector point, our models can thus be adapted by considering the transformation s Cr Figure 5. The size and shape of the crystal implies the integration along a beam of circular-arcs and not only one pair (when the crystal is assumed to be a point) anymore.
We note that depending on the shape of the specimen under study and on the desired resolution of the reconstructed image, the operator D can be negligible in comparison with the limitations of the energy resolution and the convolution with G σ . We note also that the thickness of scanning depicted in Figure 5, which affects the final resolution of the reconstructed image, corresponds to between a half and a third of the crystal depth. This is the reason why its impact may be negligible if the achieved image resolution is slightly smaller than the crystal size. This consideration shall be studied more in detail in further studies. Taking into account the stochastic nature of the emission of photons, we obtain a final model Spec LBC Co 60 (E, s, d) = Pois DSpec Co 60 (·, s, d) * G σ (E, s, d) (8) in which Pois(·) stands for the Poisson distribution and with Spec Co 60 given in Equation (7).

A General Approach for the Reconstruction of the Electron Density
In order to exploit the scattered spectrum as an imaging agent, we need to find n e from Spec LBC Co 60 in Equation (8). The main issue for solving this inverse problem is the complexity of C n , n > 1, in terms of modelling and computation time. The use of neural networks could circumvent this complexity but the lack of database for such scattered data prevents this approach. The use of standard optimization techniques such as the ROF model for total-variation regularization [43] requires however an accurate knowledge and an efficient computation of the forward model. The only suitable model for such standard approaches in our case is the first-order scattered radiation modelling, C 1 (n e ).
Following on from [41,42], the first-order scattered intensity C 1 (n e ) is less smooth than the multiple-order scattered intensity C n (n e ), n ≥ 2. Indeed, the complexity of multiple scattering tends to smooth the information about the electron density n e and thus the first-order scattering constitutes the most reliable information in the Compton spectrum. Relying upon theoretical results, this observation implies enhancing the variations of the data in order to reduce the part of multiple scattering. This can be done by considering the following problem: find n e from ∂ E Spec LBC Co 60 ≈ ∂ E C 1 (n e ) (9) in which ∂ E denotes the derivative with respect to the energy. Following the ROFmodel [43], a solution can be constructed by solving the minimization problem min n e where TV(n e ) stands for the total-variation of n e . We consider this reconstruction method in the following section.

A CT-CST Scanner
We propose combining the standard fan-beam CT scanner with CST, see Figure 6. The architecture is identical with the exception that an X-ray tube and the CCD cameras are replaced respectively by a γ-ray source (60Co) and scintillation crystals (LBC:Ce). Taking advantage of the energetical variable brought by scintillators, we measure only 12 projections. Each projection is obtained by a 30°rotation of the scanner or of the specimen.
In our setting, 40 crystals are placed on the half circle opposite the source position. For simplification of the mathematical and numerical model, every detector as well as the source is modeled as a point of no space expansion. The imaged area is a square centered in the architecture of size 150 × 150 cm 2 , see Figure 6a. The complete scanning process is depicted in Figure 6b. Furthermore, the source emits 10 6 photons per projection. We assume that the energy of the photons is measured in the energy range [355 keV, 1332 keV]. The lower bound is hereby chosen as the energy of a photon that has initial energy 1173 keV and is scattered once, changing its trajectory angle by π/2, see Equation (2). The energy range is sampled equidistantly with 128 measured levels, leading to a necessary sampling of 7.68 keV. In practice, the energy resolution of scintillation crystals is not constant; however, this information can be simply incorporated in our approach and thus we consider a constant energy resolution for the sake of simplicity.

Ballistic Data
For the simulations, we consider a simple toy object made of an annulus of aluminum with a rectangle of polyethylene inside, see Figure 7a. The details are given in the following Table 2: Table 2. Material and dimensions of the specimen and of the scintillator.

Annulus
Rectangle Crystal

Aluminium
Polyethylene LBC:Ce Surface 53 cm outer/45 cm inner radius 45 × 38 cm 2 5 × 5 cm 2 Thickness 2 cm 2 cm 5 cm The ballistic part of the spectral data, g 0 defined in Equation (1), corresponds to the standard measure in CT. Due to the scintillation crystals, the primary measurement suffers also the energy resolution of the LBC:Ce and thus can be modeled by G σ * g 0 .
Obtained by the standard TV regularization method [43], the average of the reconstructions of the attenuation map µ E at energies 1173 and 1333 keV is depicted in Figure 7b. The reconstruction suffers many artifacts which can be explained by: the energy resolution of the crystals and the very limited number of data (12 views and 40 detectors). However, this provides very important information for us. Indeed, at the energy range considered here, [355 keV, 1332 keV], the attenuation map is essentially proportional to the electron density, see [19], in which σ c (E) stands for the Compton total-cross section at energy E. Thereby, the ballistic data g 0 delivers us a first approximation for the sought-for electron density and can help us to compute a more accurate model C 1 (n e ).

Results Based on the Compton Scattered Spectrum
Due to the nature of scattering, the complete architecture has to be three dimensional. As an illustration for our general approach, we considered the much simpler fan-beam geometry CT-scan with a 2 cm thickness. This way the geometry of 3D-scattering can be approximated by a 2D-geometry. Therefore, the torus T, see Equation (5), simplifies to two circular-arcs in this setting.
The scattered data for one projection view is depicted in Figure 8c. The spectrum is here essentially composed of the first-order and second-order scattered radiation, Figure 8a photon energy measured counts photon energy measured counts In order to emphasize the general reconstruction strategy described in Section 2.4, we consider the profiles of the measured spectrum in Figure 8 for four different detector positions d 1 , d 2 , d 3 , d 4 making an angle between source and center of the circle of −20, −10, 0 and 10 degrees respectively. We observe the two characteristic peaks of 60Co and the Compton spectrum composed of first-order and multiple scattering. As demonstrated in [41,42], it is possible to exploit the smoothness properties of the different components of the spectrum. This is done simply by a discrete derivative applied on the spectral data. The derivative of the spectrum for one view and one detector is given in Figure 9. We observe that the differentiation step leads to reducing the part of multiple scattering and highlights the first-order scattering, which has the most suited structure for standard reconstruction methods. Using standard TV-regularization as for the ballistic data, we propose solving the problem (9). The result is depicted in Figure 10b. The reconstruction is satisfactory but has a relative poor accuracy. While the modeling of the crystals has a huge impact on this reconstruction, a more theoretical aspect has to be considered. Indeed, the attenuation factors A E 0 and A E 1 are the most harmful physical factors of the inverse problem as they increase dramatically the ill-posedness or ill-conditioning of the method, i.e., its instability. While this instability is controlled by the TV-regularization, it leads to a reduced quality of reconstruction. To illustrate this aspect, we consider also reducing the complete system scanner/specimen by a factor of 66%, 50% and 40%. The results are depicted in Figure 10b-d. By downsizing the object, we shorten the length of the photon travel and thus diminish the attenuation weights. As measured by the errors in Table 3, the reconstruction quality increases significantly when we reduce the size of the object. However, it also implies that the crystal size is reduced by the same factor.

Conclusions and Discussion
This work illustrates the potential of scintillation crystals in the novel Compton scattering tomography. The Compton spectrum is decomposed into multiple scattering events of order n and modeled by integral transforms. The proposed reconstruction method takes advantage of the smoothness properties of the nth-order scattered radiation by applying standard techniques combined with a differentiation step. As revealed by the simulations, the ill-posedness induced by the attenuation factors and the energy resolution of the crystals has a substantial impact on the quality and accuracy of reconstructions. To address this issue, the development of more flexible and suited reconstruction methods is the core of our future research.
A limitation of this study is the 2D aspect of the simulations. Indeed, the nature of Compton scattering is three dimensional; this is why the simulations made here will be extended to the 3D case in order to reveal the full potential of the approach. In addition, we considered that ≈10 7 photons were emitted by the γ-ray source for the total of views and were held into the slice of thickness 2 cm. In a real system, many more photons would be required in order to obtain a signal-to-noise ratio satisfactory in the data.
Furthermore, while the reconstruction of the electron density is shown possible, the resolution of reconstruction depends proportionally on the crystal size. Therefore, based on this preliminary work, the relevance and reliability of CST scanners lies upon many technological challenges: (i) the production of scintillation crystals smaller than the millimeter in order to compete with conventional CT, (ii) with the smallest FWHM possible and (iii) the production of γ-ray sources able to emit a large number of photons (10 7 photons per second) would help the method produce fast and accurate images.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript:

CeBr3
Cerium Bromide CCD Charge-coupled device CST Compton scattering tomography CT Computerized tomography LBC:Ce Lanthanum BromoChloride (La(Br x Cl 1−x )3:Ce) MSE Mean square error PSNR Peak signal-to-noise ratio