Three dimensional Compton scattering tomography

We propose a new acquisition geometry for electron density reconstruction in three dimensional X-ray Compton imaging using a monochromatic source. This leads us to a new three dimensional inverse problem where we aim to reconstruct a real valued function $f$ (the electron density) from its integrals over spindle tori. We prove injectivity of a generalized spindle torus transform on the set of smooth functions compactly supported on a hollow ball. This is obtained through the explicit inversion of a class of Volterra integral operators, whose solutions give us an expression for the harmonic coefficients of $f$. The polychromatic source case is later considered, and we prove injectivity of a new spindle interior transform, apple transform and apple interior transform on the set of smooth functions compactly supported on a hollow ball. A possible physical model is suggested for both source types. We also provide simulated density reconstructions with varying levels of added pseudo random noise and model the systematic error due to the attenuation of the incoming and scattered rays in our simulation.


Introduction
In this paper we lay the foundations for a new three dimensional imaging technique in X-ray Compton scattering tomography. Recent publications present various two dimensional scattering modalities, where a function in the plane is reconstructed from its integrals over circular arcs [2,3,4]. Three dimensional Compton tomography is also considered in the literature, where a gamma source is reconstructed from its integrals over cones with a fixed axis direction [5,6,7]. In [8], Truong and Nguyen give a history of Compton scattering tomography, from the point by point reconstruction case in earlier modalities to the circular arc transform modalities in later work. Here we present a new three dimensional scattering modality, where we aim to reconstruct the electron density (the number of electrons per unit volume) from its integrals over the surfaces of revolution of circular arcs. This work provides the theoretical basis for a new form of non invasive density determination which would be applied in fields such as fossil imaging, airport baggage screening and more generally in X-ray spectroscopic imaging. Our main goal is to show that a unique three dimensional density reconstruction is possible with knowledge of the Compton scattered intensity with our proposed acquisition geometry, and to provide an analytic expression for the density in terms of the Compton scattered data.
Compton scattering is the process which describes the inelastic scattering of photons with charged particles (usually electrons). A loss in photon energy occurs upon the collision. This is known as the Compton effect. The energy loss is dependant on the initial photon energy and the angle of scattering and is described by the following equation: Here E s is the energy of the scattered photon which had an initial energy E λ , ω is the scattering angle and E 0 ≈ 511keV is the electron rest energy. Typically Compton scattering refers to the scattering of photons in the mid energy range. That is, the scattering of X-rays and gamma rays with photon energies ranging from 1keV up to 1MeV. Forward Compton scattering is the scattering of photons with scattering angles ω ≤ π/2. Conversely, backscatter refers to the scattering of photons at angles ω > π/2.
If the photon source is monochromatic (E λ is fixed) and the detector is energy resolving (we can measure photon intensity at a given scattered energy E s ), then, in a given plane, the locus of scattering points is a circular arc connecting the source and detector points. See figure 1.
s d ω C Figure 1: A circular arc C is the locus of scattering points for a given measured energy E s for source and detector positions s and d.
A spindle torus is the surface of revolution of a circular arc. Specifically we define: to be the spindle torus, radially symmetric about the x 3 axis, with tube centre offset r ≥ 0 and tube radius √ 1 + r 2 . Let B 0,1 denote the unit ball in R 3 . Then we define the spindle S r = T r ∩ B 0,1 as the interior of the torus T r and we define the apple A r = T r \S r as the remaining exterior. See figure 2.
In three dimensions, the surface of scatterers is the surface of revolution of a circular arc C about its circle chord sd. Equivalently, the surface of scattering points is a spindle torus with an axis of revolution sd, whose tube radius and tube centre offset are determined by the distance |sd| and the scattering angle ω. In [3], Nguyen and Truong present an acquisition geometry in two dimensions for a monochromatic source (e.g. a gamma ray source) and energy resolving detector pair, where the source and detector remain at a fixed distance opposite one another and are rotated about the origin on the curve S 1 (the unit circle). Here, the dimensionality of the data is two (an energy variable and a one dimensional rotation). Taking our inspiration from Nguyen and Truong's idea, we propose a novel acquisition geometry in three dimensions for a single source and energy resolving detector pair, which are rotated opposite each other at a fixed distance about the origin on the surface S 2 (the unit sphere). Our data set is three dimensional (an energy variable and a two dimensional rotation).
As illustrated in figure 2, the forward scattered (for scattering angles ω ≤ π/2) intensity measured at the detector d for a given energy E s can be written as a weighted integral over the spindle S r (the measured energy E s determines r). With this in mind we aim to reconstruct a function supported within the unit ball from its weighted integrals over spindles. Similarly, the backscattered (ω > π/2) intensity can be given as the weighted integral over the apple A r . We also consider the exterior problem, where we aim to reconstruct a function supported on the exterior of the unit ball from its weighted integrals over apples.
In section 2 we introduce a new spindle transform for the monochromatic forward In section 3 we discuss possible approaches to the physical modelling of our problem for the case of a monochromatic and polychromatic photon source, and explain how this relates to the theory presented in section 2.
In section 4 we provide simulated density reconstructions of a test phantom via a discrete approach. We simulate data sets using the equations given in section 2 and apply our reconstruction method with varying levels of added pseudo random noise. We also simulate the added effects due to the attenuation of the incoming and scattered rays in our data and see how this systematic error effects the quality of our reconstruction.

