Since the advent of Computerized Tomography (CT), many imaging concepts have emerged and the need in imaging has grown. One can mention Single Photon Emission CT, Positron Emission Tomography or Cone-Beam CT for the standard system based on an ionising source. In these configurations, the energy has a very limited use but the idea of exploiting it in order to enhance the image quality, optimize the acquisition process or to compensate some limitations (such as limited angle issues) has led to various works [20, 2, 29, 11, 34, 24, 15, 14].
Computerized Tomography (CT) is a well-established and widely used technique which images an object by exploiting the properties of penetration of the x-rays. Due to the interactions of the photons with the atomic structure, the matter will resist to the propagation of the photon beam, denoted by its intensity at position and in direction , following the stationary transport equation
with the energy of the beam. The resistance of the matter is symbolized by the lineic attenuation coefficient . Solving this ordinary equation leads to the well-known Beer-Lambert law
with denotes the straight line . It describes the loss of intensities along the path to . In two dimensions, it is then standard to interpret the measurement of the intensity in a CT-scan using the Radon transform, which maps the attenuation map into its line integrals, i.e.
with and where and stand for the position of the source and of a camera. The task to recover from the data can be achieved by various techniques such as the filtered back-projection (FBP) which regularizes the inversion formula of the Radon transform, i.e.
with the -adjoint (back-projection) operator of and the Riesz potential of order .
When one focuses on the physics between the matter and the photons, four types of interactions come out: Thomson-Rayleigh scattering, photoelectric absorption, Compton scattering and pair production. In the classic range of applications of the x-rays or -rays, keV, the photoelectric absorption and the Compton scattering are the dominant phenomena which leads to a model for the lineic attenuation factor due to Stonestrom et al.  which writes
where is a factor depending on the materials and symbolizing the photoelectric absorption, the total-cross section of the Compton effect at energy and the electron density at .
The Compton effect stands for the collision of a photon with an electron. The photon transfers a part of its energy to the electron. The electron suffers then a recoil and the photon is then scattered of an (scattering) angle with the axis of propagation, see Fig. 1. The energy of the photon after scattering is expressed by the Compton formula ,
where keV represents the energy of an electron at rest.
The recent development of camera able to measure accurately the energy of incoming photons opens the way to an innovative 3D imaging concept, Compton scattering imaging (abbreviated here by CSI). Although the technology of such detectors has not yet reached the same level of maturity as the one used in conventional imaging (CT, SPECT, PET), scientists have proposed and studied in the last decades different bi-dimensional systems, called Compton scattering tomography (CST), see e.g. [23, 27, 7, 1, 3, 5, 6, 8, 12, 13, 16, 19, 25, 4, 17, 32]. It is also possible to consider interior sources, for instance via the insertion of radiotracers like in SPECT. Then considering collimated detectors, it is possible to model the scattered flux through conical Radon transforms, see for instance [26, 28]. However, in this work, we consider only systems with external sources.
In this paper we assume that the source is monochromatic, i.e. it emits photons with same energy . For sufficiently large , larger than 80keV in medical applications, the Compton effect represents a substantial part of the radiation as more of 70% of the emitted radiation is scattered within the whole body. Therefore, the variation in terms of energy due to the Compton scattering, eq. (3), will produce a polychromatic response to the monochromatic impulse due to the source. We decompose the spectrum measured at a detector with energy as follows
In this equation, represents the primary radiation which crossed the object without being subject to the Compton effect. It corresponds to the signal measured in CT, eq. (1). The functions corresponds to the photons that were measured at with incoming energy after scattering events.
The purpose of this paper is to propose a strategy to extract features of the electron density from the spectrum in eq. (4). The main idea arises from the smoothness properties of the operators modeling the scattered radiations at various orders. It appears that the first scattering encodes the richest information about the contours while the scattering of higher orders lead to much smoother data. To make our point but also for the sake of implementation, we focus on the first and second scattering only, which means and . is also discarded as it brings no information related to the energy although the primary radiation could be used as rich supplemetary information in the future.
The manuscript is organized as follows: we first recall the results from  which provides the modeling of the first scattering, , and a reconstruction strategy based on contour extraction. Then we develop an analytic modeling for the second-order scattered radiation. The study of these two parts of the spectrum under Fourier integral operators shows that is substantially smoother than on the Sobolev scale. It follows that the contour based strategy proposed in  can be straightforwardly applied to the more general case of the full spectrum. The various modelings as well as the reconstruction technique are validated using Monte-Carlo simulations.
2 First-order scattering and toric Radon transforms
In this section, the first order scattering is characterized by a toric Radon transform. Let consider an ”emission” point and a ”receiver” point . The photon beam emitted by travels in the direction within a solid angle , see Figure 1. The nature of the photon emission leads its intensity to reduce proportionally to the square of the travelled distance, the photometric dispersion, and exponentionally due to the resistance of the matter to the propagation of light, the attenuation factor or Beer-Lambert law. To represent these two physical factors, we define the following mapping of the attenuation map
Using this notation and considering a ionising source with energy and photon beam intensity , the variation of the number of photons scattered at and detected at with energy , see , can be expressed as
stands for the Klein-Nishina probability, is the classical radius of an electron. This formula describes the evolution of the first scattered radiation which is detected at a given energy and at a given detector position.
Due to the Compton formula eq. (3), the scattered energy corresponds to only one scattering angle and thus delivers a specific geometry when focusing on the first scattering. Indeed, all scattering points M responsible for a detected scattered photon at energy belongs to
the angle between two vectors. As depicted by Figure2, this set corresponds to a part of a spindle torus. More precisely, denotes the lemon part of the spindle torus for and to its apple part for , see Figure 3. Assuming now to be a point detector and integrating over the equation (5), one obtains an integral representation of the detected first scattered radiation, i.e.
As proven into , the quantities can be interpreted as a weighted toric Radon transform,
with , and the phase given by
To simplify the study but also for the sake of potential applications, we propose to fix the source . The induced modalities will then have the advantage to avoid any rotation/displacement of the source and thus enable fast acquistion time. In this case, for and , the phase simplifies to
and the first scattering radiation can be modelled by
We define as the domain of definition of . We voluntarily kept the same notation between both cases fixed or not for the sake of simplicity. One can also consider the adjoint operator in a weighted space, also called here a weighted back-projection operator
with a positive weight function, the space of the displacement of the detector w.r.t the source and the assoicated measure.
The inverse problem to solve is very difficult problem since first of all, is a non-linear operator w.r.t the electron density . Indeed the weight is an operator of the lineic attenuation factor which is also a function of , see eq. (2). Non-linear iterative techniques such as the Landweber iteration or the Kaczmarz’s method could provide satisfactory reconstruction of
but at the cost of a tremendous computation time. Deep learning techniques could be used but the lack of large dataset with real/ground truth measurement prevents any training step. A first strategy was proposed in for extracting the contours of from using the following filtered backprojection formula:
with a -smooth operator and in which
is assumed . The contours (or high variations) of are finally recovered by applying any differential operator (gradient, laplacian, wavelets, etc). The advantage of this approach is to not require specific geometry for the displacement/position of the detectors since no structure or invariances are assumed. We can then apply the approach for various architectures of 3D CSI, see Figure 4.
Eq. (6) assumes the weight to be -smooth in order to ensure the smoothness of . This is an issue for our problem since assuming to be a function with contours (singularities) will lead to a non-smooth weight . However, as pointed out in , this can be circumvented in practice by decorrelating and . One can first assume to be smooth and to be a -smooth approximation of a sharp-edges function . This approach is detailed and argumented in Section 4, while Section 5 provides simulations of the first scattering data compared with Monte-Carlo simulations as well as reconstructions of the electron density and its contours.
3 Modeling the second order scattered radiation
In this section we derive an integral representation of the part of second-order scattering in the measured spectrum. The higher orders are difficult to handle but we expect a similar approach for the higher-orders (larger than 2), this is why we will postulate the extension of the results for higher orders.
The derivation of the second scattering will use the modelling of the first scattering seen above. In this case, a measured photon is scattered twice instead of once before to be detected. We denote by M and N the first and second scattering points respectively and by and the first and second scattering angles respectively. The keyidea is to consider each first scattering point as a new polychromatic source with respect to the second scattering points and detectors.
Let us consider a detector with a detected energy . Due to the Compton formula, the scattering angles must satisfy
The boundaries 0 and 2 are here excluded since they correspond to the degenerated cases - primary () and backscattered () radiation - of the torus. In consequence, the second scattering angle will be expressed as a function of the first one:
The first angle means that a photon arriving at M with energy is scattered by an angle and has afterwards an energy . Such photons belong then to the cone with aperture , vertex M and direction , i.e. to the cone
This phase leads to a standard representation of the cone, i.e.
with the rotation matrix which maps into . Such a matrix can be computed using the Rodrigues formula.
The second angle means that a photon arriving at N with energy is scattered by an angle and is afterwards detected at with an energy . As seen in the previous section, such photons belong to the torus with fixed points and , , such that
A parametric representation of the spindle torus is given in  and writes
with the rotation matrix which maps into . Shifting the torus around , it can be expressed as a unit vector multiplied by the radius
Since the new source M emits photons with various energy in the corresponding cone, the Compton formula and the relationship implies that a photon detected at with energy and scattered at M with angle must belong to the intersection
This intersection and the principle described above is depicted in Figure 5.
To simplify the analysis, we consider the torus to be oriented in the direction , which is achieved by applying the reverse matrix . In this setting, corresponds to the third component of the normalized vector. Since we are interested in the intersection of the cone and of the torus, one gets
with the third row of the rotation matrix . Using the parametrisation of the cone, one gets for the intersection, the following radius
Since the torus is oriented ( to ), we discard the intersection between the cone and the opposite torus which corresponds to negative radii in order to fit the physics. Practically we can ignore them since they lead to detected energies outside the considered range. Also, the case is discarded as it is physically impossible and would correspond to two successive scattering occuring at the exact same location. We have now the tools to give the first result of this paper.
Considering an electron density function with compact support , a monochromatic source with energy as well as a detector both located outside . Then, the number of detected photons scattered twice arriving with an energy are given by
with the physical factors symbolized by
and the differential form of the intersection given by
The structure follows on from the physics of Compton scattering. Akin to the first scattering, the variation of photon scattered twice can be expressed as
Therefore, ignoring some physical constants and based on our development above, one can integrate to get the theoretical number of photons detected at with energy after two scattering events,
for given. First, one can compute
This leads to
Since the rotation matrices do not change the norm, they can be ignored in the computation of the norm. After some computations, one gets
This ends the proof. ∎
This Theorem delivers an integral representation for the second scattering in any configuration/architecture for the displacement/location of the source and detectors. The modeling will be validated and compared to Monte-Carlo simulations in section 5.
The next section discusses the smoothness properties of the different orders of Compton scattering within the spectrum.
4 Smoothness properties within the spectrum and contour reconstruction
This section sets our operators in the context of Fourier integral operators. For modalities involving Compton scattering, some cases were already studied in 2D  and 3D  using microlocal analysis. As mentioned above, the whole theory assumes the weight function to be smooth. In our case, this implies that the electron density must be considered in . We, however, use a first result on the continuity of our operator in .
For and fixed, the spindle tori are one-to-one with .
First, we discard the degenerate case of the spindle torus which occurs when (there is no scattering event, only primary radiation) and corresponds to the line .
The parameter representation of the spindle torus is given by
with and the rotation matrix which maps into . For and fixed, it is clear that
is one-to-one with the unit sphere. This is why each corresponds to only one pair . Now, since is defined for , the norm of the vector is a bijective function on .
Therefore, the spindle torus parametrized by is one-to-one with . ∎
Denoting by , the restriction of to one detector and a given source , then
From , we know that
with , and the appropriate Jacobian given in . Using the Cauchy-Schwarz inequality, one gets
in which stands for a smooth cut-off. Taking now the norm yields
which is well-defined as the integrand is a compactly supported smooth function. Due to Theorem 4.1 and since the discarded line passing through and is negligible regarding the Lebesgue measure, one can finally apply the mapping and gets
For our approach, the operator has two main drawbacks: (i) it is non-linear in as the physical factors, , depends on ; (ii) the properties of Fourier integral operators cannot apply if the weight function and thus are not smooth. To circumvent these difficulties, we propose for this study to decorrelate the dependence of the weight w.r.t. with the integrand itself. At this purpose, we assume that is smooth and that there exists a function such that
Thanks to the continuity of the operator shown in Theorem 4.2, the spectral data can thus be seen as an approximation of
with the operator acting on the electron density and standing for the physical factors. For our study, the edges (high-variations) of the electron density will be carried by the function and we will consider the operator instead of . Now we have the operator in a correct shape and we can recall some results from the theory on Fourier integral operators, see [35, 22].
Definition 4.3 ().
Let and be open subsets. A real valued function is called a phase function if
is positive homogeneous of degree 1 in , i.e. for all .
and do not vanish for all .
A Fourier integral operator is then defined as
where the symbol/ampliture and it satisfies:
For every compact set and for every multi-index , there is a constant such that
We write then that .
Theorem 4.4 ().
Given a FIO with symbol
then is a FIO of order and maps continuously the Sobolev spaces into for every .
Let given with an open subset of . Then the operator is a continuous Fourier integral operator of order .
The proof here is inspired from . Using the Fourier representation of the Dirac distribution, one gets
Defining the phase
one needs to prove that and do not vanish for all . It is clear that
The second term is never zero, so does not vanish. Also, setting and , one can prove that
The gradient of the phase is thus the sum of two vectors each one collinear with and respectively. Important cases here show up, first must be different from , which means . This is naturally excluded by taking the source outside . Then, comes the case , which means and are collinear. It corresponds to the spindle torus degenerating into the straight line to (no scattering only primary radiation). This case is also excluded since strictly or equivalently . Since and cannot be collinear, the gradient is never zero and therefore, is a FIO.
Due to the smoothness properties of the electron density, and is thus a symbol of order 0. Given Theorem 4.4 and that and do not vanish, is thus a FIO of order
and therefore maps continuously into where denotes the domain of definition of . ∎
Let given with an open subset of . Then, the scattered radiation of order 2 can be understood as a continuous mapping between and if the phase function
In order to write the second scattering data as a FIO, one needs to find the corresponding phase function of the intersection between the cone and the spindle torus. The phase is delivered by our condition on the scattering angles and the energy derived from the Compton formula eq. (7). It yields
inherits essentially the properties of phase of and of as is monotone w.r.t the energy. Although it can be empirically checked, it is unclear when analytically, this is why this is assumed here. With this symbolism, can be understood as an approximation of the operator
The case is discarded as it corresponds for the tori to a singularity that shows up in the definition of the phase (it was the case for the first scattering). Physically, this case cannot occur as it would imply that the same local point is responsible for two successive scattering events which is impossible. To interpret this as a FIO, we embed the variable into , with denoting a small neighborhood around . Taking , it yields
with and the appropriate embedded version of and . Given Theorem 4.4 and assuming the properties of phase of , is a FIO of order
which ends the proof. ∎
The assumption that is in fact almost always satisfied. Indeed, the computation of is a linear combination of three vectors , and . Since the degenerated cases for the cone, and collinear, and for the spindle torus, and collinear, are excluded, there are two possibilities: either these three vectors are linearly independant, then inherits the property of phase of (), or they belong to the same plane. In the latter case, the analysis is unclear. The desired property can however be checked empirically.
By analogy, the scattering of order () will rely on the relation
with the scattering angle. Unfortunately, the geometry of the scattering events becomes harder to model or even to implement. However, we expect the principle we developed for the second scattering to extend to higher orders since each additional scattering will add a new layer of intermediary sources. This extension is expressed in the following statement.
When focusing on the Compton scattering, the measured spectrum can be decomposed into the different scattered radiation of order,
with apprximating the operator which maps continuously
Intuitively, the more scattering events, the smoother the corresponding part within the spectrum is.
Assuming the attenuation factors to be known and to be a sharp version of the smooth electron density , the reconstruction operator with the backprojection operator
maps the spectrum onto up to an error .
The result follows on from the smoothness properties given above and from . The reconstruction operator recovers up to a -smooth function and acts as a FIO of order 1 on the second scattering data