Three dimensional Compton scattering tomography

TL;DR

This paper introduces a novel 3D X-ray Compton imaging method using a new acquisition geometry, proving injectivity of associated transforms and demonstrating its effectiveness through simulated reconstructions with noise and error modeling.

Contribution

It develops a new 3D inverse problem framework for electron density reconstruction from spindle tori integrals and proves injectivity of related transforms, including polychromatic source cases.

Findings

Injectivity of the generalized spindle torus transform is established.

Explicit inversion formulas for harmonic coefficients are derived.

Simulated reconstructions demonstrate robustness to noise and systematic errors.

Abstract

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.

Figures

Equations

$$E_s=\frac{E_{\lambda }}{1+\left(E_{\lambda }/E_0\right)\left(1-\cos \omega \right)}.$$

$$T_r=\left\{(x_1,x_2,x_3)\in {\mathbb{R}}^3:{\left(r-\sqrt{x_1^2+x_2^2}\right)}^2+x_3^2=1+r^2\right\}$$

$$R^2={\rho }^2+r^2+2r\rho \sin \phi$$

$$\rho =\sqrt{r^2{\sin }^2\phi +1}-r\sin \phi .$$

$$\mathrm{d}v=\rho \sqrt{\frac{1+r^2}{1+r^2{\sin }^2\phi }}\mathrm{d}\phi$$

$$\mathrm{d}A=\rho \sin \phi \mathrm{d}v\mathrm{d}\theta .$$

$$\mathcal{S}f(r,\alpha ,\beta )={\int }_0^{2\pi }{\int }_0^{\pi }{\rho }^2\sin \phi \sqrt{\frac{1+r^2}{1+r^2{\sin }^2\phi }}(h\cdot F)\left(\rho ,\theta ,\phi \right){\mid }_{\rho =\sqrt{r^2{\sin }^2\phi +1}-r\sin \phi }\mathrm{d}\phi \mathrm{d}\theta$$

$$\rho (r,\phi )=p({\sin }^{-1}(r\sin \phi )).$$

$${\mathcal{S}}_{w,p}f(r,\alpha ,\beta )={\int }_0^{2\pi }{\int }_0^{\pi }w(r,\phi )\rho \sin \phi \sqrt{{\rho }^2+{\left(\frac{\mathrm{d}\rho }{\mathrm{d}\phi }\right)}^2}(h\cdot F)\left(\rho ,\theta ,\phi \right){\mid }_{\rho =\rho (r,\phi )}\mathrm{d}\phi \mathrm{d}\theta .$$

$$Y_l^m(\theta ,\phi )=(-1{)}^m\sqrt{\frac{(2l+1)(l-m)!}{4\pi (l+m)!}}{P}_l^m(\cos \phi )e^{\mathrm{i}m\theta }$$

$$P_l^m(x)=(-1{)}^m(1-x^2{)}^{m/2}\frac{{\mathrm{d}}^m}{\mathrm{d}{x}^m}{P}_l(x)$$

$$P_l(x)=\frac{1}{2^l}\sum _{k=0}^l{\left(\genfrac{}{}{0ex}{}{l}{k}\right)}^2(x-1{)}^{l-k}(x+1{)}^k,$$

$$⟨f,g⟩={\int }_{S^2}f\bar{g}\mathrm{d}\tau ={\int }_0^{\pi }{\int }_0^{2\pi }f(\theta ,\phi )\bar{g}(\theta ,\phi )\sin \phi \mathrm{d}\theta \mathrm{d}\phi .$$

$${\int }_{S^2}{Y}_l^m{\bar{Y}}_{l^{\prime }}^{m^{\prime }}\mathrm{d}\tau ={\int }_0^{\pi }{\int }_0^{2\pi }{Y}_l^m(\theta ,\phi ){\bar{Y}}_{l^{\prime }}^{m^{\prime }}(\theta ,\phi )\sin \phi \mathrm{d}\theta \mathrm{d}\phi ={\delta }_{l{l}^{\prime }}{\delta }_{m{m}^{\prime }},$$

$$F_{lm}={\int }_{S^2}F{\bar{Y}}_l^m\mathrm{d}\tau ,$$

$$F_N=\sum _{0\le l\le N}\sum _{|m|\le l}{F}_{lm}{Y}_l^m$$

$${\mathcal{S}}_{w,p}{f}_{lm}(r)=2\pi {\int }_0^{\pi }w(r,\phi )\rho \sin \phi \sqrt{{\rho }^2+{\left(\frac{\mathrm{d}\rho }{\mathrm{d}\phi }\right)}^2}{F}_{lm}(\rho )P_l(\cos \phi ){\mid }_{\rho =\rho (r,\phi )}\mathrm{d}\phi$$