A spindle transform
We will now parameterize the set of points on a spindle S r in terms of spherical coordinates (ρ, θ, ϕ) and give some preliminary definitions before going on to define our spindle transform later in this section.
Consider the circular arc C as illustrated in figure 3. By the cosine rule, we have: and hence: Let B 1 , 2 = {x ∈ R 3 : 1 < |x| < 2 } denote the set of points on a hollow ball with inner radius 1 and outer radius 2 , and let S 2 denote the unit sphere. Let For a function f : R 3 → R, let F : Z → R be iits polar form The arc element for the circular arc C is given in [3] as: dv = ρ 1 + r 2 1 + r 2 sin 2 ϕ dϕ (5) and the area element for the spindle S r is: We define the spindle transform S : where h = U (α)V (β).
This transform belongs to a larger class of integral transforms as we will now show.
From this we define the generalized spindle transform S w,p : where the weighting w is dependant only on r and ϕ.
We now aim to prove injectivity of the generalized spindle transform S w,p on the set of smooth functions on a hollow ball. First we give some definitions and background theory on spherical harmonics.
For integers l ≥ 0, |m| ≤ l, we define the spherical harmonics Y m l as: where, and are Legendre polynomials of degree l. The spherical harmonics Y m l form an orthonormal basis for L 2 (S 2 ), with the inner product: So we have: where δ denotes the Kroneker delta.
From [9] we have the theorem: and let: where dτ is the surface measure on the sphere. Then the series: converges uniformly absolutely on compact subsets of Z to F .
We now show how the problem of reconstructing a density f from its integrals ) can be broken down into a set of one dimensional inverse problems to solve for the harmonic coefficients f lm : ) be a curve, then: where, and P l is a Legendre polynomial degree l.
We may write h · Y m l as a linear combination of spherical harmonics of the same degree [12, page 209]: where the block diagonal entries D (l) n,m are defined: and d n,m is given, for m = 0, by: Here the D l are the irreducible blocks which form the regular representation of SO (3).
After the expansion (20) is substituted into equation (19), we see that the inserted sum is zero unless n = 0 and we have: Since the regular representation of SO(3) is unitary, we have D m,0 (h) and from the above we have: It follows that: From which we have: and the result is proven.
We now go on to show that the set of one dimensional integral equations derived above are uniquely solvable for f lm for every l ∈ N, |m| ≤ l. Given the antisymmetry of the Legendre polynomials for odd l, we see that S w,p f lm = 0 for odd l, for any ). So for the moment we focus on the reconstruction of the coefficients f lm for even l.
We now have our main theorem: and its first order partial derivative with respect to r is continuous on Proof. Given the symmetry of the Legendre polynomials P l for even l, we have: After making the substitution ρ = r sin ϕ , equation (27) becomes: a Volterra integral equation of the first kind with weakly singular kernel, where: As F lm • g is zero for ρ close to 0, given our prior assumptions regarding W and given that g ∈ C 1 ([0, 1)), we can see that K l and its first order derivative with respect to r is continuous on the support of F lm • g.
Multiplying both sides of equation (28) by 1/ √ z − r and integrating with respect to r over the interval [0, z], yields: after changing the integration order. Making the substitution r = ρ + (z − ρ)t, gives: from which we have: where, So by our assumptions that p is non zero and W is non zero on the diagonal, Q l (z, z) = 0 on the support of F lm • g.
Differentiating both sides of equation (30) with respect to z and rearranging gives: which is a Volterra type integral equation of the second kind, where: and, Given the continuity of H l on the support of F lm • g and the continuity of F lm • g on [0, 1], the Neumann series associated with equation (34) converges and we may write our solution: where the resolvent kernel, is defined by, Since the series converges, the solution is unique and we may reconstruct F lm explicitly for ρ ∈ g([0, 1)).
Here as we can only reconstruct the coefficients F lm for even l, for a general f ∈ C ∞ 0 (R 3 ) and curve p ∈ C 1 ([0, π 2 ]), the spindle data S w,p f is insufficient to recover f uniquely for |x| ∈ p([0, π 2 )). However if we consider those functions f ∈ C ∞ 0 (R 3 ) whose support lies in the upper half space x 3 > 0 we find that uniqueness is possible, as the following theorem shows: Then the even coefficients F lm for l ∈ {2k : k ∈ N} = N e , |m| ≤ l determine F uniquely and: Proof. We can write F as the sum of its even and odd coefficients: Then from our assumption, we have: for all ρ ≥ 0, 0 ≤ θ ≤ 2π and π 2 ≤ ϕ ≤ π, where: We have: The associated Legendre polynomials P m l are symmetric when l + m is even and antisymmetric otherwise. It follows that: for ρ ≥ 0, 0 ≤ θ ≤ 2π and 0 ≤ ϕ ≤ π 2 . The result follows.

