Three dimensional Compton scattering tomography
TThree dimensional Compton scatteringtomography
By James Webber and William Lionheart
Abstract
We propose a new acquisition geometry for electron density reconstructionin three dimensional X-ray Compton imaging using a monochromatic source.This leads us to a new three dimensional inverse problem where we aim toreconstruct a real valued function f (the electron density) from its integralsover spindle tori. We prove injectivity of a generalized spindle torus transformon the set of smooth functions compactly supported on a hollow ball. This isobtained through the explicit inversion of a class of Volterra integral operators,whose solutions give us an expression for the harmonic coefficients of f . Thepolychromatic source case is later considered, and we prove injectivity of a newspindle interior transform, apple transform and apple interior transform on theset of smooth functions compactly supported on a hollow ball.A possible physical model is suggested for both source types. We also providesimulated density reconstructions with varying levels of added pseudo randomnoise and model the systematic error due to the attenuation of the incoming andscattered rays in our simulation. a r X i v : . [ m a t h . F A ] A p r Introduction
In this paper we lay the foundations for a new three dimensional imaging techniquein X-ray Compton scattering tomography. Recent publications present various twodimensional scattering modalities, where a function in the plane is reconstructed fromits integrals over circular arcs [2, 3, 4]. Three dimensional Compton tomography isalso considered in the literature, where a gamma source is reconstructed from its in-tegrals over cones with a fixed axis direction [5, 6, 7]. In [8], Truong and Nguyen givea history of Compton scattering tomography, from the point by point reconstructioncase in earlier modalities to the circular arc transform modalities in later work. Herewe present a new three dimensional scattering modality, where we aim to reconstructthe electron density (the number of electrons per unit volume) from its integrals overthe surfaces of revolution of circular arcs. This work provides the theoretical basisfor a new form of non invasive density determination which would be applied in fieldssuch as fossil imaging, airport baggage screening and more generally in X-ray spec-troscopic imaging. Our main goal is to show that a unique three dimensional densityreconstruction is possible with knowledge of the Compton scattered intensity with ourproposed acquisition geometry, and to provide an analytic expression for the densityin terms of the Compton scattered data.Compton scattering is the process which describes the inelastic scattering of pho-tons with charged particles (usually electrons). A loss in photon energy occurs uponthe collision. This is known as the Compton effect. The energy loss is dependant onthe initial photon energy and the angle of scattering and is described by the followingequation: E s = E λ E λ /E ) (1 − cos ω ) . (1)Here E s is the energy of the scattered photon which had an initial energy E λ , ω isthe scattering angle and E ≈ ω ≤ π/
2. Conversely, backscatter refers to the scattering of photons at angles ω > π/ E λ is fixed) and the detector is energy2esolving (we can measure photon intensity at a given scattered energy E s ), then, in agiven plane, the locus of scattering points is a circular arc connecting the source anddetector 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: T r = (cid:40) ( x , x , x ) ∈ R : (cid:18) r − (cid:113) x + x (cid:19) + x = 1 + r (cid:41) (2)to be the spindle torus, radially symmetric about the x axis, with tube centre offset r ≥ √ r . Let B , denote the unit ball in R . Then we definethe spindle S r = T r ∩ B , 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 circulararc C about its circle chord sd . Equivalently, the surface of scattering points is aspindle torus with an axis of revolution sd , whose tube radius and tube centre offsetare determined by the distance | sd | and the scattering angle ω . In [3], Nguyen andTruong 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 anddetector remain at a fixed distance opposite one another and are rotated about theorigin on the curve S (the unit circle). Here, the dimensionality of the data is two (anenergy variable and a one dimensional rotation). Taking our inspiration from Nguyenand Truong’s idea, we propose a novel acquisition geometry in three dimensions for asingle source and energy resolving detector pair, which are rotated opposite each otherat a fixed distance about the origin on the surface S (the unit sphere). Our data setis three dimensional (an energy variable and a two dimensional rotation).As illustrated in figure 2, the forward scattered (for scattering angles ω ≤ π/ d for a given energy E s can be written as a weighted3 sω Hx x S r A r objectscanning region rR Figure 2: A spindle torus slice with height H , tube radius R and tube centre offset r has tips at source and detector points s and d . The spindle S and the apple A arehighlighted by a solid green line and a dashed blue line respectively. The sphericalscanning region highlighted has unit radius.integral over the spindle S r (the measured energy E s determines r ). With this in mindwe aim to reconstruct a function supported within the unit ball from its weightedintegrals over spindles. Similarly, the backscattered ( ω > π/
2) intensity can be givenas 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 fromits weighted integrals over apples.In section 2 we introduce a new spindle transform for the monochromatic forward4catter problem and introduce a generalization of the spindle transform which gives theintegrals of a function over the surfaces of revolution of a particular class of symmetriccurves. We prove the injectivity of the generalized spindle transform on the domainof smooth functions compactly supported on the intersection of a hollow ball with theupper half space x >
0. We show that our problem can be decomposed as a set of onedimensional inverse problems, which we then solve to provide an explicit expressionfor the harmonic coefficients of f . In section 2.1 we introduce a new toric interiortransform for the polychromatic forward scatter problem and prove its injectivity onthe domain of smooth functions compactly supported on the intersection of a hollowball and x >
0. A new apple and apple interior transform are also introduced insection 2.2 for the monochromatic and polychromatic backscatter problem. Theirinjectivity is proven on the domain of smooth functions compactly supported on theintersection of the exterior of the unit ball and x > We will now parameterize the set of points on a spindle S r in terms of sphericalcoordinates ( ρ, θ, ϕ ) and give some preliminary definitions before going on to defineour spindle transform later in this section.Consider the circular arc C as illustrated in figure 3. By the cosine rule, we have: R = ρ + r + 2 rρ sin ϕ (3)and hence: ρ = (cid:113) r sin ϕ + 1 − r sin ϕ. (4)5 sω x x ϕ C objectscanning region rR ρ Figure 3: A circular arc C connecting source and detector points s and d . The circularscanning region highlighted has unit radius.Let B (cid:15) ,(cid:15) = { x ∈ R : (cid:15) < | x | < (cid:15) } denote the set of points on a hollow ballwith inner radius (cid:15) and outer radius (cid:15) , and let S denote the unit sphere. Let Z = R + × S . For a function f : R → R , let F : Z → R be iits polar form F ( ρ, θ, ϕ ) = f ( ρ sin ϕ cos θ, ρ sin ϕ sin θ, ρ cos ϕ ). We parameterize h ∈ SO(3) in termsof Euler angles α and β , h = U ( α ) V ( β ), where U and V are rotations about the x and x axis respectively. We define an action of the rotation group SO(3) on the setof real valued functions f on R in the natural way, and define an action on the polarform as ( h · F ) = ( h · f ).The arc element for the circular arc C is given in [3] as:d v = ρ (cid:115) r r sin ϕ d ϕ (5)6nd the area element for the spindle S r is:d A = ρ sin ϕ d v d θ. (6)We define the spindle transform S : C ∞ ( R ) → C ∞ ( Z ) as: S f ( r, α, β ) = (cid:90) π (cid:90) π ρ sin ϕ (cid:115) r r sin ϕ ( h · F ) ( ρ, θ, ϕ ) | ρ = √ r sin ϕ +1 − r sin ϕ d ϕ d θ (7)where h = U ( α ) V ( β ).This transform belongs to a larger class of integral transforms as we will now show.Let p ∈ C ([0 , π ]) be a curve parameterized by a colatitude ϕ ∈ [0 , π ]. Then we definethe class of symmetric curves ρ ∈ C ([0 , × [0 , π ]) by: ρ ( r, ϕ ) = p (sin − ( r sin ϕ )) . (8)From this we define the generalized spindle transform S w,p : C ∞ ( R ) → C ∞ ([0 , × S )as: S w,p f ( r, α, β ) = (cid:90) π (cid:90) π w ( r, ϕ ) ρ sin ϕ (cid:115) ρ + (cid:18) d ρ d ϕ (cid:19) ( h · F ) ( ρ, θ, ϕ ) | ρ = ρ ( r,ϕ ) d ϕ d θ. (9)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 setof smooth functions on a hollow ball. First we give some definitions and backgroundtheory on spherical harmonics.For integers l ≥ , | m | ≤ l , we define the spherical harmonics Y ml as: Y ml ( θ, ϕ ) = ( − m (cid:115) (2 l + 1)( l − m )!4 π ( l + m )! P ml (cos ϕ ) e i mθ (10)where, P ml ( x ) = ( − m (1 − x ) m/ d m d x m P l ( x ) (11)and P l ( x ) = 12 l l (cid:88) k =0 (cid:18) lk (cid:19) ( x − l − k ( x + 1) k , (12)are Legendre polynomials of degree l . The spherical harmonics Y ml form an orthonor-mal basis for L ( S ), with the inner product: (cid:104) f, g (cid:105) = (cid:90) S f ¯ g d τ = (cid:90) π (cid:90) π f ( θ, ϕ )¯ g ( θ, ϕ ) sin ϕ d θ d ϕ. (13)7o we have: (cid:90) S Y ml ¯ Y m (cid:48) l (cid:48) d τ = (cid:90) π (cid:90) π Y ml ( θ, ϕ ) ¯ Y m (cid:48) l (cid:48) ( θ, ϕ ) sin ϕ d θ d ϕ = δ ll (cid:48) δ mm (cid:48) , (14)where δ denotes the Kroneker delta.From [9] we have the theorem: Theorem 1.
Let f ∈ C ∞ ( R ) and let: F lm = (cid:90) S F ¯ Y ml d τ, (15) where d τ is the surface measure on the sphere. Then the series: F N = (cid:88) ≤ l ≤ N (cid:88) | m |≤ l F lm Y ml (16) converges uniformly absolutely on compact subsets of Z to F . We now show how the problem of reconstructing a density f from its integrals S w,p f for some curve p ∈ C ([0 , π ]) can be broken down into a set of one dimensionalinverse problems to solve for the harmonic coefficients f lm : Lemma 1.
Let f ∈ C ∞ ( R ) and let p ∈ C ([0 , π ]) be a curve, then: S w,p f lm ( r ) = 2 π (cid:90) π w ( r, ϕ ) ρ sin ϕ (cid:115) ρ + (cid:18) d ρ d ϕ (cid:19) F lm ( ρ ) P l (cos ϕ ) | ρ = ρ ( r,ϕ ) d ϕ (17) where, S w,p f lm = (cid:90) S ( S w,p f ) ¯ Y ml d τ (18) and P l is a Legendre polynomial degree l .Proof. Let τ ∈ S and let h ∈ SO(3) act on a spherical harmonic Y ml by ( h · Y ml )( τ ) = Y ml ( hτ ). Then we have: S w,p f ( r, α, β ) = (cid:90) π (cid:90) π w ( r, ϕ ) ρ sin ϕ (cid:115) ρ + (cid:18) d ρ d ϕ (cid:19) ( h · F ) ( ρ, θ, ϕ ) | ρ = ρ ( r,ϕ ) d ϕ d θ = (cid:88) l ∈ N (cid:88) | m |≤ l (cid:90) π (cid:90) π w ( r, ϕ ) ρ sin ϕ (cid:115) ρ + (cid:18) d ρ d ϕ (cid:19) F lm ( ρ )( h · Y ml )( θ, ϕ )d ϕ d θ (19)where h = U ( α ) V ( β ). 8e may write h · Y ml as a linear combination of spherical harmonics of the samedegree [12, page 209]: ( h · Y ml )( τ ) = (cid:88) | n | Let p ∈ C ([0 , π ]) be such that p ( ϕ ) (cid:54) = 0 for all ϕ ∈ [0 , π ] , and let g ∈ C ([0 , be defined as g ( x ) = p (sin − x ) . Let A = { ρ ∈ [0 , 1) : g ( ρ ) ∈ [0 , (cid:15) ] ∪ [ (cid:15) , (cid:15) ] } and let w : [0 , × [0 , π ] → R be a weighting such that W ( r, ρ ) = w (cid:0) r, sin − ( ρ/r ) (cid:1) and its first order partial derivative with respect to r is continuous on { ( r, ρ ) : r ∈ A, min( A ) ≤ ρ ≤ r } and W ( ρ, ρ ) (cid:54) = 0 on A . Let f ∈ C ∞ ( B ,(cid:15) ∪ B (cid:15) ,(cid:15) ) , where (cid:15) < p (0) < (cid:15) . Then S w,p f determines the harmonic coefficients F lm uniquely for ρ ∈ g ([0 , for all l ∈ { k : k ∈ N } , | m | ≤ l .Proof. Given the symmetry of the Legendre polynomials P l for even l , we have:14 π S w,p f lm ( r ) = (cid:90) π w ( r, ϕ ) ρ sin ϕ (cid:115) ρ + (cid:18) d ρ d ϕ (cid:19) F lm ( ρ ) P l (cos ϕ ) | ρ = ρ ( r,ϕ ) d ϕ = (cid:90) π w ( r, ϕ ) g ( ρ ) sin ϕ (cid:112) g ( ρ ) + r cos ϕg (cid:48) ( ρ ) F lm ( g ( ρ )) P l (cos ϕ ) | ρ = r sin ϕ d ϕ. (27)After making the substitution ρ = r sin ϕ , equation (27) becomes:14 π S w,p f lm ( r ) = (cid:90) r F lm ( g ( ρ )) K l ( r, ρ ) √ r − ρ d ρ (28)a Volterra integral equation of the first kind with weakly singular kernel, where: K l ( r, ρ ) = W ( r, ρ ) ρg ( ρ ) (cid:112) g ( ρ ) + ( r − ρ ) g (cid:48) ( ρ ) r √ r + ρ P l (cid:32)(cid:114) − ρ r (cid:33) . (29)As F lm ◦ g is zero for ρ close to 0, given our prior assumptions regarding W and giventhat g ∈ C ([0 , 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 respectto r over the interval [0 , z ], yields:14 π (cid:90) z S w,p f lm ( r ) √ z − r d r = (cid:90) z F lm ( g ( ρ )) (cid:20)(cid:90) zρ K l ( r, ρ ) √ z − r √ r − ρ d r (cid:21) d ρ (30)10fter changing the integration order. Making the substitution r = ρ + ( z − ρ ) t , gives: Q l ( z, ρ ) = (cid:90) zρ K l ( r, ρ ) √ z − r √ r − ρ d r = (cid:90) K l ( ρ + ( z − ρ ) t, ρ ) √ t √ − t d t (31)from which we have: Q l ( z, z ) = πK l ( z, z ) = W ( z, z ) πc l g ( z ) √ z (32)where, c l = P l (0) = ( − ( l/ l (cid:18) ll/ (cid:19) . (33)So by our assumptions that p is non zero and W is non zero on the diagonal, Q l ( z, z ) (cid:54) =0 on the support of F lm ◦ g .Differentiating both sides of equation (30) with respect to z and rearranging gives: v lm ( z ) = (cid:90) z F lm ( g ( ρ )) H l ( z, ρ )d ρ + F lm ( g ( z )) (34)which is a Volterra type integral equation of the second kind, where: v lm ( z ) = 14 π K l ( z, z ) dd z (cid:90) z S w,p f lm ( r ) √ z − r d r (35)and, H l ( z, ρ ) = 1 πK l ( z, z ) dd z Q l ( z, ρ ) . (36)Given the continuity of H l on the support of F lm ◦ g and the continuity of F lm ◦ g on[0 , F lm ( g ( z )) = (cid:90) z R l ( z, ρ ) v lm ( ρ )d ρ + v lm ( z ) (37)where the resolvent kernel, R l ( z, ρ ) = ∞ (cid:88) i =1 H l,i ( z, ρ ) (38)is defined by, H l, ( z, ρ ) = H l ( z, ρ ) , H l,i ( z, ρ ) = (cid:90) zρ H l ( z, x ) H l,i − ( x, ρ )d x ∀ i ≥ . (39)Since the series converges, the solution is unique and we may reconstruct F lm explicitlyfor ρ ∈ g ([0 , F lm for even l , for a general f ∈ C ∞ ( R ) and curve p ∈ C ([0 , π ]), the spindle data S w,p f is insufficient to recover f uniquely for | x | ∈ p ([0 , π )). However if we consider those functions f ∈ C ∞ ( R ) whosesupport lies in the upper half space x > Theorem 3. Let U = { ( x , x , x ) ∈ R : x > } and let f ∈ C ∞ ( U ) . Then the evencoefficients F lm for l ∈ { k : k ∈ N } = N e , | m | ≤ l determine F uniquely and: F ( ρ, θ, ϕ ) = 2 (cid:88) l ∈ N e (cid:88) | m |≤ l F lm ( ρ ) Y ml ( θ, ϕ ) (40) for ρ ≥ , ≤ θ ≤ π and ≤ ϕ ≤ π .Proof. We can write F as the sum of its even and odd coefficients: F ( ρ, θ, ϕ ) = F o ( ρ, θ, ϕ )+ F e ( ρ, θ, ϕ ) = (cid:88) odd l ∈ N (cid:88) | m |≤ l F lm ( ρ ) Y ml ( θ, ϕ )+ (cid:88) even l ∈ N (cid:88) | m |≤ l F lm ( ρ ) Y ml ( θ, ϕ ) . (41)Then from our assumption, we have: − F e ( ρ, θ, ϕ ) = F o ( ρ, θ, ϕ ) = (cid:88) odd l ∈ N (cid:88) | m |≤ l F lm ( ρ ) Y ml ( θ, ϕ )= (cid:88) m ∈ Z (cid:88) odd l ≥| m | c ( l, m ) F lm ( ρ ) P ml (cos ϕ ) e imθ = (cid:88) m ∈ Z F mo ( ρ, ϕ ) e imθ (42)for all ρ ≥ 0, 0 ≤ θ ≤ π and π ≤ ϕ ≤ π , where: c ( l, m ) = ( − m (cid:115) (2 l + 1)( l − m )!4 π ( l + m )! . (43)We have: − F me ( ρ, ϕ ) = − π (cid:90) π F e ( ρ, θ, ϕ ) e − imθ d θ = F mo ( ρ, ϕ ) for ρ ≥ , π ≤ ϕ ≤ π. (44)The associated Legendre polynomials P ml are symmetric when l + m is even andantisymmetric otherwise. It follows that: F o ( ρ, θ, ϕ ) = (cid:88) m ∈ Z ( − m +1 F mo ( ρ, π − ϕ ) e imθ = (cid:88) m ∈ Z ( − m F me ( ρ, π − ϕ ) e imθ = (cid:88) m ∈ Z F me ( ρ, ϕ ) e imθ (45)12or ρ ≥ 0, 0 ≤ θ ≤ π and 0 ≤ ϕ ≤ π . The result follows. Corollary 1. Let f ∈ C ∞ ( B (cid:15) ,(cid:15) ∩ U ) for some 0 < (cid:15) < (cid:15) < 1. Let (cid:15) ≤ (cid:15) ≤ (cid:15) and let δ = − (cid:15) (cid:15) . Then S f known for 0 ≤ r ≤ δ and for all ( α, β ) ∈ S determines f uniquely for (cid:15) ≤ | x | ≤ Proof. Let p ( ϕ ) = (cid:112) δ sin ϕ − δ sin ϕ and let w ≡ 1. Then S f ( δr, α, β ) = S w,p f ( r, α, β ) for r ∈ [0 , α, β ) ∈ S . The result follows from Theorems 2 and3. In the previous section we considered the three dimensional Compton scatter tomog-raphy problem for a monochromatic source and energy sensitive detector pair. Thepolychromatic source case is covered here and we consider the modifications neededin 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 atarget material and photons are emitted. Due to conservation of energy, the emittedphotons have energy no greater than E m keV. So for a given measured energy E s , where E m E m /E < E s < E m , the set of scatterers is the union of spindle tori correspondingto scattering angles ω in the range 0 < ω < cos − (cid:16) − E ( E m − E s ) E s E m (cid:17) (corresponding toenergies E in the range E s < E < E m ). That is, the set of scatterers is a torus interior: I r = (cid:40) ( x , x , x ) ∈ R : (cid:18) r − (cid:113) x + x (cid:19) + x < r (cid:41) (46)where r is determined by the scattered energy measured E s . See [1] for an explanationof the two dimensional case.With this in mind we define the spindle interior transform I : C ∞ ( R ) → C ∞ ( Z )as: I f ( r, α, β ) = (cid:90) π (cid:90) π (cid:90) √ r sin ϕ +1 − r sin ϕ ρ sin ϕ ( h · F )( ρ, θ, ϕ )d ρ d ϕ d θ. (47)and we have the uniqueness theorem for a weighted spindle interior transform:13 heorem 4. Let f ∈ C ∞ ( B (cid:15) ,(cid:15) ∩ U ) for some < (cid:15) < (cid:15) < and and let δ = − (cid:15) (cid:15) and δ = − (cid:15) (cid:15) . Let ˜ f be defined as: (cid:113) (1 −| x | ) | x | + 1 ˜ f ( x ) = − | x | | x | − (1 −| x | ) | x | (cid:113) (1 −| x | ) | x | + 1 f ( x ) . (48) Define the weighted interior transform: I w f ( r (cid:48) , α, β ) = (cid:90) r (cid:48) (cid:90) π (cid:90) π w ( r (cid:48) , t, ϕ ) ρ sin ϕt (cid:113) sin ϕt + 1 ( h · ˜ F ) ( ρ, θ, ϕ ) | ρ = (cid:113) sin2 ϕt +1 − sin ϕt d ϕ d θ d t, (49) where h = U ( α ) V ( β ) and r (cid:48) = 1 /r . Let us suppose that the weighting w satisfies thefollowing:1. w can be decomposed as w ( r (cid:48) , t, ϕ ) = w ( r (cid:48) , t ) w ( t, ϕ ) .2. W ( t, ρ ) = w (cid:0) /t, sin − ( ρ/t ) (cid:1) and its first order partial derivative with respectto t is continuous on the triangle T = { ( t, ρ ) ∈ R : (cid:15) ≤ t ≤ (cid:15) , (cid:15) ≤ ρ ≤ t } and W ( ρ, ρ ) (cid:54) = 0 on [ (cid:15) , (cid:15) ] .3. w ( r (cid:48) , t ) and its first order partial derivatives are bounded on ≤ r (cid:48) < ∞ and w ( r (cid:48) , r (cid:48) ) (cid:54) = 0 for ≤ r (cid:48) < ∞ .Then I w f known for all r (cid:48) ∈ [0 , ∞ ) and for all ( α, β ) ∈ S determines f uniquely for < | x | < .Proof. We have: I w f ( r (cid:48) , α, β ) = (cid:90) r (cid:48) w ( r (cid:48) , t ) G ( t, α, β )d t, (50)where, G ( t, α, β ) = (cid:90) π (cid:90) π w ( t, ϕ ) ρ sin ϕt (cid:113) sin ϕt + 1 ( h · ˜ F ) ( ρ, θ, ϕ ) | ρ = (cid:113) sin2 ϕt +1 − sin ϕt d ϕ d θ. (51)Differentiating both sides of equation (50) with respect to r (cid:48) and rearranging yields: g ( r (cid:48) , α, β ) = (cid:90) r (cid:48) L ( r (cid:48) , t ) G ( t, α, β )d t + G ( r (cid:48) , α, β ) (52)where, g ( r (cid:48) , α, β ) = − dd r (cid:48) I w f ( r (cid:48) , α, β ) w ( r (cid:48) , r (cid:48) ) (53)14nd, L ( r (cid:48) , t ) = dd r (cid:48) w ( r (cid:48) , t ) w ( r (cid:48) , r (cid:48) ) . (54)Given our prior assumptions regarding w , we can solve the Volterra equation of thesecond kind (52) uniquely for G lm ( t ) for 0 < t < ∞ . From which we have: G (cid:18) δ t , α, β (cid:19) = (cid:90) π (cid:90) π w (cid:18) δ t , ϕ (cid:19) δ tρ sin ϕ (cid:112) δ t sin ϕ + 1 ( h · ˜ F ) ( ρ, θ, ϕ ) | ρ = √ δ t sin ϕ +1 − δ t sin ϕ d ϕ d θ = δ t (cid:112) δ t S w ,p ˜ f ( t, α, β ) (55)for t ∈ (0 , w ( t, ϕ ) = w (cid:16) δ t , ϕ (cid:17) and p ( ϕ ) = (cid:112) δ sin ϕ − δ sin ϕ . Byour assumptions regarding w , the weighting w satisfies the conditions of Theorem 2,as does the curve p . The result follows from Theorem 2. Corollary 2. Let f ∈ C ∞ ( B (cid:15) ,(cid:15) ∩ U ) for some 0 < (cid:15) < (cid:15) < 1. Then I f as definedin equation (47) known for all r ∈ [0 , ∞ ) and for all ( α, β ) ∈ S determines f uniquelyfor 0 < | x | < Proof. Let ˜ f be defined as in Theorem 4. Then, after making the substitution ρ = (cid:113) sin ϕt + 1 − sin ϕt in equation (47), we have: I f ( r, α, β ) = (cid:90) r (cid:90) π (cid:90) π ρ sin ϕt (cid:113) sin ϕt + 1 ( h · ˜ F ) ( ρ, θ, ϕ ) | ρ = (cid:113) sin2 ϕt +1 − sin ϕt d ϕ d θ d t = I w f (1 /r, α, β ) (56)for all r ∈ (0 , ∞ ) and for all ( α, β ) ∈ S , where the weighting w ≡ 1. The resultfollows from Theorem 4.The advantage of using a polychromatic source (e.g. an X-ray tube) over a monochro-matic source, which would most commonly be some type of gamma ray source, is thatthey 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 imag-ing. The downside, when compared to using a monochromatic source, would be adecrease in efficiency of our reconstruction algorithm and the added differentiationstep required in the inversion process. This makes the polychromatic problem more illposed, and so small errors in our measurements would be more greatly amplified in the15econstruction. We should consider both source types for further testing to determinean optimal imaging technique. 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 compactlysupported density on the exterior of the unit ball. For backscattered photons (for scattering angles ω > π/ 2) the surface of scatterers isan apple A r . Refer back to figure 2. The spherical coordinates ( ρ, θ, ϕ ) of the pointson A r can be parameterized as follows: ρ = (cid:113) r sin ϕ + 1 + r sin ϕ, ≤ θ ≤ π, ≤ ϕ ≤ π. (57)We define the apple transform A : C ∞ ( R ) → C ∞ ( Z ) as: A f ( r, α, β ) = (cid:90) π (cid:90) π ρ sin ϕ (cid:115) r r sin ϕ ( h · F )( ρ, θ, ϕ ) | ρ = √ r sin ϕ +1+ r sin ϕ d ϕ d θ, (58)where h = U ( α ) V ( β ). We have the uniqueness theorem for the apple transform forfunctions supported on the exterior of the unit ball: Theorem 5. Let f ∈ C ∞ ( B (cid:15) ,(cid:15) ∩ U ) for some < (cid:15) < (cid:15) < ∞ . Let (cid:15) ≤ (cid:15) ≤ (cid:15) and let δ = (cid:15) − (cid:15) . Then A f known for ≤ r ≤ δ and for all ( α, β ) ∈ S determines f uniquely for ≤ | x | ≤ (cid:15) .Proof. Let p ( ϕ ) = (cid:112) δ sin ϕ + δ sin ϕ and let w ≡ 1. Then A f ( δr, α, β ) = S w,p f ( r, α, β ) for r ∈ [0 , α, β ) ∈ S . The result follows from Theorems 2 and3. Similar to our discussion at the start of section 2.1, if the source is polychromatic theset of backscatterers is an apple interior. Hence we define the apple interior transform16 I : C ∞ ( R ) → C ∞ ( Z ): AI f ( r, α, β ) = (cid:90) π (cid:90) π (cid:90) √ r sin ϕ +1+ r sin ϕ ρ sin ϕ ( h · F )( ρ, θ, ϕ )d ρ d ϕ d θ. (59)and we have a theorem for its injectivity: Theorem 6. Let f ∈ C ∞ ( B (cid:15) ,(cid:15) ∩ U ) for some < (cid:15) < (cid:15) < ∞ . Let (cid:15) ≤ (cid:15) ≤ (cid:15) andlet δ = (cid:15) − (cid:15) . Then AI f known for ≤ r ≤ δ and for all ( α, β ) ∈ S determines f uniquely for ≤ | x | ≤ (cid:15) .Proof. Let ˜ f be defined as:1 (cid:113) ( | x | − | x | + 1 ˜ f ( x ) = ( | x | − | x | (cid:113) ( | x | − | x | + 1 + | x | − | x | f ( x ) . (60)Then by Leibniz rule we have: r dd r AI f ( r, α, β ) = (cid:90) π (cid:90) π ρ sin ϕ (cid:112) r sin ϕ ( h · ˜ F )( ρ, θ, ϕ ) | ρ = √ r sin ϕ +1+ r sin ϕ d ϕ d θ = 1 √ r A ˜ f ( r, α, β ) (61)So AI f known for 0 ≤ r ≤ δ and for ( α, β ) ∈ S determines A ˜ f on the same set afterdifferentiating. The result follows from Theorem 5. Here we explain how the theory presented in the previous section relates to what wemeasure in a practical setting.We consider an intensity of photons scattering from a point x as illustrated in figure4. 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: I ( E s ) = I ( E λ ) e − (cid:82) Lsx µ ( E λ ,Z ) f ( x ) d V × d σ dΩ ( E s , ω ) e − (cid:82) Lxd µ ( E s ,Z ) S ( q, Z ) dΩ . (62)where Z denotes the atomic number, ω the scattering angle and L sx and L xd are theline segments connecting s to x and x to d respectively. f ( x ) denotes the electron17 sϕ E λ E s x CrR ρ Figure 4: A scattering event occurs on the circular arc C with initial photon energy E λ from the source s to the detector d with energy E s .density (number of electrons per unit volume) at the scattering point x and d V is thevolume measure. f is the quantity to be reconstructed.Let W k ( E λ ) denote the incident photon flux (number of photons per unit area perunit time), energy E λ , measured at a fixed distance D from the source. Then theincident photon intensity I can be written: I ( E λ ) = tD W k ( E λ )( ρ cos ϕ + 1) + ρ sin ϕ , (63)where t is the emission time.The Klein-Nishina differential cross section d σ/ dΩ, is defined by:d σ dΩ ( E s , ω ) = r (cid:18) E s E λ (cid:19) (cid:18) E s E λ + E λ E s − ω (cid:19) , (64)where r is the classical electron radius and cos ω = r/ √ r . This predicts thescattering distribution for a photon off a free electron at rest. Given that the atomicelectrons typically are neither free nor at rest, a correction factor is included, namelythe incoherent scattering function S ( q, Z ). Here q = E λ hc sin ( ω/ 2) is the momentumtransferred by a photon with initial energy: E λ = E s − ( E s /E ) (1 − cos ω ) (65)scattering at an angle ω , where h is Planck’s constant and c is the speed of light. Asthe scattering function S is dependant on the atomic number Z , we set Z = Z avg tosome average atomic number as an approximation and interpolate values of S ( q, Z avg )from the tables given in [10]. 18he solid angle subtended by x and d may be approximated as:dΩ = A π × − ρ cos ϕ (1 − ρ cos ϕ ) + ρ sin ϕ , (66)where A is the detector area.The exponential terms in equation (62) account for the attenuation of the incomingand scattered rays. We approximate: e − (cid:82) Lsx µ ( E λ ,Z ) e − (cid:82) Lxd µ ( E s ,Z ) ≈ ≈ Z densities (e.g.nail varnish (Acetone), water, some plastic (polyethylene) etc.) and the rest may beclothes or air. Here also, as we are not worried about dosage, we can scan at highenergies (e.g. using a high voltage X-ray tube or high emission energy gamma raysource). We will simulate the error due to attenuation later, in section 4, and showthe effects of neglecting the attenuation in our reconstruction. Let our density f be supported on B (cid:15) ,(cid:15) ∩ U for some 0 < (cid:15) < (cid:15) < δ = − (cid:15) (cid:15) . If the source s is monochromatic ( E λ remains fixed), then the forwardCompton scattered intensity measured is: I ( r, α, β ) = ctD W k ( E λ ) d σ c dΩ ( r ) S w,p f ( r, α, β ) , (68)where p ∈ C ([0 , π ]) is defined by p ( ϕ ) = (cid:112) δ sin ϕ − δ sin ϕ , c is a constantthickness and: d σ c dΩ ( r ) = d σ dΩ ( E s , ω ) S ( q, Z avg ) (69)19epends only on r . Here the variable r ∈ [0 , 1] determines the scattered energy E s and( α, β ) ∈ S determines the source and detector position. The weighting w is given by: w ( r, ϕ ) = A (1 − ρ cos ϕ )4 π (cid:2) (1 − ρ cos ϕ ) + ρ sin ϕ (cid:3) (cid:2) ( ρ cos ϕ + 1) + ρ sin ϕ (cid:3) . (70)where ρ = (cid:112) r sin ϕ − 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 physicalmodelling 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 Comptonscattered intensity I . Here we have a spectrum of incoming photon energies. See figure 5. The ComptonFigure 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 (cid:48) for a source and detectorposition ( α, β ) ∈ S may be written: I ( r (cid:48) , α, β ) = stD I w f ( r (cid:48) , α, β ) , (71)20ith the weighting: w ( r (cid:48) , t, ϕ ) = (cid:20) W k ( r (cid:48) , t ) d σ c dΩ ( r (cid:48) , t ) (cid:21) w ( t, ϕ ) = w ( r (cid:48) , t ) w ( t, ϕ ) , (72)where w = w as in equation (70). Here the scattering probabilty d σ c dΩ and the photonflux W k are dependant on r (cid:48) and the integration variable t as in subsection 2.1. Whilethe weighting w is separable as in Theorem 4, the incident photon flux W k ( r (cid:48) , r (cid:48) ) = W k ( E m ) = 0 for all r (cid:48) , where E m is the maximum spectrum energy. Hence w ( r (cid:48) , r (cid:48) ) ≡ w fails to meet the conditions of Theorem 4. To deal with this, let: w avg ( r (cid:48) ) = 1 r (cid:48) (cid:90) r (cid:48) w ( r (cid:48) , t )d t (73)and let w app ( r (cid:48) , t, ϕ ) = w avg ( r (cid:48) ) w ( t, ϕ ). Then, although we sacrifice some accuracyin our forward model, w app would satisfy the conditions given in Theorem 4 and wecan obtain an analytic expression for the density. The error in our approximation isbounded by: | w avg ( r (cid:48) ) − w ( r (cid:48) , t ) | ≤ (cid:12)(cid:12)(cid:12)(cid:12) max t ∈ [0 ,r (cid:48) ] w ( r (cid:48) , t ) − min t ∈ [0 ,r (cid:48) ] w ( r (cid:48) , t ) (cid:12)(cid:12)(cid:12)(cid:12) (74)for all r (cid:48) ≥ 0, 0 ≤ t ≤ r (cid:48) . 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 (cid:48) ], the error in w app wouldbe small. Here we provide reconstructions of a test phantom density at a low resolution usingsimulated datasets of the unweighted spindle and spindle interior transforms, and sim-ulate noise as additive psuedo Gaussian noise. To simulate data for each transform wediscretize the integrals in equations (7) and (47), and solve the least squares problem:arg min x (cid:107) Ax − b (cid:107) + λ (cid:107) x (cid:107) (75)for some regularisation parameter λ > 0, where A is the discrete forward operator ofthe transform considered, x is the vector of density pixel values and b = Ax is thesimulated transform data. We simulate perturbed data b (cid:15) via: b (cid:15) = b + (cid:15) G (cid:107) b (cid:107)√ n (76)21here G is a pseudo random vector of samples from the standard Gaussian distibutionand n is the number of entries in b . The proposed noise model has the property that (cid:107) b − b (cid:15) (cid:107) / (cid:107) b (cid:107) ≈ (cid:15) , so a noise level of (cid:15) is (cid:15) × × × 50 pixels and a hollow ball, some stairs and a low density block with a metalsheet passing through it are placed in upper unit hemisphere. Slice profiles are givenin the right hand figure to display the metal sheet and the hollow shell. We sample25 r values (corresponding to spindles with heights H ∈ { . 02 + 0 . i : 0 ≤ i < } ),45 α values α ∈ { πi : 0 ≤ i < } and 45 β values β ∈ { π (1 + i ) : 0 ≤ i < } .So we have an underdetermined sparse system matrix A with 50625 rows and 125000columns. To reconstruct we solve the regularised problem (75) using the conjugategradient least squares algorithm (CGLS) and pick our regularisation parameter λ viaa manual approach. In each reconstruction presented, any negative values are set tozero and the iteration number and noise level are given in the figure caption.In figure 7 we have presented a reconstruction of the test phantom in the absenceof added noise. Given the symmetries involved in our geometry, the density has beenrotated and reflected in the xy plane in the reconstruction. To better visualise thereconstruction we set our reconstructed images to 0 in the lower half space. In figure8 we display the same zero noise reconstruction but with the voxel values set to zeroin the lower hemisphere.Figures 8–13 show test phantom reconstructions from spindle and spindle interiortransform data with varying levels of added noise. In the absence of noise (figures 8and 11) the reconstructions are ideal. However in the presence of noise we notice aharsher degradation in image quality when reconstructing with spindle interior data.This is because the spindle interior problem is inherently more ill posed due to theextra differentiation step required in the inversion process. At a low noise level (figure12), 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 metalsheet start to become lost in the reconstruction. At a higher noise level (figure 13), theartefacts in the reconstruction are severe, the ball no longer appears to be hollow andthe metal sheet fails to reconstruct. So, although the practical advantages of using apolychromatic source such as an X-ray tube or linear accelerator are clear (i.e. reduced22ata acquisition time, easy and safe to use etc.), a low noise level would need to bemaintained to achieve a satisfactory image quality. The spindle transform inversionperforms better when there is added noise. At a low noise level (figure 9), the size andshape 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 hardto identify. We notice, in the presence of noise, that the thinner object (namely themetal sheet) is hardest to reconstruct. This is as we’d expect since smaller densitiesare determined by the higher frequency harmonic components which degrade fasterwith noise.To simulate data for the reconstructions presented in figures 6–13 we applied thediscrete 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 inour simulated data and see how this effects the quality of our reconstruction. For ourfirst example, we set the size of the scanning cube to be 20 × × (each pixelis 4 × × ) and the phantom materials are polyethylene (low density block),water (the stairs), rubber (the ball) and Aluminium (the metal sheet). We model thesource as Co60, a monochromatic gamma ray source widely used in security screeningapplications (e.g. screening of freight shipping containers), which has an emissionenergy of 2824keV. We also add 1% Gaussian noise as in equation (76) to simulaterandom error as well as the systematic error due to attenuation effects. Our resultsare presented in figure 14. Here the phantom reconstruction is satisfactory with onlyminor artefacts appearing in the image. So if we scan a low effective Z target of asmall enough size with a high energy source, the effects due to attenuation can beneglected while maintaining a satisfactory image quality. For our next example, weset the scanning cube size to be 50 × × (each pixel is 1 × × ) and thephantom 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 againwe add a further 1% Gaussian noise to the simulated data after we have accountedfor 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 artefactsin the reconstruction are more severe. 23igure 6: Test phantom.Figure 7: Test phantom reconstruction from spindle transform data no noise, 2000iterations. 24igure 8: Test phantom reconstruction from spindle transform data no noise, set tozero in lower hemisphere, 2000 iterations.Figure 9: Test phantom reconstruction from spindle transform data, noise level (cid:15) =0 . 01, 2000 iterations. 25igure 10: Test phantom reconstruction from spindle transform data, noise level (cid:15) =0 . 05, 2000 iterations.Figure 11: Test phantom reconstruction from spindle interior data no noise, 2000iterations. 26igure 12: Test phantom reconstruction from spindle interior data, noise level (cid:15) = 0 . (cid:15) = 0 . Z test phantom reconstruction from spindle transformdata with attenuation effects added and noise level (cid:15) = 0 . 01, 2000 iterations.Figure 15: Large, high effective Z test phantom reconstruction from spindle transformdata with attenuation effects added and noise level (cid:15) = 0 . 01, 2000 iterations.28 Conclusions and further work We have presented a new acquisition geometry for three dimensional density recon-struction in Compton imaging with a monochromatic source and introduced a newspindle transform and a generalization of the spindle transform for the surfaces ofrevolution of a class of symmetric C curves. The generalized spindle transform wasshown to be injective on the domain of smooth functions f supported on a the inter-section of a hollow ball with the upper half space x > 0. In section 2 it was shownthat our problem could be decomposed into a set of one dimensional inverse problemsto solve for the harmonic coefficients of a given density, which we then went on to solvevia the explicit inversion of a class of Volterra integral operators. Later in section 2.1we considered the problem for a polychromatic source and introduced a new spindleinterior transform and proved its injectivity on the set of smooth functions compactlysupported on the intersection of a hollow ball and x > 0. We also considered theexterior problem for backscattered photons with a monochromatic and polychromaticphoton source, where in section 2.2 we introduced a new apple and apple interiortransform. Their injectivity was proven on the set of smooth functions compactlysupported on the intersection of the exterior of the unit ball and x > 0. We note thatin Palamodov’s paper on generalized Funk transforms [14], although he provides anexplicit inverse for a fairly general family of integral transforms over surfaces in threedimensional 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 for-ward problem for both a monochromatic and polychromatic source. In the monochro-matic source case, we found that the Compton scattered intensity resembled a weightedspindle transform (as in Theorem 2) which could be solved explicitly via repeated ap-proximations. In the polychromatic source case, with a more accurate forward model,we found that an analytic reconstruction was not possible. So we suggested a simpli-fied model in order to obtain an analytic expression for the density, and gave somesimple error estimates for our approximation.Test phantom reconstructions were presented in section 4 using simulated datasetsfrom spindle and spindle interior transform data with varying levels of added pseudorandom Gaussian noise. When reconstructing with spindle transform data the recon-29tructions were satifactory for noise levels up to 5%. We saw a harsher reduction inimage quality in the presence of added noise when reconstructing from spindle inte-rior data. This was as expected given the increased instability of the spindle interiorproblem.In the latter part of section 4, we provided further reconstructions of our testphantom image in the presence of a systematic error due to the attenuative effectsof 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 neglectedwhile maintaining a good image quality. We also gave an example where the targetmaterials were larger and higher effective Z . Here the effects due to attenuation weresignificant 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 workwe aim to analyse the spindle transform from a microlocal perspective to investigatewhether this can shed some light on its stability. Although the injectivity of theexterior problem is covered here, simulations of an exterior density reconstruction areleft for future work. Acknowledgements The authors would like to thank Rapiscan and the EPSRC for the CASE award fundingthis project, and WL would also like to thank the Royal Society for the WolfsonResearch Merit Award that also contributed to this project.