$${\mathcal{S}}_{w,p}{f}_{lm}={\int }_{S^2}({\mathcal{S}}_{w,p}f){\bar{Y}}_l^m\mathrm{d}\tau$$

$$\begin{array}{rl}{\mathcal{S}}_{w,p}f(r,\alpha ,\beta )& ={\int }_0^{2\pi }{\int }_0^{\pi }w(r,\phi )\rho \sin \phi \sqrt{{\rho }^2+{\left(\frac{\mathrm{d}\rho }{\mathrm{d}\phi }\right)}^2}(h\cdot F)\left(\rho ,\theta ,\phi \right){\mid }_{\rho =\rho (r,\phi )}\mathrm{d}\phi \mathrm{d}\theta \\ & =\sum _{l\in \mathbb{N}}\sum _{|m|\le l}{\int }_0^{2\pi }{\int }_0^{\pi }w(r,\phi )\rho \sin \phi \sqrt{{\rho }^2+{\left(\frac{\mathrm{d}\rho }{\mathrm{d}\phi }\right)}^2}{F}_{lm}(\rho )(h\cdot Y_l^m)(\theta ,\phi )\mathrm{d}\phi \mathrm{d}\theta \end{array}$$

$$(h\cdot Y_l^m)(\tau )=\sum _{|n|<l}{Y}_l^n(\tau )D_{n,m}^{(l)}(h^{-1})$$

$$D_{n,m}^{(l)}(h)=D_{n,m}^{(l)}(U(\alpha )V(\beta )U(\gamma ))=e^{-in\gamma }{d}_{n,m}^{(l)}(\cos \beta )e^{-im\alpha }$$

$$d_{n,0}^{(l)}(\cos \beta )=(-1{)}^n\sqrt{\frac{(l-n)!}{4\pi (l+n)!}}{P}_l^n(\cos \beta ).$$

$$\begin{array}{rl}{\mathcal{S}}_{w,p}f(r,\alpha ,\beta )& =2\pi \sum _{l\in \mathbb{N}}\sum _{|m|\le l}\sqrt{\frac{(2l+1)}{4\pi }}{D}_{0,m}^{(l)}(h^{-1}){\int }_0^{\pi }w(r,\phi )\rho \sin \phi \sqrt{{\rho }^2+{\left(\frac{\mathrm{d}\rho }{\mathrm{d}\phi }\right)}^2}{F}_{lm}(\rho )P_l(\cos \phi )\mathrm{d}\phi .\end{array}$$

$$Y_l^m(\alpha ,\beta )=\sqrt{\frac{(2l+1)}{4\pi }}{\bar{D}}_{m,0}^{(l)}(h).$$

$$\begin{array}{rl}{\mathcal{S}}_{w,p}f(r,\alpha ,\beta )& =2\pi \sum _{l\in \mathbb{N}}\sum _{|m|\le l}{Y}_l^m(\alpha ,\beta ){\int }_0^{\pi }w(r,\phi )\rho \sin \phi \sqrt{{\rho }^2+{\left(\frac{\mathrm{d}\rho }{\mathrm{d}\phi }\right)}^2}{F}_{lm}(\rho )P_l(\cos \phi )\mathrm{d}\phi .\end{array}$$

$$\begin{array}{rl}{\mathcal{S}}_{w,p}{f}_{lm}(r)& ={\int }_{S^2}{\mathcal{S}}_{w,p}f{\bar{Y}}_l^m\mathrm{d}\tau \\ & =2\pi {\int }_0^{\pi }w(r,\phi )\rho \sin \phi \sqrt{{\rho }^2+{\left(\frac{\mathrm{d}\rho }{\mathrm{d}\phi }\right)}^2}{F}_{lm}(\rho )P_l(\cos \phi ){\mid }_{\rho =\rho (r,\phi )}\mathrm{d}\phi \end{array}$$

$$\begin{array}{rl}\frac{1}{4\pi }{\mathcal{S}}_{w,p}{f}_{lm}(r)& ={\int }_0^{\frac{\pi }{2}}w(r,\phi )\rho \sin \phi \sqrt{{\rho }^2+{\left(\frac{\mathrm{d}\rho }{\mathrm{d}\phi }\right)}^2}{F}_{lm}(\rho )P_l(\cos \phi ){\mid }_{\rho =\rho (r,\phi )}\mathrm{d}\phi \\ & ={\int }_0^{\frac{\pi }{2}}w(r,\phi )g(\rho )\sin \phi \sqrt{g(\rho {)}^2+r^2{\cos }^2\phi g^{\prime }(\rho {)}^2}{F}_{lm}(g(\rho ))P_l(\cos \phi ){\mid }_{\rho =r\sin \phi }\mathrm{d}\phi .\end{array}$$