A toric interior transform
In the previous section we considered the three dimensional Compton scatter tomography problem for a monochromatic source and energy sensitive detector pair. The polychromatic source case is covered here and we consider the modifications needed in our model to describe a full spectrum of initial photon energies.
In an X-ray tube electrons are accelerated by a large voltage (E m keV) towards a target material and photons are emitted. Due to conservation of energy, the emitted photons have energy no greater than E m keV. So for a given measured energy E s , where Em 1+2Em/E 0 < E s < E m , the set of scatterers is the union of spindle tori corresponding to scattering angles ω in the range 0 < ω < cos −1 1 − E 0 (Em−Es) EsEm (corresponding to energies E in the range E s < E < E m ). That is, the set of scatterers is a torus interior: where r is determined by the scattered energy measured E s . See [1] for an explanation of the two dimensional case.
Proof. We have: where, Differentiating both sides of equation (50) with respect to r and rearranging yields: where, and, Given our prior assumptions regarding w 1 , we can solve the Volterra equation of the second kind (52) uniquely for G lm (t) for 0 < t < ∞. From which we have: for t ∈ (0, 1], where w 3 (t, ϕ) = w 2 1 δ 1 t , ϕ and p(ϕ) = 1 + δ 2 1 sin 2 ϕ − δ 1 sin ϕ. By our assumptions regarding w 2 , the weighting w 3 satisfies the conditions of Theorem 2, as does the curve p. The result follows from Theorem 2.
Proof. Letf be defined as in Theorem 4. Then, after making the substitution ρ = for all r ∈ (0, ∞) and for all (α, β) ∈ S 2 , where the weighting w ≡ 1. The result follows from Theorem 4.
The advantage of using a polychromatic source (e.g. an X-ray tube) over a monochromatic source, which would most commonly be some type of gamma ray source, is that they have a significantly higher output intensity and so the data acquisition is faster.
They are also safer to handle and use and are already in use in many fields of imaging. The downside, when compared to using a monochromatic source, would be a decrease in efficiency of our reconstruction algorithm and the added differentiation step required in the inversion process. This makes the polychromatic problem more ill posed, and so small errors in our measurements would be more greatly amplified in the reconstruction. We should consider both source types for further testing to determine an optimal imaging technique.