$$\frac{1}{4\pi }{\mathcal{S}}_{w,p}{f}_{lm}(r)={\int }_0^r\frac{F_{lm}(g(\rho ))K_l(r,\rho )}{\sqrt{r-\rho }}\mathrm{d}\rho$$

$$K_l(r,\rho )=W(r,\rho )\frac{\rho g(\rho )\sqrt{g(\rho {)}^2+(r^2-{\rho }^2)g^{\prime }(\rho {)}^2}}{r\sqrt{r+\rho }}{P}_l\left(\sqrt{1-\frac{{\rho }^2}{r^2}}\right).$$

$$\frac{1}{4\pi }{\int }_0^z\frac{{\mathcal{S}}_{w,p}{f}_{lm}(r)}{\sqrt{z-r}}\mathrm{d}r={\int }_0^z{F}_{lm}(g(\rho ))\left[{\int }_{\rho }^z\frac{K_l(r,\rho )}{\sqrt{z-r}\sqrt{r-\rho }}\mathrm{d}r\right]\mathrm{d}\rho$$

Full text

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.

1 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:

[TABLE]

Here $E_{s}$ is the energy of the scattered photon which had an initial energy $E_{\lambda}$, $\omega$ is the scattering angle and $E_{0}\approx 511$keV 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 $\omega\leq\pi/2$. Conversely, backscatter refers to the scattering of photons at angles $\omega>\pi/2$.