The exterior problem
Here we consider the exterior problem, and show that a full set of apple integrals (which represent the backscattered intensity) are sufficient to reconstruct a compactly supported density on the exterior of the unit ball.

An apple transform
For backscattered photons (for scattering angles ω > π/2) the surface of scatterers is an apple A r . Refer back to figure 2. The spherical coordinates (ρ, θ, ϕ) of the points on A r can be parameterized as follows: We define the apple transform A : where h = U (α)V (β). We have the uniqueness theorem for the apple transform for functions supported on the exterior of the unit ball: Theorem 5. Let f ∈ C ∞ 0 (B 1 , 2 ∩ U ) for some 1 < 1 < 2 < ∞. Let 1 ≤ ≤ 2 and let δ = 2 −1 2 . Then Af known for 0 ≤ r ≤ δ and for all (α, β) ∈ S 2 determines f uniquely for 1 ≤ |x| ≤ .

An apple interior transform
Similar to our discussion at the start of section 2.1, if the source is polychromatic the set of backscatterers is an apple interior. Hence we define the apple interior transform and we have a theorem for its injectivity: Then AIf known for 0 ≤ r ≤ δ and for all (α, β) ∈ S 2 determines f uniquely for 1 ≤ |x| ≤ .
Proof. Letf be defined as: Then by Leibniz rule we have: So AIf known for 0 ≤ r ≤ δ and for (α, β) ∈ S 2 determines Af on the same set after differentiating. The result follows from Theorem 5.