If the photon source is monochromatic ($E_{\lambda}$ 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.

A spindle torus is the surface of revolution of a circular arc. Specifically we define:

[TABLE]

to be the spindle torus, radially symmetric about the $x_{3}$ axis, with tube centre offset $r\geq 0$ and tube radius $\sqrt{1+r^{2}}$. Let $B_{0,1}$ denote the unit ball in $\mathbb{R}^{3}$. Then we define the spindle $S_{r}=T_{r}\cap B_{0,1}$ as the interior of the torus $T_{r}$ and we define the apple $A_{r}=T_{r}\backslash 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 $\omega$. 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 $\omega\leq\pi/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 ($\omega>\pi/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 scatter problem and introduce a generalization of the spindle transform which gives the integrals of a function over the surfaces of revolution of a particular class of symmetric curves. We prove the injectivity of the generalized spindle transform on the domain of smooth functions compactly supported on the intersection of a hollow ball with the upper half space $x_{3}>0$. We show that our problem can be decomposed as a set of one dimensional inverse problems, which we then solve to provide an explicit expression for the harmonic coefficients of $f$. In section 2.1 we introduce a new toric interior transform for the polychromatic forward scatter problem and prove its injectivity on the domain of smooth functions compactly supported on the intersection of a hollow ball and $x_{3}>0$. A new apple and apple interior transform are also introduced in section 2.2 for the monochromatic and polychromatic backscatter problem. Their injectivity is proven on the domain of smooth functions compactly supported on the intersection of the exterior of the unit ball and $x_{3}>0$.

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.

2 A spindle transform

We will now parameterize the set of points on a spindle $S_{r}$ in terms of spherical coordinates $(\rho,\theta,\varphi)$ 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:

[TABLE]

and hence:

[TABLE]

Let $B_{\epsilon_{1},\epsilon_{2}}=\{x\in\mathbb{R}^{3}:\epsilon_{1}<|x|<\epsilon_{2}\}$ denote the set of points on a hollow ball with inner radius $\epsilon_{1}$ and outer radius $\epsilon_{2}$, and let $S^{2}$ denote the unit sphere. Let $Z=\mathbb{R}^{+}\times S^{2}$. For a function $f:\mathbb{R}^{3}\to\mathbb{R}$, let $F:Z\to\mathbb{R}$ be iits polar form $F(\rho,\theta,\varphi)=f(\rho\sin\varphi\cos\theta,\rho\sin\varphi\sin\theta,\rho\cos\varphi)$. We parameterize $h\in\text{SO}(3)$ in terms of Euler angles $\alpha$ and $\beta$, $h=U(\alpha)V(\beta)$, where $U$ and $V$ are rotations about the $x_{3}$ and $x_{2}$ axis respectively. We define an action of the rotation group $\text{SO}(3)$ on the set of real valued functions $f$ on $\mathbb{R}^{3}$ in the natural way, and define an action on the polar form as $(h\cdot F)=(h\cdot f)$.

The arc element for the circular arc $C$ is given in [3] as:

[TABLE]

and the area element for the spindle $S_{r}$ is:

[TABLE]

We define the spindle transform $\mathcal{S}:C_{0}^{\infty}(\mathbb{R}^{3})\to C_{0}^{\infty}(Z)$ as:

[TABLE]

where $h=U(\alpha)V(\beta)$.

This transform belongs to a larger class of integral transforms as we will now show. Let $p\in C^{1}([0,\frac{\pi}{2}])$ be a curve parameterized by a colatitude $\varphi\in[0,\frac{\pi}{2}]$. Then we define the class of symmetric curves $\rho\in C^{1}([0,1]\times[0,\pi])$ by:

[TABLE]

From this we define the generalized spindle transform $\mathcal{S}_{w,p}:C_{0}^{\infty}(\mathbb{R}^{3})\to C_{0}^{\infty}([0,1]\times S^{2})$ as:

[TABLE]

where the weighting $w$ is dependant only on $r$ and $\varphi$.

We now aim to prove injectivity of the generalized spindle transform $\mathcal{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\geq 0,\ |m|\leq l$, we define the spherical harmonics $Y_{l}^{m}$ as:

[TABLE]

where,

[TABLE]

and

[TABLE]

are Legendre polynomials of degree $l$. The spherical harmonics $Y_{l}^{m}$ form an orthonormal basis for $L^{2}(S^{2})$, with the inner product:

[TABLE]

So we have:

[TABLE]

where $\delta$ denotes the Kroneker delta.

From [9] we have the theorem:

Theorem 1**.**

Let $f\in C^{\infty}(\mathbb{R}^{3})$ and let:

[TABLE]

where $\mathrm{d}\tau$ is the surface measure on the sphere. Then the series:

[TABLE]

converges uniformly absolutely on compact subsets of $Z$ to $F$.

We now show how the problem of reconstructing a density $f$ from its integrals $\mathcal{S}_{w,p}f$ for some curve $p\in C^{1}([0,\frac{\pi}{2}])$ can be broken down into a set of one dimensional inverse problems to solve for the harmonic coefficients $f_{lm}$:

Lemma 1**.**

Let $f\in C_{0}^{\infty}(\mathbb{R}^{3})$ and let $p\in C^{1}([0,\frac{\pi}{2}])$ be a curve, then:

[TABLE]

where,

[TABLE]

and $P_{l}$ is a Legendre polynomial degree $l$.

Proof.

Let $\tau\in S^{2}$ and let $h\in\text{SO}(3)$ act on a spherical harmonic $Y_{l}^{m}$ by $(h\cdot Y_{l}^{m})(\tau)=Y_{l}^{m}(h\tau)$. Then we have:

[TABLE]

where $h=U(\alpha)V(\beta)$.

We may write $h\cdot Y_{l}^{m}$ as a linear combination of spherical harmonics of the same degree [12, page 209]:

[TABLE]

where the block diagonal entries $D_{n,m}^{(l)}$ are defined:

[TABLE]

and $d_{n,m}^{(l)}$ is given, for $m=0$, by:

[TABLE]

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:

[TABLE]

Since the regular representation of $\text{SO}(3)$ is unitary, we have $D_{0,m}^{(l)}(h^{-1})=\bar{D}_{m,0}^{(l)}(h)$ and from the above we have:

[TABLE]

It follows that:

[TABLE]

From which we have:

[TABLE]

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\in\mathbb{N}$, $|m|\leq l$. Given the antisymmetry of the Legendre polynomials for odd $l$, we see that $\mathcal{S}_{w,p}f_{lm}=0$ for odd $l$, for any given $f\in C_{0}^{\infty}(\mathbb{R}^{3})$ and curve $p\in C^{1}([0,\frac{\pi}{2}])$. So for the moment we focus on the reconstruction of the coefficients $f_{lm}$ for even $l$.

We now have our main theorem:

Theorem 2**.**

Let $p\in C^{1}([0,\frac{\pi}{2}])$ be such that $p(\varphi)\neq 0$ for all $\varphi\in[0,\frac{\pi}{2}]$, and let $g\in C^{1}([0,1))$ be defined as $g(x)=p(\sin^{-1}x)$. Let $A=\{\rho\in[0,1):g(\rho)\in[0,\epsilon_{1}]\cup[\epsilon_{2},\epsilon_{3}]\}$ and let $w:[0,1]\times[0,\pi]\to\mathbb{R}$ be a weighting such that $W(r,\rho)=w\left(r,\sin^{-1}(\rho/r)\right)$ and its first order partial derivative with respect to $r$ is continuous on $\{(r,\rho):r\in A,\textnormal{min}(A)\leq\rho\leq r\}$ and $W(\rho,\rho)\neq 0$ on $A$. Let $f\in C_{0}^{\infty}(B_{0,\epsilon_{1}}\cup B_{\epsilon_{2},\epsilon_{3}})$, where $\epsilon_{1}<p(0)<\epsilon_{2}$. Then $\mathcal{S}_{w,p}f$ determines the harmonic coefficients $F_{lm}$ uniquely for $\rho\in g([0,1))$ for all $l\in\{2k:k\in\mathbb{N}\},|m|\leq l$.

Proof.

Given the symmetry of the Legendre polynomials $P_{l}$ for even $l$, we have:

[TABLE]

After making the substitution $\rho=r\sin\varphi$ , equation (27) becomes:

[TABLE]

a Volterra integral equation of the first kind with weakly singular kernel, where:

[TABLE]

As $F_{lm}\circ g$ is zero for $\rho$ close to [math], given our prior assumptions regarding $W$ and given that $g\in 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}\circ g$.

Multiplying both sides of equation (28) by $1/\sqrt{z-r}$ and integrating with respect to $r$ over the interval $[0,z]$, yields:

[TABLE]

after changing the integration order. Making the substitution $r=\rho+(z-\rho)t$, gives:

[TABLE]

from which we have:

[TABLE]

where,

[TABLE]

So by our assumptions that $p$ is non zero and $W$ is non zero on the diagonal, $Q_{l}(z,z)\neq 0$ on the support of $F_{lm}\circ g$.

Differentiating both sides of equation (30) with respect to $z$ and rearranging gives:

[TABLE]

which is a Volterra type integral equation of the second kind, where:

[TABLE]

and,

[TABLE]

Given the continuity of $H_{l}$ on the support of $F_{lm}\circ g$ and the continuity of $F_{lm}\circ g$ on $[0,1]$, the Neumann series associated with equation (34) converges and we may write our solution:

[TABLE]

where the resolvent kernel,

[TABLE]

is defined by,

[TABLE]

Since the series converges, the solution is unique and we may reconstruct $F_{lm}$ explicitly for $\rho\in g([0,1))$. ∎

Here as we can only reconstruct the coefficients $F_{lm}$ for even $l$, for a general $f\in C_{0}^{\infty}(\mathbb{R}^{3})$ and curve $p\in C^{1}([0,\frac{\pi}{2}])$, the spindle data $\mathcal{S}_{w,p}f$ is insufficient to recover $f$ uniquely for $|x|\in p([0,\frac{\pi}{2}))$. However if we consider those functions $f\in C_{0}^{\infty}(\mathbb{R}^{3})$ whose support lies in the upper half space $x_{3}>0$ we find that uniqueness is possible, as the following theorem shows:

Theorem 3**.**

Let $U=\{(x_{1},x_{2},x_{3})\in\mathbb{R}^{3}:x_{3}>0\}$ and let $f\in C_{0}^{\infty}(U)$. Then the even coefficients $F_{lm}$ for $l\in\{2k:k\in\mathbb{N}\}=\mathbb{N}_{e},|m|\leq l$ determine $F$ uniquely and:

[TABLE]

for $\rho\geq 0$, $0\leq\theta\leq 2\pi$ and $0\leq\varphi\leq\frac{\pi}{2}$.

Proof.

We can write $F$ as the sum of its even and odd coefficients:

[TABLE]

Then from our assumption, we have:

[TABLE]

for all $\rho\geq 0$, $0\leq\theta\leq 2\pi$ and $\frac{\pi}{2}\leq\varphi\leq\pi$, where:

[TABLE]

We have:

[TABLE]

The associated Legendre polynomials $P_{l}^{m}$ are symmetric when $l+m$ is even and antisymmetric otherwise. It follows that:

[TABLE]

for $\rho\geq 0$, $0\leq\theta\leq 2\pi$ and $0\leq\varphi\leq\frac{\pi}{2}$. The result follows. ∎

Corollary 1**.**

Let $f\in C_{0}^{\infty}(B_{\epsilon_{1},\epsilon_{2}}\cap U)$ for some $0<\epsilon_{1}<\epsilon_{2}<1$. Let $\epsilon_{1}\leq\epsilon\leq\epsilon_{2}$ and let $\delta=\frac{1-\epsilon^{2}}{2\epsilon}$. Then $\mathcal{S}f$ known for $0\leq r\leq\delta$ and for all $(\alpha,\beta)\in S^{2}$ determines $f$ uniquely for $\epsilon\leq|x|\leq 1$.

Proof.

Let $p(\varphi)=\sqrt{1+\delta^{2}\sin^{2}\varphi}-\delta\sin\varphi$ and let $w\equiv 1$. Then $\mathcal{S}f(\delta r,\alpha,\beta)=\mathcal{S}_{w,p}f(r,\alpha,\beta)$ for $r\in[0,1]$, $(\alpha,\beta)\in S^{2}$. The result follows from Theorems 2 and 3. ∎

2.1 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 $\frac{E_{m}}{1+2E_{m}/E_{0}}<E_{s}<E_{m}$, the set of scatterers is the union of spindle tori corresponding to scattering angles $\omega$ in the range $0<\omega<\cos^{-1}\left(1-\frac{E_{0}\left(E_{m}-E_{s}\right)}{E_{s}E_{m}}\right)$ (corresponding to energies $E$ in the range $E_{s}<E<E_{m}$). That is, the set of scatterers is a torus interior:

[TABLE]

where $r$ is determined by the scattered energy measured $E_{s}$. See [1] for an explanation of the two dimensional case.

With this in mind we define the spindle interior transform $\mathcal{I}:C_{0}^{\infty}(\mathbb{R}^{3})\to C_{0}^{\infty}(Z)$ as:

[TABLE]

and we have the uniqueness theorem for a weighted spindle interior transform:

Theorem 4**.**

Let $f\in C_{0}^{\infty}(B_{\epsilon_{1},\epsilon_{2}}\cap U)$ for some $0<\epsilon_{1}<\epsilon_{2}<1$ and and let $\delta_{2}=\frac{1-\epsilon_{2}^{2}}{2\epsilon_{2}}$ and $\delta_{1}=\frac{1-\epsilon_{1}^{2}}{2\epsilon_{1}}$. Let $\tilde{f}$ be defined as:

[TABLE]

Define the weighted interior transform:

[TABLE]

where $h=U(\alpha)V(\beta)$ and $r^{\prime}=1/r$. Let us suppose that the weighting $w$ satisfies the following:

$w$* can be decomposed as $w(r^{\prime},t,\varphi)=w_{1}(r^{\prime},t)w_{2}(t,\varphi)$.* 2. 2.

$W_{2}(t,\rho)=w_{2}\left(1/t,\sin^{-1}(\rho/t)\right)$* and its first order partial derivative with respect to $t$ is continuous on the triangle $T=\{(t,\rho)\in\mathbb{R}^{2}:\epsilon_{1}\leq t\leq\epsilon_{2},\epsilon_{1}\leq\rho\leq t\}$ and $W_{2}(\rho,\rho)\neq 0$ on $[\epsilon_{1},\epsilon_{2}]$.* 3. 3.

$w_{1}(r^{\prime},t)$* and its first order partial derivatives are bounded on $0\leq r^{\prime}<\infty$ and $w_{1}(r^{\prime},r^{\prime})\neq 0$ for $0\leq r^{\prime}<\infty$.*

Then $\mathcal{I}_{w}f$ known for all $r^{\prime}\in[0,\infty)$ and for all $(\alpha,\beta)\in S^{2}$ determines $f$ uniquely for $0<|x|<1$.

Proof.

We have:

[TABLE]

where,

[TABLE]

Differentiating both sides of equation (50) with respect to $r^{\prime}$ and rearranging yields:

[TABLE]

where,

[TABLE]

and,

[TABLE]

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<\infty$. From which we have:

[TABLE]

for $t\in(0,1]$, where $w_{3}(t,\varphi)=w_{2}\left(\frac{1}{\delta_{1}t},\varphi\right)$ and $p(\varphi)=\sqrt{1+\delta_{1}^{2}\sin^{2}\varphi}-\delta_{1}\sin\varphi$. 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. ∎

Corollary 2**.**

Let $f\in C_{0}^{\infty}(B_{\epsilon_{1},\epsilon_{2}}\cap U)$ for some $0<\epsilon_{1}<\epsilon_{2}<1$. Then $\mathcal{I}f$ as defined in equation (47) known for all $r\in[0,\infty)$ and for all $(\alpha,\beta)\in S^{2}$ determines $f$ uniquely for $0<|x|<1$.

Proof.

Let $\tilde{f}$ be defined as in Theorem 4. Then, after making the substitution $\rho=\sqrt{\frac{\sin^{2}\varphi}{t^{2}}+1}-\frac{\sin\varphi}{t}$ in equation (47), we have:

[TABLE]

for all $r\in(0,\infty)$ and for all $(\alpha,\beta)\in S^{2}$, where the weighting $w\equiv 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.

2.2 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.

2.2.1 An apple transform

For backscattered photons (for scattering angles $\omega>\pi/2$) the surface of scatterers is an apple $A_{r}$. Refer back to figure 2. The spherical coordinates $(\rho,\theta,\varphi)$ of the points on $A_{r}$ can be parameterized as follows:

[TABLE]

We define the apple transform $\mathcal{A}:C_{0}^{\infty}(\mathbb{R}^{3})\to C_{0}^{\infty}(Z)$ as:

[TABLE]

where $h=U(\alpha)V(\beta)$. We have the uniqueness theorem for the apple transform for functions supported on the exterior of the unit ball:

Theorem 5**.**

Let $f\in C_{0}^{\infty}(B_{\epsilon_{1},\epsilon_{2}}\cap U)$ for some $1<\epsilon_{1}<\epsilon_{2}<\infty$. Let $\epsilon_{1}\leq\epsilon\leq\epsilon_{2}$ and let $\delta=\frac{\epsilon^{2}-1}{2\epsilon}$. Then $\mathcal{A}f$ known for $0\leq r\leq\delta$ and for all $(\alpha,\beta)\in S^{2}$ determines $f$ uniquely for $1\leq|x|\leq\epsilon$.

Proof.

Let $p(\varphi)=\sqrt{1+\delta^{2}\sin^{2}\varphi}+\delta\sin\varphi$ and let $w\equiv 1$. Then $\mathcal{A}f(\delta r,\alpha,\beta)=\mathcal{S}_{w,p}f(r,\alpha,\beta)$ for $r\in[0,1]$, $(\alpha,\beta)\in S^{2}$. The result follows from Theorems 2 and 3. ∎

2.2.2 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 $\mathcal{AI}:C_{0}^{\infty}(\mathbb{R}^{3})\to C_{0}^{\infty}(Z)$:

[TABLE]

and we have a theorem for its injectivity:

Theorem 6**.**

Let $f\in C_{0}^{\infty}(B_{\epsilon_{1},\epsilon_{2}}\cap U)$ for some $1<\epsilon_{1}<\epsilon_{2}<\infty$. Let $\epsilon_{1}\leq\epsilon\leq\epsilon_{2}$ and let $\delta=\frac{\epsilon^{2}-1}{2\epsilon}$. Then $\mathcal{AI}f$ known for $0\leq r\leq\delta$ and for all $(\alpha,\beta)\in S^{2}$ determines $f$ uniquely for $1\leq|x|\leq\epsilon$.

Proof.

Let $\tilde{f}$ be defined as:

[TABLE]

Then by Leibniz rule we have:

[TABLE]

So $\mathcal{AI}f$ known for $0\leq r\leq\delta$ and for $(\alpha,\beta)\in S^{2}$ determines $\mathcal{A}\tilde{f}$ on the same set after differentiating. The result follows from Theorem 5. ∎

3 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:

[TABLE]

where $Z$ denotes the atomic number, $\omega$ the scattering angle and $L_{sx}$ and $L_{xd}$ are the line segments connecting $s$ to $x$ and $x$ to $d$ respectively. $f(x)$ denotes the electron density (number of electrons per unit volume) at the scattering point $x$ and $\mathrm{d}V$ is the volume measure. $f$ is the quantity to be reconstructed.

Let $W_{k}(E_{\lambda})$ denote the incident photon flux (number of photons per unit area per unit time), energy $E_{\lambda}$, measured at a fixed distance $D$ from the source. Then the incident photon intensity $I_{0}$ can be written:

[TABLE]

where $t$ is the emission time.

The Klein-Nishina differential cross section $\mathrm{d}\sigma/\mathrm{d}\Omega$, is defined by:

[TABLE]

where $r_{0}$ is the classical electron radius and $\cos\omega=r/\sqrt{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\left(q,Z\right)$. Here $q=\frac{E_{\lambda}}{hc}\sin\left(\omega/2\right)$ is the momentum transferred by a photon with initial energy:

[TABLE]

scattering at an angle $\omega$, 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_{\text{avg}}$ to some average atomic number as an approximation and interpolate values of $S(q,Z_{\text{avg}})$ from the tables given in [10].

The solid angle subtended by $x$ and $d$ may be approximated as:

[TABLE]

where $A$ is the detector area.

The exponential terms in equation (62) account for the attenuation of the incoming and scattered rays. We approximate:

[TABLE]

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 ($\approx 1$MeV), 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.

3.0.1 The monochromatic case

Let our density $f$ be supported on $B_{\epsilon_{1},\epsilon_{2}}\cap U$ for some $0<\epsilon_{1}<\epsilon_{2}<1$ and let $\delta_{2}=\frac{1-\epsilon^{2}_{2}}{2\epsilon_{2}}$. If the source $s$ is monochromatic ($E_{\lambda}$ remains fixed), then the forward Compton scattered intensity measured is:

[TABLE]

where $p\in C^{1}([0,\pi])$ is defined by $p(\varphi)=\sqrt{1+\delta^{2}_{2}\sin^{2}\varphi}-\delta_{2}\sin\varphi$, $c$ is a constant thickness and:

[TABLE]

depends only on $r$. Here the variable $r\in[0,1]$ determines the scattered energy $E_{s}$ and $(\alpha,\beta)\in S^{2}$ determines the source and detector position. The weighting $w$ is given by:

[TABLE]

where $\rho=\sqrt{1+r^{2}\sin^{2}\varphi}-r\sin\varphi$. 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 $\mathcal{S}_{w,p}f$ as in Theorem 2 to obtain an analytic expression for $f$ in terms of the Compton scattered intensity $I$.

3.0.2 The polychromatic case

Here we have a spectrum of incoming photon energies. See figure 5.

The Compton scattered intensity measured for a tube centre offset $r=1/r^{\prime}$ for a source and detector position $(\alpha,\beta)\in S^{2}$ may be written:

[TABLE]

with the weighting:

[TABLE]

where $w_{2}=w$ as in equation (70). Here the scattering probabilty $\frac{\mathrm{d}\sigma_{c}}{\mathrm{d}\Omega}$ and the photon flux $W_{k}$ are dependant on $r^{\prime}$ 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^{\prime},r^{\prime})=W_{k}(E_{m})=0$ for all $r^{\prime}$, where $E_{m}$ is the maximum spectrum energy. Hence $w_{1}(r^{\prime},r^{\prime})\equiv 0$ and $w$ fails to meet the conditions of Theorem 4. To deal with this, let:

[TABLE]

and let $w_{\text{app}}(r^{\prime},t,\varphi)=w_{\text{avg}}(r^{\prime})w_{2}(t,\varphi)$. Then, although we sacrifice some accuracy in our forward model, $w_{\text{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:

[TABLE]

for all $r^{\prime}\geq 0$, $0\leq t\leq r^{\prime}$. So provided the changes in the incident photon energy $E_{s}$ ($E_{s}$ is determined by $t$) are negligible over the range $t\in[0,r^{\prime}]$, the error in $w_{\text{app}}$ would be small.

4 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:

[TABLE]

for some regularisation parameter $\lambda>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^{\epsilon}$ via:

[TABLE]

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^{\epsilon}\|/\|b\|\approx\epsilon$, so a noise level of $\epsilon$ is $\epsilon\times 100\%$ relative error.

We consider the test phantom displayed in figure 7. The unit cube is discretised into $50\times 50\times 50$ pixels and a hollow ball, some stairs and a low density block with a metal sheet passing through it are placed in upper unit hemisphere. Slice profiles are given in the right hand figure to display the metal sheet and the hollow shell. We sample 25 $r$ values (corresponding to spindles with heights $H\in\{0.02+0.04i:0\leq i<25\}$), 45 $\alpha$ values $\alpha\in\{\frac{2\pi i}{45}:0\leq i<45\}$ and 45 $\beta$ values $\beta\in\{\frac{\pi}{90}(1+i):0\leq i<45\}$. So we have an underdetermined sparse system matrix $A$ with 50625 rows and 125000 columns. To reconstruct we solve the regularised problem (75) using the conjugate gradient least squares algorithm (CGLS) and pick our regularisation parameter $\lambda$ via a manual approach. In each reconstruction presented, any negative values are set to zero 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 absence of added noise. Given the symmetries involved in our geometry, the density has been rotated and reflected in the $xy$ plane in the reconstruction. To better visualise the reconstruction we set our reconstructed images to 0 in the lower half space. In figure 9 we display the same zero noise reconstruction but with the voxel values set to zero in the lower hemisphere.

Figures 9–13 show test phantom reconstructions from spindle and spindle interior transform data with varying levels of added noise. In the absence of noise (figures 9 and 11) the reconstructions are ideal. However in the presence of noise we notice a harsher degradation in image quality when reconstructing with spindle interior data. 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 13), 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 7–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\times 20\times 20\text{cm}^{3}$ (each pixel is $4\times 4\times 4$mm3) 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 15. 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\times 50\times 50\text{cm}^{3}$ (each pixel is $1\times 1\times 1$cm3) 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.

5 Conclusions and further work

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 reconstructions 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.

Acknowledgements

The authors would like to thank Rapiscan and the EPSRC for the CASE award funding this project, and WL would also like to thank the Royal Society for the Wolfson Research Merit Award that also contributed to this project.