A physical model
Here we explain how the theory presented in the previous section relates to what we measure in a practical setting.
We consider an intensity of photons scattering from a point x as illustrated in figure   4. The points s and d are the centre points of the source and detector respectively.
The intensity of photons scattered from x to d with energy E s is: Let W k (E λ ) denote the incident photon flux (number of photons per unit area per unit time), energy E λ , measured at a fixed distance D from the source. Then the incident photon intensity I 0 can be written: where t is the emission time.
The Klein-Nishina differential cross section dσ/dΩ, is defined by: where r 0 is the classical electron radius and cos ω = r/ √ 1 + r 2 . This predicts the scattering distribution for a photon off a free electron at rest. Given that the atomic electrons typically are neither free nor at rest, a correction factor is included, namely the incoherent scattering function S (q, Z). Here q = E λ hc sin (ω/2) is the momentum transferred by a photon with initial energy: scattering at an angle ω, where h is Planck's constant and c is the speed of light. As the scattering function S is dependant on the atomic number Z, we set Z = Z avg to some average atomic number as an approximation and interpolate values of S(q, Z avg ) from the tables given in [10].
The solid angle subtended by x and d may be approximated as: where A is the detector area.
The exponential terms in equation (62) account for the attenuation of the incoming and scattered rays. We approximate: In general this approximation is unrealistic. For example in medical CT, if we were scanning a relatively large (e.g. the size of someone's head) mass of organic material at a low energy (e.g. 50keV), the absorption would play a significant role and the above model would over approximate the data. However in an application where the objects are smaller (centimetres in diameter) and where we can scan at a higher energy (≈ 1MeV), the effects due to absorption would be less prevalent and Compton scattering would be the dominant interaction. For example, in airport baggage screening typical hand luggage would be a small bag containing a few low effective Z densities (e.g. nail varnish (Acetone), water, some plastic (polyethylene) etc.) and the rest may be clothes or air. Here also, as we are not worried about dosage, we can scan at high energies (e.g. using a high voltage X-ray tube or high emission energy gamma ray source). We will simulate the error due to attenuation later, in section 4, and show the effects of neglecting the attenuation in our reconstruction.

The monochromatic case
Let our density f be supported on B 1 , 2 ∩ U for some 0 < 1 < 2 < 1 and let If the source s is monochromatic (E λ remains fixed), then the forward Compton scattered intensity measured is: where p ∈ C 1 ([0, π]) is defined by p(ϕ) = 1 + δ 2 2 sin 2 ϕ − δ 2 sin ϕ, c is a constant thickness and: dσ c dΩ (r) = dσ dΩ (E s , ω) S(q, Z avg ) depends only on r. Here the variable r ∈ [0, 1] determines the scattered energy E s and (α, β) ∈ S 2 determines the source and detector position. The weighting w is given by: where ρ = 1 + r 2 sin 2 ϕ − r sin ϕ. It is left to the reader to show that the weighting w satisfies the conditions given in Theorem 2. After dividing through by the physical modelling terms in equation (68), we can invert the weighted spindle transform S w,p f as in Theorem 2 to obtain an analytic expression for f in terms of the Compton scattered intensity I.

The polychromatic case
Here we have a spectrum of incoming photon energies. See figure 5. The Compton Figure 5: A polychromatic Tungsten target spectrum with an X-ray tube voltage of V = 200keV and a 1mm Copper filter. Spectrum calculated using SpekCalc [16]. scattered intensity measured for a tube centre offset r = 1/r for a source and detector position (α, β) ∈ S 2 may be written: with the weighting: where w 2 = w as in equation (70). Here the scattering probabilty dσc dΩ and the photon flux W k are dependant on r and the integration variable t as in subsection 2.1. While the weighting w is separable as in Theorem 4, the incident photon flux W k (r , r ) = W k (E m ) = 0 for all r , where E m is the maximum spectrum energy. Hence w 1 (r , r ) ≡ 0 and w fails to meet the conditions of Theorem 4. To deal with this, let: and let w app (r , t, ϕ) = w avg (r )w 2 (t, ϕ). Then, although we sacrifice some accuracy in our forward model, w app would satisfy the conditions given in Theorem 4 and we can obtain an analytic expression for the density. The error in our approximation is bounded by: for all r ≥ 0, 0 ≤ t ≤ r . So provided the changes in the incident photon energy E s (E s is determined by t) are negligible over the range t ∈ [0, r ], the error in w app would be small.

Simulations
Here we provide reconstructions of a test phantom density at a low resolution using simulated datasets of the unweighted spindle and spindle interior transforms, and simulate noise as additive psuedo Gaussian noise. To simulate data for each transform we discretize the integrals in equations (7) and (47), and solve the least squares problem: arg min for some regularisation parameter λ > 0, where A is the discrete forward operator of the transform considered, x is the vector of density pixel values and b = Ax is the simulated transform data. We simulate perturbed data b via: where G is a pseudo random vector of samples from the standard Gaussian distibution and n is the number of entries in b. The proposed noise model has the property that b − b / b ≈ , so a noise level of is × 100% relative error.
We consider the test phantom displayed in figure 6. The unit cube is discretised into 50 × 50 × 50 pixels and a hollow ball, some stairs and a low density block with a metal  This is because the spindle interior problem is inherently more ill posed due to the extra differentiation step required in the inversion process. At a low noise level ( figure   12), while the shape of the objects can be deciphered and the ball appears to be hollow, the reconstructions are not clear, in particular the lower density block and the metal sheet start to become lost in the reconstruction. At a higher noise level (figure 13), the artefacts in the reconstruction are severe, the ball no longer appears to be hollow and the metal sheet fails to reconstruct. So, although the practical advantages of using a polychromatic source such as an X-ray tube or linear accelerator are clear (i.e. reduced data acquisition time, easy and safe to use etc.), a low noise level would need to be maintained to achieve a satisfactory image quality. The spindle transform inversion performs better when there is added noise. At a low noise level (figure 9), the size and shape of the objects is clear and the image contrast is good. At a higher noise level, the background noise in the image is amplified and the ends of metal sheet are hard to identify. We notice, in the presence of noise, that the thinner object (namely the metal sheet) is hardest to reconstruct. This is as we'd expect since smaller densities are determined by the higher frequency harmonic components which degrade faster with noise.
To simulate data for the reconstructions presented in figures 6-13 we applied the discrete operator A to the vector of test phantom pixels x and added Gaussian noise.
We now include the effects due the attenuation of the incoming and scattered rays in our simulated data and see how this effects the quality of our reconstruction. For our first example, we set the size of the scanning cube to be 20 × 20 × 20cm 3 (each pixel is 4 × 4 × 4mm 3 ) and the phantom materials are polyethylene (low density block), water (the stairs), rubber (the ball) and Aluminium (the metal sheet). We model the source as Co60, a monochromatic gamma ray source widely used in security screening applications (e.g. screening of freight shipping containers), which has an emission energy of 2824keV. We also add 1% Gaussian noise as in equation (76) to simulate random error as well as the systematic error due to attenuation effects. Our results are presented in figure 14. Here the phantom reconstruction is satisfactory with only minor artefacts appearing in the image. So if we scan a low effective Z target of a small enough size with a high energy source, the effects due to attenuation can be neglected while maintaining a satisfactory image quality. For our next example, we set the scanning cube size to be 50 × 50 × 50cm 3 (each pixel is 1 × 1 × 1cm 3 ) and the phantom materials are Teflon (low density block), PVC (stairs), polyoxymethylene (ball) and steel (the metal sheet). The source is Co60 as in our last example and again we add a further 1% Gaussian noise to the simulated data after we have accounted for the attenuative effects. Our results are presented in figure 15. Here, with a larger, higher density target, the effects due to attenuation are significant and the artefacts in the reconstruction are more severe.          We have presented a new acquisition geometry for three dimensional density reconstruction in Compton imaging with a monochromatic source and introduced a new spindle transform and a generalization of the spindle transform for the surfaces of revolution of a class of symmetric C 1 curves. The generalized spindle transform was shown to be injective on the domain of smooth functions f supported on a the intersection of a hollow ball with the upper half space x 3 > 0. In section 2 it was shown that our problem could be decomposed into a set of one dimensional inverse problems to solve for the harmonic coefficients of a given density, which we then went on to solve via the explicit inversion of a class of Volterra integral operators. Later in section 2.1 we considered the problem for a polychromatic source and introduced a new spindle interior transform and proved its injectivity on the set of smooth functions compactly supported on the intersection of a hollow ball and x 3 > 0. We also considered the exterior problem for backscattered photons with a monochromatic and polychromatic photon source, where in section 2.2 we introduced a new apple and apple interior transform. Their injectivity was proven on the set of smooth functions compactly supported on the intersection of the exterior of the unit ball and x 3 > 0. We note that in Palamodov's paper on generalized Funk transforms [14], although he provides an explicit inverse for a fairly general family of integral transforms over surfaces in three dimensional space, the spindle transform is excluded from this family of transforms.
In section 3 we discussed a possible approach to the physical modelling of our forward problem for both a monochromatic and polychromatic source. In the monochromatic source case, we found that the Compton scattered intensity resembled a weighted spindle transform (as in Theorem 2) which could be solved explicitly via repeated approximations. In the polychromatic source case, with a more accurate forward model, we found that an analytic reconstruction was not possible. So we suggested a simplified model in order to obtain an analytic expression for the density, and gave some simple error estimates for our approximation.
Test phantom reconstructions were presented in section 4 using simulated datasets from spindle and spindle interior transform data with varying levels of added pseudo random Gaussian noise. When reconstructing with spindle transform data the recon-structions were satifactory for noise levels up to 5%. We saw a harsher reduction in image quality in the presence of added noise when reconstructing from spindle interior data. This was as expected given the increased instability of the spindle interior problem.
In the latter part of section 4, we provided further reconstructions of our test phantom image in the presence of a systematic error due to the attenuative effects of the incoming and scattered rays. We found, when scanning a small low effective Z target with a high energy source, that the attenuative effects could be neglected while maintaining a good image quality. We also gave an example where the target materials were larger and higher effective Z. Here the effects due to attenuation were significant and we saw a more severe reduction in the image quality.
As the stability of the spindle transform has not yet been addressed, in future work we aim to analyse the spindle transform from a microlocal perspective to investigate whether this can shed some light on its stability. Although the injectivity of the exterior problem is covered here, simulations of an exterior density reconstruction are left for future work.