A joint reconstruction and lambda tomography regularization technique for energy-resolved x-ray imaging

In this paper we introduce new joint reconstruction and regularization techniques based on ideas in microlocal analysis and lambda tomography [13, 14] (see also [32] for related work). We consider the simultaneous reconstruction of the attenuation coefficient μ_E and electron density n_e from joint x-ray CT (transmission) and Compton scattered data, with particular focus on the parallel line segment x-ray scanner displayed in figures 1 and 2. The acquisition geometry in question is based on a new airport baggage scanner currently in development, and has the ability to measure x-ray CT and Compton data simultaneously. The line segment geometry was first considered in [54], where injectivity results are derived in Compton scattering tomography (CST). We provide a stability analysis of the CST problem of [54] here, from a microlocal perspective. The scanner depicted in figure 1 consists of a row of fixed, switched, monochromatic fan beam sources (S), a row of detectors (D_A) to measure the transmitted photons, and a second (slightly out of plane) row of detectors (D_C) to measure Compton scatter. The detectors are assumed to energy-resolved, a common assumption in CST [35, 41, 43, 51, 55], and the sources are fan-beam (in the plane) with opening angle π (so there is no restriction due to cropped fan-beams).

The attenuation coefficient relates to the x-ray transmission data by the Beer–Lambert law, $\mathrm{log}\left(\frac{{I}_{0}}{{I}_{A}}\right){=\int }_{L}{\mu }_{E}\;\mathrm{d}l$ [31, p 2] where I_A is the photon intensity measured at the detector, I₀ is the initial source intensity and μ_E is the attenuation coefficient at energy E. Here L is a line through S and D_A, with arc measure dl. Hence the transmission data determines a set of integrals of μ_E over lines, and the problem of reconstructing μ_E is equivalent to inversion of the line Radon transform with limited data (e.g., [31, 33]). Note that we need not account for the energy dependence of μ_E in this case as the detectors are energy-resolved, and hence there are no issues due to beam-hardening. See figure 1.

When the attenuation of the incoming and scattered rays is ignored, the Compton scattered intensity in two-dimensions can be modelled as integrals of n_e over toric sections [35, 51, 55]

$$\begin{equation}{I}_{\mathrm{C}}={\int }_{T}{n}_{\mathrm{e}}\;\mathrm{d}t,\end{equation} \tag{ 1.1 }$$

where I_C is the Compton scattered intensity measured at a point on D_C. A toric section T = C₁ ∪ C₂ is the union of two intersecting circles of the same radii (as displayed in figure 2), and dt is the arc measure on T. The recovery of n_e is equivalent to inversion of the toric section Radon transform [34, 35, 51, 55]. See figure 2. See also [41, 43] for alternative reconstruction methods. We now discuss the approximation made above to neglect the attenuative effects from the CST model. When the attenuation effects are included, the inverse scattering problem becomes nonlinear [42]. We choose to focus on the analysis of the idealised linear case here, as this allows us to apply the well established theory on linear Fourier integral operators (FIO) and microlocal analysis to derive expression for the image artefacts. Such analysis will likely give valuable insight into the expected artefacts in the nonlinear case. The nonlinear models and their inversion properties are left for future work.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The line and circular arc Radon transforms with full data, are known [30, 35] to have inverses that are continuous in some range of Sobolev norms. Hence with adequate regularization we can reconstruct an image free of artefacts. With limited data however [31, 51], the solution is unstable and the image wavefront set (see definition 2.2) is not recovered stably in all directions. We will see later in section 3 through simulation that such data limitations in the parallel line geometry cause a blurring artefact over a cone in the reconstruction. There may also be nonlocal artefacts specific to the geometry (as in [55]), which we shall discover later in section 3 in the geometry of figure 2.

The main goal of this paper is to combine limited datasets in x-ray CT and CST with new lambda tomography regularization techniques, to recover the image edges stably in all directions. We focus particularly on the geometry of figure 1. In lambda tomography the image reconstruction is carried out by filtered backprojection of the Radon projections, where the filter is chosen to emphasize boundaries. This means that the jump singularities in the lambda reconstruction have the same location and direction to those of the target function, but the smooth parts are undetermined. A common choice of filter is a second derivative in the linear variable [7, 43]. The application of the derivative filter emphasizes the singularities in the Radon projections, and this is a key idea behind lambda tomography [13, 14, 52], and the microlocal view on lambda CT (e.g., [7, 39, 43]). The regularization penalty we propose aims to minimize the difference ${\Vert}\frac{{\mathrm{d}}^{m}}{\mathrm{d}{s}^{m}}R\left({\mu }_{E}-{n}_{\mathrm{e}}\right){{\Vert}}_{{L}^{2}\left(\mathbb{R}{\times}{S}^{1}\right)}$ for some m ⩾ 1, where R denotes the Radon line transform. Therefore, with a full set of Radon projections, the lambda penalties enforce a similarity in the locations and direction of the image singularities (edges) of μ_E and n_e. Further $\frac{{\mathrm{d}}^{m}}{\mathrm{d}{s}^{m}}R$ for m ⩾ 1 is equivalent to taking m − 1/2 derivatives of the object (this operation is continuous of positive order m − 1/2 in Sobolev scales), and hence its inverse is a smoothing operation, which we expect to be of aid in combatting the measurement noise. In addition, the regularized inverse problem we propose is linear (similarly to the Tikhonov regularized inverse [21, p 99]), which (among other benefits of linearity) allows for the fast application of iterative least squares solvers in the solution.

The literature considers joint image reconstruction and regularization in for example, [1, 4, 6, 10–12, 19, 20, 40, 45–48, 50]. See also the special issue [3] for a more general review of joint reconstruction techniques. In [12] the authors consider the joint reconstruction from Positron emission tomography (PET) and magnetic resonance imaging (MRI) data and use a parallel level set (PLS) prior for the joint regularization. The PLS approach (first introduced in [10]) imposes soft constraints on the equality of the image gradient location and direction, thus enforcing structural similarity in the image wavefront sets. This follows a similar intuition to the 'Nambu' functionals of [48] and the 'cross-gradient' methods of [16, 17] in seismic imaging, the latter of which specify hard constraints that the gradient cross products are zero (i.e. parallel image gradients). The methods of [12] use linear and quadratic formulations of PLS, denoted by linear PLS (LPLS) and quadratic PLS (QPLS). The LPLS method will be a point of comparison with the proposed method. We choose to compare with LPLS as it is shown to offer greater performance than QPLS in the experiments conducted in [12].

In [20] the authors consider a class of techniques in joint reconstruction and regularization, including inversion through correspondence mapping, mutual information and joint total variation (JTV). In addition to LPLS, we will compare against JTV as the intuition of JTV is similar to that of lambda regularization (and LPLS), in the sense that a structural similarity is enforced in the image wavefronts. Similar to standard total variation (TV), which favours sparsity in the (single) image gradient, the JTV penalties (first introduced in [45] for colour imaging) favour sparsity in the joint gradient. Thus the image gradients are more likely to occur in the same location and direction upon minimization of JTV. The JTV penalties also have generalizations in colour imaging and vector-valued imaging [4].

In [54] the authors introduce a new toric section transform $\mathcal{T}$ in the geometry of figure 2. Here explicit inversion formulae are derived, but the stability analysis is lacking. We aim to address the stability of $\mathcal{T}$ in this work from a microlocal perspective. Through an analysis of the canonical relations of $\mathcal{T}$, we discover the existence of nonlocal artefacts in the inversion, similarly to [55]. In [40] the joint reconstruction of μ_E and n_e is considered in a pencil beam scanner geometry. Here gradient descent solvers are applied to nonlinear objectives, derived from the physical models, and a weighted, iterative Tikhonov type penalty is applied. The works of [6] improve the wavefront set recovery in limited angle CT using a partially learned, hybrid reconstruction scheme, which adopts ideas in microlocal analysis and neural networks. The fusion with Compton data is not considered however. In our work we assume an equality in the wavefront sets of n_e and μ_E (in a similar vein to [47]), and we investigate the microlocal advantages of combining Compton and transmission data, as such an analysis is lacking in the literature.

The remainder of this paper is organized as follows. In section 2, we recall some definitions and theorems from microlocal analysis. In section 3 we present a microlocal analysis of $\mathcal{T}$ and explain the image artefacts in the n_e reconstruction. Here we prove our main theorem (theorem 3.2), where we show that the canonical relation $\mathcal{C}$ of $\mathcal{T}$ is 2–1. This implies the existence of nonlocal image artefacts in a reconstruction from toric section integral data. Further we find explicit expressions for the nonlocal artefacts and simulate these by applying the normal operations ${\mathcal{T}}^{{\ast}}\mathcal{T}$ to a delta function. In section 4 we consider the microlocal artefacts from x-ray (transmission) data. This yields a limited dataset for the Radon transform, whereby we have knowledge the line integrals for all L which intersect S and D_A (see figure 1). We use the results in [5] to describe the resulting artefacts in the x-ray CT reconstruction. In section 5, we detail our joint reconstruction method for the simultaneous reconstruction of μ_E and n_e. Later in section 5.4 we present simulated reconstructions of μ_E and n_e using the proposed methods and compare against JTV [20] and LPLS [12] from the literature. We also give a comparison to a separate reconstruction using TV.

We next provide some notation and definitions. Let X and Y be open subsets of ${\mathbb{R}}^{n}$. Let $\mathcal{D}\left(X\right)$ be the space of smooth functions compactly supported on X with the standard topology and let ${\mathcal{D}}^{\prime }\left(X\right)$ denote its dual space, the vector space of distributions on X. Let $\mathcal{E}\left(X\right)$ be the space of all smooth functions on X with the standard topology and let ${\mathcal{E}}^{\prime }\left(X\right)$ denote its dual space, the vector space of distributions with compact support contained in X. Finally, let $\mathcal{S}\left({\mathbb{R}}^{n}\right)$ be the space of Schwartz functions, that are rapidly decreasing at ∞ along with all derivatives. See [44] for more information.

Definition 2.1. ([25, definition 7.1.1]). For a function f in the Schwartz space $\mathcal{S}\left({\mathbb{R}}^{n}\right)$, we define the Fourier transform and its inverse as

$$\begin{equation}\mathcal{F}f\left(\boldsymbol{\xi }\right)={\int }_{{\mathbb{R}}^{n}}{\mathrm{e}}^{-\mathrm{i}\mathbf{\text{x}}\cdot \boldsymbol{\xi }}f\left(\mathbf{\text{x}}\right)\mathrm{d}\mathbf{\text{x}},\quad {\mathcal{F}}^{-1}f\left(\mathbf{\text{x}}\right)={\left(2\pi \right)}^{-n}{\int }_{{\mathbb{R}}^{n}}{\mathrm{e}}^{\mathrm{i}\mathbf{\text{x}}\cdot \boldsymbol{\xi }}f\left(\boldsymbol{\xi }\right)\mathrm{d}\boldsymbol{\xi }.\end{equation} \tag{ 2.1 }$$

We use the standard multi-index notation: if $\alpha =\left({\alpha }_{1},{\alpha }_{2},\dots ,\;{\alpha }_{n}\right)\in {\left\{0,1,2,\dots \;\right\}}^{n}$ is a multi-index and f is a function on ${\mathbb{R}}^{n}$, then

$$\begin{equation*}{\partial }^{\alpha }f={\left(\frac{\partial }{\partial {x}_{1}}\right)}^{{\alpha }_{1}}{\left(\frac{\partial }{\partial {x}_{2}}\right)}^{{\alpha }_{2}}\cdots {\left(\frac{\partial }{\partial {x}_{n}}\right)}^{{\alpha }_{n}}f.\end{equation*}$$

If f is a function of (y, x, σ) then ${\partial }_{\mathbf{\text{y}}}^{\alpha }f$ and ${\partial }_{\boldsymbol{\sigma }}^{\alpha }f$ are defined similarly.

We identify cotangent spaces on Euclidean spaces with the underlying Euclidean spaces, so we identify T*(X) with $X{\times}{\mathbb{R}}^{n}$.

If ϕ is a function of $\left(\mathbf{\text{y}},\mathbf{\text{x}},\boldsymbol{\sigma }\right)\in Y{\times}X{\times}{\mathbb{R}}^{N}$ then we define ${\mathrm{d}}_{\mathbf{\text{y}}}\phi =\left(\frac{\partial \phi }{\partial {y}_{1}},\frac{\partial \phi }{\partial {y}_{2}},\dots ,\;\frac{\partial \phi }{\partial {y}_{n}}\right)$, and d_xϕ and d_σϕ are defined similarly. We let $\mathrm{d}\phi =\left({\mathrm{d}}_{\mathbf{\text{y}}}\phi ,{\mathrm{d}}_{\mathbf{\text{x}}}\phi ,{\mathrm{d}}_{\boldsymbol{\sigma }}\phi \right)$.

The singularities of a function and the directions in which they occur are described by the wavefront set [9, p 16]:

Definition 2.2. Let X an open subset of ${\mathbb{R}}^{n}$ and let f be a distribution in ${\mathcal{D}}^{\prime }\left(X\right)$. Let $\left({\mathbf{\text{x}}}_{0},{\boldsymbol{\xi }}_{0}\right)\in X{\times}\left({\mathbb{R}}^{n}{\backslash}\left\{\mathbf{\text{0}}\right\}\right)$. Then f is smooth at x₀ in direction ξ₀ if exists a neighbourhood U of x₀ and V of ξ₀ such that for every $\phi \in \mathcal{D}\left(U\right)$ and $N\in \mathbb{R}$ there exists a constant C_N such that

$$\begin{equation}\left\vert \mathcal{F}\left(\phi f\right)\left(\lambda \boldsymbol{\xi }\right)\right\vert {\leqslant}{C}_{N}{\left(1+\left\vert \lambda \right\vert \right)}^{-N}.\end{equation} \tag{ 2.2 }$$

The pair (x₀, ξ₀) is in the wavefront set, WF(f), if f is not smooth at x₀ in direction ξ₀.

This definition follows the intuitive idea that the elements of WF(f) are the point–normal vector pairs above points of X where f has singularities. For example, if f is the characteristic function of the unit ball in ${\mathbb{R}}^{3}$, then its wavefront set is WF(f) = {(x, tx) : x ∈ S², t ≠ 0}, the set of points on a sphere paired with the corresponding normal vectors to the sphere.

The wavefront set of a distribution on X is normally defined as a subset the cotangent bundle T*(X) so it is invariant under diffeomorphisms, but we will continue to identify ${T}^{{\ast}}\left(X\right)=X{\times}{\mathbb{R}}^{n}$ and consider WF(f) as a subset of $X{\times}{\mathbb{R}}^{n}{\backslash}\left\{\mathbf{\text{0}}\right\}$.

Definition 2.3. ([25, definition 7.8.1]). We define ${S}^{m}\left(Y{\times}X{\times}{\mathbb{R}}^{N}\right)$ to be the set of $a\in \mathcal{E}\left(Y{\times}X{\times}{\mathbb{R}}^{N}\right)$ such that for every compact set K ⊂ Y × X and all multi-indices α, β, γ the bound

$$\begin{equation*}\left\vert {\partial }_{\mathbf{\text{y}}}^{\gamma }{\partial }_{\mathbf{\text{x}}}^{\beta }{\partial }_{\boldsymbol{\sigma }}^{\alpha }a\left(\mathbf{\text{y}},\mathbf{\text{x}},\boldsymbol{\sigma }\right)\right\vert {\leqslant}{C}_{K,\alpha ,\beta }{\left(1+\vert \boldsymbol{\sigma }\vert \right)}^{m-\vert \alpha \vert },\quad \left(\mathbf{\text{y}},\mathbf{\text{x}}\right)\in K,\;\boldsymbol{\sigma }\in {\mathbb{R}}^{N},\end{equation*}$$

holds for some constant C_K,α,β > 0. The elements of S^m are called symbols of order m.

Note that these symbols are sometimes denoted ${S}_{1,0}^{m}$.

Definition 2.4. ([26, definition 21.2.15]). A function $\phi =\phi \left(\mathbf{\text{y}},\mathbf{\text{x}},\boldsymbol{\sigma }\right)\in \mathcal{E}\left(Y{\times}X{\times}{\mathbb{R}}^{N}{\backslash}\left\{\mathbf{\text{0}}\right\}\right)$ is a phase function if ϕ(y, x, λσ) = λϕ(y, x, σ), ∀λ > 0 and dϕ is nowhere zero. A phase function is clean if the critical set Σ_ϕ = {(y, x, σ) : d_σϕ(y, x, σ) = 0} is a smooth manifold with tangent space defined by $\mathrm{d}\left({\mathrm{d}}_{\boldsymbol{\sigma }}\phi \right)=0$.

By the implicit function theorem the requirement for a phase function to be clean is satisfied if $\mathrm{d}\left({\mathrm{d}}_{\boldsymbol{\sigma }}\phi \right)$ has constant rank.

Definition 2.5. ([26, definition 21.2.15] and [27, section 25.2]). Let X and Y be open subsets of ${\mathbb{R}}^{n}$. Let $\phi \in \mathcal{E}\left(Y{\times}X{\times}{\mathbb{R}}^{N}\right)$ be a clean phase function. In addition, we assume that ϕ is nondegenerate in the following sense:

$$\begin{equation*}{\mathrm{d}}_{\mathbf{\text{y}},\boldsymbol{\sigma }}\phi \;\text{and}\;{\mathrm{d}}_{\mathbf{\text{x}},\boldsymbol{\sigma }}\phi \;\text{are}\;\text{never}\;\text{zero.}\end{equation*}$$

The critical set of ϕ is

$$\begin{equation*}{{\Sigma}}_{\phi }=\left\{\left(\mathbf{\text{y}},\mathbf{\text{x}},\boldsymbol{\sigma }\right)\in Y{\times}X{\times}{\mathbb{R}}^{N}{\backslash}\left\{\mathbf{\text{0}}\right\}:{\mathrm{d}}_{\boldsymbol{\sigma }}\phi =0\right\}.\end{equation*}$$

The canonical relation parametrised by ϕ is defined as

$$\begin{equation}\mathcal{C}=\left\{\left(\left(\mathbf{\text{y}},{\mathrm{d}}_{\mathbf{\text{y}}}\phi \left(\mathbf{\text{y}},\mathbf{\text{x}},\boldsymbol{\sigma }\right)\right);\left(\mathbf{\text{x}},-{\mathrm{d}}_{\mathbf{\text{x}}}\phi \left(\mathbf{\text{y}},\mathbf{\text{x}},\boldsymbol{\sigma }\right)\right)\right):\left(\mathbf{\text{y}},\mathbf{\text{x}},\boldsymbol{\sigma }\right)\in {{\Sigma}}_{\phi }\right\},\end{equation} \tag{ 2.3 }$$

Definition 2.6. Let X and Y be open subsets of ${\mathbb{R}}^{n}$. A Fourier integral operator (FIO) of order m + N/2 − n/2 is an operator $A:\mathcal{D}\left(X\right)\to {\mathcal{D}}^{\prime }\left(Y\right)$ with Schwartz kernel given by an oscillatory integral of the form

$$\begin{equation}{K}_{A}\left(\mathbf{\text{y}},\mathbf{\text{x}}\right)={\int }_{{\mathbb{R}}^{N}}{\mathrm{e}}^{\mathrm{i}\phi \left(\mathbf{\text{y}},\mathbf{\text{x}},\boldsymbol{\sigma }\right)}a\left(\mathbf{\text{y}},\mathbf{\text{x}},\boldsymbol{\sigma }\right)\;\mathrm{d}\boldsymbol{\sigma },\end{equation} \tag{ 2.4 }$$

where ϕ is a clean nondegenerate phase function and $a\in {S}^{m}\left(Y{\times}X{\times}{\mathbb{R}}^{N}\right)$ is a symbol. The canonical relation of A is the canonical relation of ϕ defined in (2.3).

This is a simplified version of the definition of FIO in [8, section 2.4] or [27, section 25.2] that is suitable for our purposes since our phase functions are global. For general information about FIOs see [8, 26, 27].

Definition 2.7. Let $\mathcal{C}\subset {T}^{{\ast}}\left(Y{\times}X\right)$ be the canonical relation associated to the FIO $A:{\mathcal{E}}^{\prime }\left(X\right)\to {\mathcal{D}}^{\prime }\left(Y\right)$. Then we let π_L and π_R denote the natural left- and right-projections of $\mathcal{C}$, ${\pi }_{\mathrm{L}}:\mathcal{C}\to {T}^{{\ast}}\left(Y\right)$ and ${\pi }_{\mathrm{R}}:\mathcal{C}\to {T}^{{\ast}}\left(X\right)$.

Because ϕ is nondegenerate, the projections do not map to the zero section. If an FIO $\mathcal{F}$ satisfies our next definition, then ${\mathcal{F}}^{{\ast}}\mathcal{F}$ (or ${\mathcal{F}}^{{\ast}}\phi \mathcal{F}$ if $\mathcal{F}$ does not map to ${\mathcal{E}}^{\prime }\left(Y\right)$) is a pseudodifferential operator [18, 37].

Definition 2.8. Let $\mathcal{F}:{\mathcal{E}}^{\prime }\left(X\right)\to {\mathcal{D}}^{\prime }\left(Y\right)$ be an FIO with canonical relation $\mathcal{C}$ then $\mathcal{F}$ (or $\mathcal{C}$) satisfies the semi-global Bolker assumption if the natural projection ${\pi }_{Y}:\mathcal{C}\to {T}^{{\ast}}\left(Y\right)$ is an embedding (injective immersion).

Here we present a microlocal analysis of the toric section transform in the translational (parallel line) scanning geometry. Through an analysis of two separate limited data problems for the circle transform (where the integrals over circles with centres on a straight line are known) and using microlocal analysis, we show that the canonical relation of the toric section transform is 2–1. The analysis follows in a similar way to the work of [55]. We discuss the nonlocal artefacts inherent to the toric section inversion in section 3.1, and then go on to explain the artefacts due to limited data in section 3.2.

We first define our geometry and formulate the toric section transform of [54] in terms of δ functions, before proving our main microlocal theory.

Let r_m > 1 and define the set of points to be scanned as

$$\begin{equation*}X{:=}\left\{\left({x}_{1},{x}_{2}\right)\in {\mathbb{R}}^{2}:2-{r}_{m}{< }{x}_{2}{< }1\right\}.\end{equation*}$$

Note that r_m controls the depth of the scanning tunnel as in figures 1 and 2. Let

$$\begin{equation*}Y{:=}\left(0,\infty \right){\times}\mathbb{R}\end{equation*}$$

then for j = 1, 2, and (s, x₀) ∈ Y, we define the circles C_j and their centres c_j

$$\begin{equation}\begin{aligned}\hfill & r=\sqrt{{s}^{2}+1},\quad {\mathbf{\text{c}}}_{j}\left(s,{x}_{0}\right)=\left({\left(-1\right)}^{j}s+{x}_{0},2\right)\hfill \\ \hfill & {C}_{j}\left(s,{x}_{0}\right)=\left\{\mathbf{\text{x}}\in {\mathbb{R}}^{2}:\vert \mathbf{\text{x}}-{\mathbf{\text{c}}}_{j}\left(s,{x}_{0}\right){\vert }^{2}-{s}^{2}-1=0\right\}.\hfill \end{aligned}\end{equation} \tag{ 3.1 }$$

Note that $r=\sqrt{{s}^{2}+1}$ is the radius of the circle C_j. The union of the reflected circles C₁ ∪ C₂ is called a toric section. Let $f\in {L}_{0}^{2}\left(X\right)$ be the electron charge density. To define the toric section transform we first introduce two circle transforms

$$\begin{equation}{\mathcal{T}}_{1}f\left(s,{x}_{0}\right)={\int }_{{C}_{1}}f\mathrm{d}s,\quad {\mathcal{T}}_{2}f\left(s,{x}_{0}\right)={\int }_{{C}_{2}}f\mathrm{d}s.\end{equation} \tag{ 3.2 }$$

Now we have the definition of the toric section transform [54]

$$\begin{equation}\mathcal{T}f\left(s,{x}_{0}\right)={\int }_{{C}_{1}\cup {C}_{2}}f\mathrm{d}s={\mathcal{T}}_{1}\left(f\right)\left(s,{x}_{0}\right)+{\mathcal{T}}_{2}\left(f\right)\left(s,{x}_{0}\right)\end{equation} \tag{ 3.3 }$$

where ds denotes the arc element on a circle and (s, x₀) ∈ Y.

We express $\mathcal{T}$ in terms of delta functions as is done for the generalized Funk–Radon transforms studied by Palamodov [36].

$$\begin{align}\hfill \mathcal{T}f\left(s,{x}_{0}\right)& ={\mathcal{T}}_{1}f\left(s,{x}_{0}\right)+{\mathcal{T}}_{2}f\left(s,{x}_{0}\right)=2r\sum _{j=1}^{2}{\int }_{{\mathbb{R}}^{2}}\delta \left(\vert \mathbf{\text{x}}-{\mathbf{\text{c}}}_{j}\left(s,{x}_{0}\right){\vert }^{2}-{s}^{2}-1\right)f\left(\mathbf{\text{x}}\right)\mathrm{d}\mathbf{\text{x}}\hfill \\ \hfill & =2r\sum _{j=1}^{2}{\int }_{-\infty }^{\infty }{\int }_{{\mathbb{R}}^{2}}{\mathrm{e}}^{-\mathrm{i}\sigma \left(\vert \mathbf{\text{x}}-\left({\left(-1\right)}^{j}s+{x}_{0},2\right){\vert }^{2}-{s}^{2}-1\right)}f\left(\mathbf{\text{x}}\right)\mathrm{d}\mathbf{\text{x}}\,\mathrm{d}\sigma .\hfill \end{align} \tag{ 3.4 }$$

Note that the factor in front of the integrals comes about using the change of variables formula and that ${\mathcal{T}}_{j}f=\int \delta \left(\vert \mathbf{\text{x}}-{\mathbf{\text{c}}}_{j}\left(s,{x}_{0}\right)\vert -\sqrt{{s}^{2}+1}\right)f\left(x\right)\;\mathrm{d}\mathbf{\text{x}}$. So the toric section transform is the sum of two FIO's with phase functions

$$\begin{equation*}{\phi }_{j}\left(s,{x}_{0},\mathbf{\text{x}},\sigma \right)=\sigma \left(\vert \mathbf{\text{x}}-\left({\left(-1\right)}^{j}s+{x}_{0},2\right){\vert }^{2}-{s}^{2}-1\right)\end{equation*}$$

for j = 1, 2. Our distributions f are supported away from the intersection points of C₁ and C₂, and hence we can consider the microlocal properties of ${\mathcal{T}}_{1}$ and ${\mathcal{T}}_{2}$ separately to describe the microlocal properties of $\mathcal{T}$.

Proposition 3.1. For j = 1, 2, the circle transform ${\mathcal{T}}_{j}$ is an FIO or order −1/2 with canonical relation

$$\begin{align}\hfill {\mathcal{C}}_{j}=& \left\{\left(\left(s,{x}_{0},{\left(-1\right)}^{j-1}\sigma \left({x}_{1}-{x}_{0}\right),-\sigma \left({\left(-1\right)}^{j-1}s+{x}_{1}-{x}_{0}\right)\right);\left(\mathbf{\text{x}},-\sigma \left(\mathbf{\text{x}}-{\mathbf{\text{c}}}_{j}\left(s,{x}_{0}\right)\right)\right)\right):\right.\hfill \\ \hfill & \;\;\left.\left(s,{x}_{0}\right)\in Y,\sigma \in \mathbb{R}{\backslash}\left\{\mathbf{\text{0}}\right\},\mathbf{\text{x}}\in {C}_{j}\left(s,{x}_{0}\right)\cap \left\{{x}_{2}{< }1\right\}\right\}.\hfill \end{align} \tag{ 3.5 }$$

Furthermore ${\mathcal{C}}_{j}$ satisfies the semi-global Bolker assumption for j = 1, 2.

Proof. First, one can check that ϕ_j and ${\mathcal{T}}_{j}$ both satisfy the restrictions in definition 2.6 so ${\mathcal{T}}_{j}$ is an FIO. Using this definition again and the fact that its symbol is order zero [37], one sees that it has order −1/2.

A straightforward calculation using definition 2.5 shows that the canonical relation of ${\mathcal{T}}_{j}$ is as given in (3.5). Note that we have absorbed a factor of 2 into σ in this calculation. Global coordinates on ${\mathcal{C}}_{j}$ are given by

$$\begin{align*}\hfill \left(s,{x}_{0},{x}_{1},\sigma \right){\mapsto}& \left(s,{x}_{0},{\left(-1\right)}^{j-1}\sigma \left({x}_{1}-{x}_{0}\right),-\sigma \left({\left(-1\right)}^{j-1}s+{x}_{1}-{x}_{0}\right);\right.\hfill \\ \hfill & \;\;\left.\left({x}_{1},{x}_{2}\right),-\sigma \left(\left({x}_{1},{x}_{2}\right)-{\mathbf{\text{c}}}_{j}\left(s,{x}_{0}\right)\right)\right)\hfill \end{align*}$$

where

$$\begin{equation}{x}_{2}=2-\sqrt{{s}^{2}+1-{\left({x}_{1}-\left({x}_{0}+{\left(-1\right)}^{j}s\right)\right)}^{2}}\end{equation} \tag{ 3.6 }$$

because x₂ < 1. Recall that c_j is given in (3.1).

We now show that ${\mathcal{C}}_{j}$ satisfies the semiglobal Bolker assumption by finding a smooth inverse in these coordinates to the projection ${{\Pi}}_{\mathrm{L}}:{\mathcal{C}}_{j}\to {T}^{{\ast}}\left(Y\right)$. Let $\lambda =\left(s,{x}_{0},{\tau }_{1},{\tau }_{2}\right)\in {{\Pi}}_{\mathrm{L}}\left({\mathcal{C}}_{j}\right)$. We solve for x₁ and σ in the equation Π_L(s, x₀, x₁, σ) = λ. Then, s and x₀ are known as are

$$\begin{equation}{\tau }_{1}={\left(-1\right)}^{j-1}\sigma \left({x}_{1}-{x}_{0}\right)\quad {\tau }_{2}=-\sigma \left({\left(-1\right)}^{j-1}s+{x}_{1}-{x}_{0}\right).\end{equation} \tag{ 3.7 }$$

A straightforward linear algebra exercise shows that the unique solutions for σ and x₁ are

$$\begin{equation}\sigma =\frac{{\left(-1\right)}^{j}{\tau }_{2}-{\tau }_{1}}{s},\quad {x}_{1}=\frac{s{\tau }_{1}}{{\left(-1\right)}^{j}{\tau }_{1}-{\tau }_{2}}+{x}_{0}\end{equation} \tag{ 3.8 }$$

This gives a smooth inverse to Π_L on the image ${{\Pi}}_{\mathrm{L}}\left({\mathcal{C}}_{j}\right)$ and finishes the proof. □

Because ${\mathcal{C}}_{j}$ satisfies the Bolker assumption, the composition ${\mathcal{C}}_{j}^{{\ast}}{\circ}{\mathcal{C}}_{j}\subset {\Delta}$, where Δ is the diagonal in T*(X). Hence in a reconstruction from circular integral data with centres on a line we would not expect to see image artefacts for functions supported in x − 2 > 0 unless one uses a sharp cutoff on the data.

The canonical relation $\mathcal{C}$ of $\mathcal{T}$ can be written as the disjoint union $\mathcal{C}={\mathcal{C}}_{1}\cup {\mathcal{C}}_{2}$ since (C₁(s, x₀) ∩ C₂(s, x₀)) ∩ supp(f) = ∅ for any (s, x₀) ∈ Y.

For convenience, we will sometimes label the coordinate x₀ in (3.6) as ${\left({x}_{0}\right)}_{1}$ it is associated with ${\mathcal{C}}_{1}$ and ${\left({x}_{0}\right)}_{2}$ if it is associated with ${\mathcal{C}}_{2}$.

Theorem 3.2. For j = 1, 2, the projection ${\pi }_{\mathrm{R}}:{\mathcal{C}}_{j}\to {T}^{{\ast}}\left(X\right)$ is bijective onto the set

$$\begin{equation}D=\left\{\left(\mathbf{\text{x}},\boldsymbol{\xi }\right)\in {T}^{{\ast}}\left(X\right):{\xi }_{2}\ne 0\right\}.\end{equation} \tag{ 3.9 }$$

In addition, ${\pi }_{\mathrm{R}}:\mathcal{C}\to {T}^{{\ast}}\left(X\right)$ is two-to one onto D.

Proof. Let $\mu =\left(\mathbf{\text{x}};\boldsymbol{\xi }\right)\in {T}^{{\ast}}\left(X\right){\backslash}\left\{\mathbf{\text{0}}\right\}$ and let x = (x₁, x₂) and ξ = (ξ₁, ξ₂). If $\mu \in {\pi }_{\mathrm{R}}\left({\mathcal{C}}_{j}\right)$ for either j = 1 or j = 2, then ξ₂ ≠ 0 by (3.5) since x₂ < 2. For the rest of the proof, assume μ is in the set D given by (3.9).

We will now describe the preimage of μ in ${\mathcal{C}}_{j}$. The covector μ is conormal to a unique circle centred on x₂ = 2, and its centre is on the line through x and parallel ξ. If the centre has coordinates (c, 2), then a calculation shows that c is given by

$$\begin{equation}c=c\left(\mathbf{\text{x}},\boldsymbol{\xi }\right)={x}_{1}-\frac{{\xi }_{1}\left({x}_{2}-2\right)}{{\xi }_{2}}.\end{equation} \tag{ 3.10 }$$

Using this calculation, one sees that the radius of the circle and coordinate s are given by

$$\begin{equation}r=r\left(\mathbf{\text{x}},\boldsymbol{\xi }\right)=\frac{\left(2-{x}_{2}\right)\left\vert \boldsymbol{\xi }\right\vert }{\left\vert {\xi }_{2}\right\vert },\quad s=s\left(\mathbf{\text{x}},\boldsymbol{\xi }\right)=\sqrt{{r}^{2}-1}\end{equation} \tag{ 3.11 }$$

and the coordinate ${\left({x}_{0}\right)}_{j}$ is given by

$$\begin{equation}{\left({x}_{0}\right)}_{j}={\left({x}_{0}\right)}_{j}\left(\mathbf{\text{x}},\boldsymbol{\xi }\right)={x}_{1}+\frac{{\xi }_{1}\left(2-{x}_{2}\right)}{{\xi }_{2}}+{\left(-1\right)}^{j-1}s\quad \text{for} j=1,2.\end{equation} \tag{ 3.12 }$$

A straightforward calculation shows that

$$\begin{equation}\sigma =\sigma \left(\mathbf{\text{x}},\boldsymbol{\xi }\right)=\frac{-{\xi }_{2}}{2-{x}_{2}}.\end{equation} \tag{ 3.13 }$$

This gives the coordinates (3.6) on ${\mathcal{C}}_{j}$ and shows that ${\pi }_{\mathrm{R}}:{\mathcal{C}}_{j}\to D$ is injective with smooth inverse.

Now, we consider the projection from $\mathcal{C}$. Given (x, ξ) ∈ D, our calculations show that the preimage in $\mathcal{C}$, in coordinates (3.6) is given by two distinct points

$$\begin{equation*}\left(s\left(\mathbf{\text{x}},\boldsymbol{\xi }\right),{\left({x}_{0}\right)}_{j}\left(\mathbf{\text{x}},\boldsymbol{\xi }\right),{x}_{1},\sigma \left(\mathbf{\text{x}},\boldsymbol{\xi }\right)\right)\quad \text{for} j=1,2.\end{equation*}$$

The coordinates are given by (3.11)–(3.13) respectively. □

The abstract adjoint ${\mathcal{T}}_{j}^{t}$ cannot be composed with ${\mathcal{T}}_{i}$ for i = 1, 2, because the support of ${\mathcal{T}}_{i}f$ can be unbounded in r, even for $f\in {\mathcal{E}}^{\prime }\left(X\right)$ and ${\mathcal{T}}_{j}^{t}$ is not defined for such distributions. Therefore, we introduce a smooth cutoff function. Choose r_M > 2 and let ψ(s) be a smooth compactly supported function equal to one for $s\in \left[1,\sqrt{1-{r}_{M}^{2}}\right]$ and define

$$\begin{equation}{\mathcal{T}}_{j}^{{\ast}}g={\mathcal{T}}_{j}^{t}\left(\psi g\right)\end{equation} \tag{ 3.14 }$$

for all $g\in {\mathcal{D}}^{\prime }\left(Y\right)$ because our bound on r introduces a bound on x₀ so the integral is over a bounded set for each x ∈ X.

3.1. The nonlocal artefacts

Now, we can state our next theorem, which describes the artefacts that can be added to the reconstruction using the normal operator, ${\mathcal{T}}^{{\ast}}\mathcal{T}$.

Theorem 3.3. If $f\in {\mathcal{E}}^{\prime }\left(X\right)$ then

$$\begin{equation}\mathrm{W}\mathrm{F}\left({\mathcal{T}}^{{\ast}}\mathcal{T}f\right)\subset \left(\mathrm{W}\mathrm{F}\left(f\right)\cap D\right)\cup {{\Lambda}}_{12}\left(f\right)\cup {{\Lambda}}_{21}\left(f\right)\end{equation} \tag{ 3.15 }$$

where D is given by (3.9), and the sets Λ_ij are given for (x, ξ) ∈ D by

$$\begin{equation}{{\Lambda}}_{ij}\left(f\right)=\left\{{\lambda }_{ij}\left(\mathbf{\text{x}},\boldsymbol{\xi }\right):\left(\mathbf{\text{x}},\boldsymbol{\xi }\right)\in \mathrm{W}\mathrm{F}\left(f\right)\cap D\right\}\end{equation} \tag{ 3.16 }$$

where the functions λ₁₂ and λ₂₁ are given by (3.19) and (3.20) respectively. Note that the functions λ_ij are defined for only some (x, ξ) ∈ D and singularities at other points do not generate artefacts.

Therefore, ${\mathcal{T}}^{{\ast}}\mathcal{T}$ recovers most singularities of f, as indicated in the first term in (3.15), but it adds two sets of nonlocal singularities, as given by Λ₁₂(f) and Λ₂₁(f). Note that, even if ${\mathcal{T}}_{j}^{{\ast}}$ and ${\mathcal{T}}_{j}$ are both elliptic above a covector (x, ξ), artefacts caused by other points could mask singularities of f that 'should' be visible in ${\mathcal{T}}^{{\ast}}\mathcal{T}f$.

Proof. Let $f\in {\mathcal{E}}^{\prime }\left(X\right)$. By the Hörmander–Sato lemma [25, theorem 8.2.13]. We have the expansion

$$\begin{align}\hfill & \mathrm{W}\mathrm{F}\left({\mathcal{T}}^{{\ast}}\mathcal{T}\left(f\right)\right)\subset \left({\mathcal{C}}^{{\ast}}{\circ}\mathcal{C}\right){\circ}\mathrm{W}\mathrm{F}\left(f\right)\hfill \\ \hfill & \;\;\;\;\quad =\left[\left({\mathcal{C}}_{1}^{{\ast}}{\circ}{\mathcal{C}}_{1}\right)\cup \left({\mathcal{C}}_{2}^{{\ast}}{\circ}{\mathcal{C}}_{2}\right)\right]{\circ}\mathrm{W}\mathrm{F}\left(f\right)\;\cup \;\left({\mathcal{C}}_{2}^{{\ast}}{\circ}{\mathcal{C}}_{1}\right){\circ}\mathrm{W}\mathrm{F}\left(f\right)\;\cup \;\left({\mathcal{C}}_{1}^{{\ast}}{\circ}{\mathcal{C}}_{2}\right){\circ}\mathrm{W}\mathrm{F}\left(f\right)\hfill \end{align} \tag{ 3.17 }$$

The first term in brackets in (3.17) is $\left\{\left(\mathbf{\text{x}},\boldsymbol{\xi };\mathbf{\text{x}},\boldsymbol{\xi }\right):\left(\mathbf{\text{x}},\boldsymbol{\xi }\right)\in D\right\}{\circ}\mathrm{W}\mathrm{F}\left(f\right)=\mathrm{W}\mathrm{F}\left(f\right)\cap D$. This proves the first part of the inclusion (3.15).

We now analyse the other two terms to define the functions λ_ij and finish the proof. Let (x, ξ) ∈ WF(f) ∩ D. First, consider ${\lambda }_{12}\left(\mathbf{\text{x}},\boldsymbol{\xi }\right)={\mathcal{C}}_{2}^{{\ast}}{\circ}{\mathcal{C}}_{1}{\circ}\left(\mathbf{\text{x}},\boldsymbol{\xi }\right)$.³ Using the calculations in the proof of theorem 3.2 one sees that ${\mathcal{C}}_{1}{\circ}\left(\mathbf{\text{x}},\boldsymbol{\xi }\right)$ is given by

$$\begin{equation}\begin{aligned}\hfill & \left(s,{x}_{0},{\tau }_{1},{\tau }_{2}\right)\quad \;\text{where}\;s=s\left(\mathbf{\text{x}},\boldsymbol{\xi }\right)=\sqrt{\frac{{\left(2-{x}_{2}\right)}^{2}{\left\vert \boldsymbol{\xi }\right\vert }^{2}}{{\xi }_{2}^{2}}-1}\hfill \\ \hfill & {x}_{0}={\left({x}_{0}\right)}_{1}\left(\mathbf{\text{x}},\boldsymbol{\xi }\right)={x}_{1}+\frac{{\xi }_{1}\left(2-{x}_{2}\right)}{{\xi }_{2}}+s\quad \sigma =\sigma \left(\mathbf{\text{x}},\boldsymbol{\xi }\right)=\frac{-{\xi }_{2}}{2-{x}_{2}}\hfill \\ \hfill & {\tau }_{1}=\sigma \left({x}_{1}-{x}_{0}\right)\quad {\tau }_{2}=-\sigma \left(s+{x}_{1}-{x}_{0}\right)\hfill \end{aligned}\end{equation} \tag{ 3.18 }$$

where we have taken these from the proof of theorem 3.2. To find ${\mathcal{C}}_{2}^{{\ast}}{\circ}{\mathcal{C}}_{1}{\circ}\left(\mathbf{\text{x}},\boldsymbol{\xi }\right)$ we calculate the composition of the covector described in (3.18) with ${\mathcal{C}}_{2}^{{\ast}}$. Note that the values of x₀ and s are the same in both calculations and are given by (3.18). After using (3.8) and that $\frac{{\xi }_{1}\left(2-{x}_{2}\right)}{{\xi }_{2}}={x}_{0}-s-{x}_{1}$, one sees that

$$\begin{equation}\begin{aligned}\hfill {\lambda }_{12}\left(\mathbf{\text{x}},\boldsymbol{\xi }\right)& =\left(\left({y}_{1},{y}_{2}\right),\boldsymbol{\eta }\right)\quad \text{where}\hfill \\ \hfill {y}_{1}={y}_{1}\left(\mathbf{\text{x}},\boldsymbol{\xi }\right)& =\frac{s\left({x}_{1}-{x}_{0}\right)}{2\left({x}_{1}-{x}_{0}\right)+s}+{x}_{0}\hfill \\ \hfill {y}_{2}={y}_{2}\left(\mathbf{\text{x}},\boldsymbol{\xi }\right)& =2-\sqrt{\frac{{\left(2-{x}_{2}\right)}^{2}{\left\vert \boldsymbol{\xi }\right\vert }^{2}}{{\xi }_{2}^{2}}-{\left({y}_{1}-\left({x}_{0}+s\right)\right)}^{2}}\hfill \\ \hfill \boldsymbol{\eta }& =\left(-2{\xi }_{1}-\frac{s{\xi }_{2}}{2-{x}_{2}}\right)\left(\mathbf{\text{y}}-{\mathbf{\text{c}}}_{2}\left(s\left(\mathbf{\text{x}},\boldsymbol{\xi }\right){\left({x}_{0}\right)}_{1}\left(\mathbf{\text{x}},\boldsymbol{\xi }\right)\right)\right)\hfill \end{aligned}\end{equation} \tag{ 3.19 }$$

where ${x}_{0}={\left({x}_{0}\right)}_{1}\left(\mathbf{\text{x}},\boldsymbol{\xi }\right)$ and s = s(x, ξ) are given in (3.18) and η is calculated using the expression (3.8) with j = 2.

Note that the function λ₁₂ is defined for only some (x, ξ) ∈ D; for example if the argument for the square root defining y₂(x, ξ) is negative, then y₂(x, ξ) is not defined and the point (x, ξ) will not generate artefacts in Λ₁₂.

A similar calculation shows for (y, η) ∈ D that

$$\begin{equation}\begin{aligned}\hfill {\lambda }_{21}\left(\mathbf{\text{y}},\boldsymbol{\eta }\right)& =\left(\left({x}_{1},{x}_{2}\right),\boldsymbol{\xi }\right)\quad \text{where}\;\hfill \\ \hfill {x}_{1}={x}_{1}\left(\mathbf{\text{y}},\boldsymbol{\eta }\right)& =\frac{s\left({y}_{1}-{x}_{0}\right)}{-2\left({y}_{1}-{x}_{0}\right)+s}+{x}_{0}\hfill \\ \hfill {x}_{2}={x}_{2}\left(\mathbf{\text{y}},\boldsymbol{\eta }\right)& =2-\sqrt{{s}^{2}\left(\mathbf{\text{y}},\boldsymbol{\eta }\right)+1-{\left({x}_{1}-\left({x}_{0}-s\right)\right)}^{2}}\hfill \\ \hfill \boldsymbol{\xi }& =\left(2{\eta }_{1}-\frac{s{\eta }_{2}}{2-{y}_{2}}\right)\left(\mathbf{\text{x}}-{\mathbf{\text{c}}}_{1}\left(s\left(\mathbf{\text{y}},\boldsymbol{\eta }\right),{\left({x}_{0}\right)}_{2}\left(\mathbf{\text{y}},\boldsymbol{\eta }\right)\right)\right)\hfill \end{aligned}\end{equation} \tag{ 3.20 }$$

where

$$\begin{equation*}s=s\left(\mathbf{\text{y}},\boldsymbol{\eta }\right)=\sqrt{\frac{{\left(2-{y}_{2}\right)}^{2}{\left\vert \boldsymbol{\eta }\right\vert }^{2}}{{\eta }_{2}^{2}}-1}\quad {x}_{0}={\left({x}_{0}\right)}_{2}\left(\mathbf{\text{y}},\boldsymbol{\eta }\right)={y}_{1}+\frac{{\eta }_{1}\left(2-{y}_{2}\right)}{{\eta }_{2}}-s.\end{equation*}$$

Note that the function λ₂₁ is not defined for all (y, η) ∈ D, and other points (y, η) do not generate artefacts. This is for the same reason as for λ₁₂. □

Remark 3.4. The artefacts caused by a singularity of f are as strong as the reconstruction of that singularity. To see this, first note that each ${\mathcal{T}}_{j}^{{\ast}}{\mathcal{T}}_{i}$ smooths of order one in Sobolev scale since it an FIO of order −1 [24, theorem 4.3.1].

The visible singularities in the reconstruction come from the compositions ${\mathcal{T}}_{1}^{{\ast}}{\mathcal{T}}_{1}$ and ${\mathcal{T}}_{2}^{{\ast}}{\mathcal{T}}_{2}$ since these are pseudodifferential operators of order −1. The artefacts come from the 'cross' compositions ${\mathcal{T}}_{2}^{{\ast}}{\mathcal{T}}_{1}$ and ${\mathcal{T}}_{1}^{{\ast}}{\mathcal{T}}_{2}$, and they are FIO of order −1. Therefore, since the terms that preserve the real singularities of f, ${\mathcal{T}}_{i}^{{\ast}}{\mathcal{T}}_{i}$, i = 1, 2, are also of order −1, ${\mathcal{T}}^{{\ast}}\mathcal{T}$ smooths each singularity of f by one order in Sobolev scale and the composition ${\mathcal{T}}_{2}^{{\ast}}{\mathcal{T}}_{1}$ (corresponding to the artefact λ₁₂, if defined at this covector) can create an artefact from that singularity that are also one order smoother than that singularity, and similarly with the composition ${\mathcal{T}}_{1}^{{\ast}}{\mathcal{T}}_{2}$.

Second, our results are valid, not only for the normal operator ${\mathcal{T}}^{{\ast}}\mathcal{T}$ but for any filtered backprojection method ${\mathcal{T}}^{{\ast}}P\mathcal{T}$ where P is a pseudodifferential operator. This is true since pseudodifferential operators have canonical relation Δ and they do not move singularities, so our microlocal calculations are the same. If P has order k, then ${\mathcal{T}}^{{\ast}}P\mathcal{T}$ decreases the Sobolev order of each singularity of f by order (k − 1) in Sobolev norm and can create an artefact from that singularity of the same order.

3.2. Artefacts for ${\mathcal{T}}^{{\ast}}\mathcal{T}$ due to limited data

In practice we do not have access to $\mathcal{T}f\left(s,{x}_{0}\right)$ for all s ∈ (0, ∞) (or r ∈ (1, ∞)) and ${x}_{0}\in \mathbb{R}$, and will have knowledge of x₀ ∈ (−a, a) and r ∈ (1, r_M) for some a > 0 (see figures 1 and 2) and maximum radius r_M > 1.

We now evaluate which wavefront directions (x, ξ) will be visible from this limited data. Let us consider the pair (x, ξ) ∈ C₂(s, x₀) × S¹ and let β be the angle of ξ from the vertical as depicted in figure 2. Then c₂(s, x₀) = ((2 − x₂)tan β + x₁, 2) and

$$\begin{equation*}\vert \mathbf{\text{x}}-{\mathbf{\text{c}}}_{2}\left(s,{x}_{0}\right){\vert }^{2}={r}^{2}\;\Longrightarrow\;\left(1+{\mathrm{tan}}^{2}\;\beta \right){\left(2-{x}_{2}\right)}^{2}={r}^{2}\;\Longrightarrow\;\mathrm{tan}\;\beta =\sqrt{\frac{{r}^{2}}{{\left(2-{x}_{2}\right)}^{2}}-1}.\end{equation*}$$

Let β_m = β_m(x) ∈ (0, π/2) be defined by

$$\begin{equation}\mathrm{tan}\;{\beta }_{m}=\sqrt{\frac{{r}_{M}^{2}}{{\left(2-{x}_{2}\right)}^{2}}-1}\end{equation} \tag{ 3.21 }$$

(noting that we only consider x such that 1 > x₂ > 2 − r_M). Then the maximum directional coverage of the singularities (wavefront set) at a given x ∈ X which are resolved by the Compton data are described by the open cone of ξ ∈ S¹

$$\begin{equation}{C}_{T}=\left\{{\pm}\left(\mathrm{sin}\;\beta ,-\mathrm{cos}\;\beta \right):-{\beta }_{m}{< }\beta {< }{\beta }_{m}\right\},\end{equation} \tag{ 3.22 }$$

and the opening angle of the cone depends on the depth of x (i.e. x₂). See figure 2. The cone C_T illustrated corresponds to the case when β = β_m.

In all of our numerical experiments, we set the tunnel height as r_m − 1 = 6 and the detector line width is 2a = 8. We let r_M > r_m be large enough to penetrate the entire scanning tunnel (up to the line {x₂ = 2 − r_m} as highlighted in figures 1 and 2), so as to imply a unique reconstruction [55]. Specifically we set the maximum radius r_M = 9 and simulate $\mathcal{T}\left(r,{x}_{0}\right)$ for r ∈ {1 + 0.02j : 1 ⩽ j ⩽ 400} and x₀ ∈ {−4 + 0.04j : 1 ⩽ j ⩽ 200}. Further the densities considered are represented on [−2, 2] × [−3, 1] (200 × 200 pixel grid) in the reconstructions shown. The machine design considered is such that for any x ∈ [−2, 2] × [−1.5, 1] we have the maximal directional coverage in C_T allowed for the limited r < r_M (see figure 7). With the exception of the horizontal bar phantom depicted in figure 14, all objects considered for reconstruction are approximately in this region.

To demonstrate the artefacts, we apply a discrete form of ${\mathcal{T}}^{{\ast}}\mathcal{T}$ to a delta function. We have the expansion

$$\begin{equation}{\mathcal{T}}^{{\ast}}\mathcal{T}\;=\;{\left({\mathcal{T}}_{1}+{\mathcal{T}}_{2}\right)}^{{\ast}}\left({\mathcal{T}}_{1}+{\mathcal{T}}_{2}\right)\;=\;{\mathcal{T}}_{1}^{{\ast}}{\mathcal{T}}_{1}+{\mathcal{T}}_{2}^{{\ast}}{\mathcal{T}}_{2}+{\mathcal{T}}_{1}^{{\ast}}{\mathcal{T}}_{2}+{\mathcal{T}}_{2}^{{\ast}}{\mathcal{T}}_{1}.\end{equation} \tag{ 3.23 }$$

Using equations (3.18)–(3.20) one can show for g ∈ L²(Y) that the backprojection operators ${\mathcal{T}}_{j}^{{\ast}}$, j = 1, 2 can be written

$$\begin{equation}{\mathcal{T}}_{j}^{{\ast}}g\left(\mathbf{\text{x}}\right)={\left.{\int }_{-{\beta }_{m}}^{{\beta }_{m}}g\left(\sqrt{{r}^{2}-1},{x}_{1}+{\left(-1\right)}^{j-1}\sqrt{{r}^{2}-1}+r\;\mathrm{sin}\;\beta \right)\right\vert }_{r=\frac{\left(2-{x}_{2}\right)}{\mathrm{cos}\beta }}\mathrm{d}\beta .\end{equation} \tag{ 3.24 }$$

Note that we are not restricting x₀ to [−a, a] but we are restricting s to $\left(0,\sqrt{{r}_{M}^{2}-1}\right)$, and hence the cutoff function ψ of equation (3.14) is equal to one on the bounds of integration.

Now, let f be a delta function at y. We calculate the artefacts

$$\begin{equation}{\mathcal{T}}_{1}^{{\ast}}{\mathcal{T}}_{2}f\left(\mathbf{\text{x}}\right)\ne 0\;\Longleftrightarrow\;\exists \beta \in \left[-{\beta }_{m},{\beta }_{m}\right]\;\text{s.t.}\;\vert \mathbf{\text{y}}-{\mathbf{\text{c}}}_{2}\left(s,{x}_{0}\right)\vert =r,\end{equation} \tag{ 3.25 }$$

where $r=\frac{2-{x}_{2}}{\mathrm{cos}\;\beta }$, $s=\sqrt{{r}^{2}-1}$ and x₀ = x₁ + s + r sin β. Similarly

$$\begin{equation}{\mathcal{T}}_{2}^{{\ast}}{\mathcal{T}}_{1}f\left(\mathbf{\text{x}}\right)\ne 0\;\Longleftrightarrow\;\exists \beta \in \left[-{\beta }_{m},{\beta }_{m}\right]\;\text{s.t.}\;\vert \mathbf{\text{y}}-{\mathbf{\text{c}}}_{1}\left(s,{x}_{0}\right)\vert =r,\end{equation} \tag{ 3.26 }$$

where $r=\frac{2-{x}_{2}}{\mathrm{cos}\;\beta }$, $s=\sqrt{{r}^{2}-1}$ and x₀ = x₁ − s + r sin β. Hence the only contributions to the backprojection from ${\mathcal{T}}_{1}^{{\ast}}{\mathcal{T}}_{2}$ and ${\mathcal{T}}_{2}^{{\ast}}{\mathcal{T}}_{1}$ are on the following sets:

$$\begin{equation}{S}_{12}=\left\{\mathbf{\text{x}}:\exists \beta \in \left[-{\beta }_{m},{\beta }_{m}\right]\;\text{s.t.}\;\vert \mathbf{\text{y}}-{\mathbf{\text{c}}}_{2}\left(s,{x}_{0}\right)\vert =r\right\}\end{equation} \tag{ 3.27 }$$

where $r=\frac{2-{x}_{2}}{\mathrm{cos}\;\beta }$ and x₀ = x₁ + s + r sin β and

$$\begin{equation}{S}_{21}=\left\{\mathbf{\text{x}}:\exists \beta \in \left[-{\beta }_{m},{\beta }_{m}\right]\;\text{s.t.}\;\vert \mathbf{\text{y}}-{\mathbf{\text{c}}}_{1}\left(s,{x}_{0}\right)\vert =r\right\}.\end{equation} \tag{ 3.28 }$$

where $r=\frac{2-{x}_{2}}{\mathrm{cos}\;\beta }$ and x₀ = x₁ − s + r sin β. This means that all Λ_ij artefacts will be in these sets. Note that besides the Λ_ij artefacts shown in figures 4(e) and (f) there are limited data artefacts caused by circles meeting y of radius r_M (figures 4(a)–(c)) and these are of higher strength in Sobolev norm.

To simulate a δ function discretely we assign a value of one to nine neighbouring pixels in a 200–200 grid (which will represent [−2, 2] × [−3, 1]) and set all other pixel values to zero. Letting our discrete delta function be denoted by x_δ, we approximate ${\mathcal{T}}^{{\ast}}\mathcal{T}\delta \approx {A}^{T}A{x}_{\delta }$, where A is the discrete form of $\mathcal{T}$. For comparison we show images of

$$\begin{equation*}\left({\mathcal{T}}_{1}^{{\ast}}{\mathcal{T}}_{2}+{\mathcal{T}}_{2}^{{\ast}}{\mathcal{T}}_{1}\right)\delta \approx \left({A}_{1}^{T}{A}_{2}+{A}_{2}^{T}{A}_{1}\right){x}_{\delta },\end{equation*}$$

a characteristic function on the set S₁₂ ∪ S₂₁, and the artefacts induced by Λ₁₂ and Λ₂₁. Here A_j is the discrete form of ${\mathcal{T}}_{j}$, for j = 1, 2. See figure 3. For example, in figure 3(b) we see a butterfly wing type artefact in $\left({A}_{1}^{T}{A}_{1}+{A}_{2}^{T}{A}_{2}\right){x}_{\delta }$. This is due to the limited r and x₀ data inherent to our acquisition geometry (there are unresolved wavefront directions). In the $\left({A}_{1}^{T}{A}_{2}+{A}_{1}^{T}{A}_{2}\right){x}_{\delta }$ image of figure 3(c) we see artefacts appearing on the set S₁₂ ∪ S₂₁ as shown in figure 3(d). This is as predicted by our theory. The artefacts induced by the Λ_ij in this case lie outside the scanning region ([−2, 2] × [−3, 1]), and hence they are not observed in the reconstruction. See figures 3(e) and (f). In figure 4 the artefact curves intersect [−2, 2] × [−3, 1] in the top left and right-hand corners respectively. See figures 4(e) and (f). In this case the artefacts are observed faintly in the reconstructions (their magnitude is small compared to the delta function), and it is unclear whether they align with our predictions.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To show the artefacts induced by the Λ_ij more clearly, we repeat the analysis above using filtered backprojection, and a second derivative filter ${\Phi}=\frac{{\mathrm{d}}^{2}}{\mathrm{d}{r}^{2}}$. That is we show images of ${\mathcal{T}}^{{\ast}}{\Phi}\mathcal{T}\delta$. Note that Φ is applied in the variable r (the torus radius). The application of derivative filters is a common idea in lambda tomography [13, 14], and is known to highlight the image contours (singularities or edges) in the reconstruction [43, theorem 3.5]. See figure 5. As the artefacts induced by Λ_ij appeared to be largely outside the scanning region ([−2, 2] × [−3, 1]) in our previous simulations, we have increased the scanning region size to [−3, 3] × [−4, 2], to show more the effects of the Λ_ij in the observed reconstruction. Here Φ suppresses the artefacts due to limited data, and the Λ_ij artefacts appear as additional contours in the reconstruction. The observed artefacts appear most clearly in figures 5(b) and (e), and align exactly with our predictions in figures 5(c) and (f).

Zoom In Zoom Out Reset image size

Remark 3.5. With precise knowledge of the locations of the artefacts induced by the Λ_ij we can assist in the design of the proposed parallel line scanner. That is we can choose a, r_M and the scanning tunnel size to minimize the presence of the nonlocal artefacts in the reconstruction (i.e., those from Λ_ij(f)). Such advice would be of benefit to our industrial partners in airport screening to remove the concern for nonlocal artefacts in the image reconstruction of baggage. Indeed the machine dimensions we have chosen seem to be suitable as the artefacts appear largely outside the reconstruction space (see figures 3 and 4).

The detector row D_C collects Compton (back) scattered data, which determines $\mathcal{T}f\left(s,{x}_{0}\right)$ for a range of s and x₀, where f = n_e is the electron charge density. The forward detector array D_A collects transmission (standard x-ray CT) data, which determine a set of straight line integrals over the attenuation coefficient f = μ_E, for some photon energy E. The data is limited to the set of lines which intersect S (the source array) and D_A. This limited data can cause artefacts in the x-ray reconstruction, and we will analyse these artefacts using the theory in [5]. Let ${L}_{s,\theta }=\left\{\mathbf{\text{x}}\in {\mathbb{R}}^{2}:\mathbf{\text{x}}\cdot {\Theta}=s\right\}$ be the line parameterized by a rotation θ ∈ [0, π] and a directed distance from the origin $s\in \mathbb{R}$. Here Θ = (cos θ, sin θ) and L_s,θ is the line containing sΘ and perpendicular to Θ.

In the scanning geometry of this article, the set S of x-ray transmitters is the segment between (−4, 3) and (4, 3) and the set of x-ray detectors, D_A, is the segment between (−4, −5) and (4, −5) as in figure 1. For this reason, the cutoff in the sinogram space is described by the set

$$\begin{equation}H=\left\{\left(s,\theta \right)\in \mathbb{R}{\times}\left[-\pi /2,\pi /2\right]:{L}_{s,\theta }\cap S\ne \varnothing \quad \text{and}\quad {L}_{s,\theta }\cap {D}_{\mathrm{A}}\ne \varnothing \right\}.\end{equation} \tag{ 4.1 }$$

The characteristic function of H appears as a skewed diamond shape in sinogram space.

To illustrate the added artefacts inherent in this incomplete data problem, we simulate reconstructions of delta functions with transmission CT data on H. That is, we apply the backprojection operator ${R}_{L}^{{\ast}}{R}_{L}$ to δ functions, where R_Lf denotes Rf for (s, θ) ∈ H. See figure 6. By the theory in [5], artefacts caused by the incomplete data occur on lines L_s,θ for (s, θ) in the boundary of H. Each delta function in figure 6 is at a point (t, −1) for some t with 0 < t < 2, so the lines that meet the support of the delta function, (t, −1) that are in the boundary of H must also contain either (4, 3) or (4, −5). This is true because S and D_A are mirror images about the line y = −1 and t ∈ (0, 2).

Zoom In Zoom Out Reset image size

Furthermore, by symmetry (the δ functions are on the centre line of the scanning region), these artefact lines will be reflections of each other in the vertical line x₁ = t. This is illustrated in our reconstructions in figure 6. The opening angle of the cone in the delta reconstructions decreases (fewer wavefront set directions are stably resolved) as we translate δ to the right on the line x₁ = −1.

Example 4.1. We now use these ideas to analyse the visible wavefront directions for the joint problem. Let $\mathcal{S}=\left[-4,4\right]{\times}\left[-5,3\right]$ be the square between S and D_A and let O = (0, −1), the centre of $\mathcal{S}$. We consider wavefronts at points (x₁, x₂) ∈ [−2, 2] × [−3, 1] which is a square centred at O and the region in which our simulated reconstructions are done.

By (4.1), lines in the data set must intersect both S and D_A, so lines through O in the data set are all lines through O that are more vertical than the diagonals of $\mathcal{S}$. Because visible wavefront directions are normal to lines in the x-ray CT data set [38], the wavefront directions which are resolved lie in the horizontal open cone between normals to these diagonals. Therefore, they are in the cone

$$\begin{equation*}{C}_{R}=\left\{{\pm}\left(\mathrm{cos}\;\alpha ,\mathrm{sin}\;\alpha \right):-\pi /4{< }\alpha {< }\pi /4\right\},\end{equation*}$$

which is shown in figure 1.

An analysis of the singularities that are visible by the Compton data was done in section 3.2. For the point O, the angle defined by (3.21) gives β_m = 1.23 and the cone of visible directions given by (3.22) is the vertical cone with angles from the vertical between −β_m and β_m since the parameter r_M = 9. A calculation shows that 2(π/4 + β_m) > π and this implies C_R ∪ C_T = S¹ and we have a full resolution of the image singularities at O.

However, for points near the bottom x₂ = −3, there are invisible singularities that are not visible from either the Compton or x-ray data. For example, the vertical direction (0, 1) is not normal to any circle in the Compton data set at any point (0, x₂) for x₂ ∈ [−2.5, −3]. Figure 7 shows the points for which all wavefront directions are visible at those points (yellow color-roughly for points (x₁, x₂) for x₂ > −2) and near the bottom of the reconstruction region, there are more missing directions.

Zoom In Zoom Out Reset image size

In this section we detail our joint reconstruction scheme and lambda tomography regularization technique, and show the effectiveness of our methods in combatting the artefacts observed and predicted by our microlocal theory. We first explain the physical relationship between μ_E and n_e, which will be needed later in the formulation of our regularized inverse problem.

5.1. Relating μ_E and n_e

The attenuation coefficient and electron density satisfy the formula [53, p 36]

$$\begin{equation}{\mu }_{E}\left(Z\right)={n}_{\mathrm{e}}\left(Z\right){\sigma }_{E}\left(Z\right),\end{equation} \tag{ 5.1 }$$

where σ_E denotes the electron cross section, at energy E. Here Z denotes the effective atomic number. In the proposed application in airport baggage screening (among many other applications such as medical CT) we are typically interested in the materials with low effective Z. Hence we consider the materials with Z < 20 in this paper. For large enough E and Z < 20, σ_E(Z) is approximately constant as a function of Z. Equivalently μ_E and n_E are approximately proportional for low Z and high E by equation (5.1). See figure 8. We see a strong correlation between n_e and μ_E when E = 100 keV and Z < 20, and even more so when E is increased to E = 1 MeV. The sample of materials considered consists of 153 compounds (e.g. wax, soap, salt, sugar, the elements) taken from the NIST database [29]. In this case σ_E ≈ ν for some $\nu \in \mathbb{R}$ is approximately constant and we have μ_E ≈ νn_e. For a given energy E, ν is the slope of the straight line fit as in figure 8. Throughout the rest of this paper, we set ν as the slope of the straight line in the left hand of figure 8 (i.e. ν ≈ 0.57), and present reconstructions of μ_E for E = 100 keV.

Zoom In Zoom Out Reset image size

5.2. Lambda regularization; the idea

In sections 3 and 4 we discovered that the R_Lμ_E and $\mathcal{T}{n}_{\mathrm{e}}$ data provide complementary information regarding the detection and resolution of edges in an image. More specifically the line integral data resolved singularities in an open cone C_R with central axis x₁ and the toric section integral data resolved singularities in a cone C_T with central axis x₂. So the overlapping cones C_R ∪ C_T give a greater coverage of S¹ than when considered separately. In figures 3, 4 and 6, this theory was later verified through reconstructions of a delta functions by (un)filtered backprojection.

For a further example, let us consider a more complicated phantom than a delta function, one which is akin to densities considered later for testing our joint reconstruction and lambda regularization method. In figure 9 we have presented reconstructions of an image phantom f (with no noise) from R_Lf (transmission data—middle figure) using FBP, and from $\mathcal{T}f$ (Compton data—right figure) by an application of ${\mathcal{T}}^{{\ast}}\frac{{\mathrm{d}}^{2}}{\mathrm{d}{r}^{2}}$ (a contour reconstruction). In the reconstruction from Compton data, we see that the image singularities are well resolved in the vertical direction (x₂), and conversely in the horizontal direction (x₁) in the reconstruction from transmission data. In the middle picture (reconstruction from R_L), the visible singularities of the object are tangent to lines in the data set (normal wavefront set) and the artefacts are along lines at the end of the data set that are tangent to boundaries of the objects. In the right-hand reconstruction from Compton data, the visible boundaries are tangent to circles in the data set and the streaks are along circles at the end of the data set. Note that the visible boundaries in each picture complement each other and together, image the full objects. This is all as predicted by the theory of sections 3 and 4 (and is consistent with the theory in [5, 15]) and highlights the complementary nature of the Compton and transmission data in their ability to detect and resolve singularities.

Zoom In Zoom Out Reset image size

Given the complementary edge resolution capabilities of R_Lμ_E and $\mathcal{T}{n}_{\mathrm{e}}$, and given the approximate linear relationship between μ_E and n_e, we can devise a joint linear least squares reconstruction scheme with the aim to recover the image singularities stably in all directions in the n_e and μ_E images simultaneously. To this end we employ ideas in lambda tomography and microlocal analysis (figure 10).

Let $f\in \mathcal{E}\left({\mathbb{R}}^{n}\right)$ and let $Rf\left(s,\theta \right)=R{f}_{\theta }\left(s\right){=\int }_{{L}_{s,\theta }}f\mathrm{d}l$ denote the hyperplane Radon transform of f, where L_s,θ is as defined in section 4. The Radon projections Rf_θ detect singularities in f in the direction Θ = (cos θ, sin θ) (i.e. the elements (x, Θ) ∈ WF(f)). Applying a derivative filter $\frac{{\mathrm{d}}^{m}}{\mathrm{d}{s}^{m}}{R}_{\theta }$, for some m ⩾ 1, increases the strength of the singularities in the Θ direction by order m in Sobolev scale. Given $f,g\in \mathcal{E}\left({\mathbb{R}}^{n}\right)$, we aim to enforce a similarity in WF(f) and WF(g) through the addition of the penalty term ${\Vert}\frac{{\mathrm{d}}^{m}}{\mathrm{d}{s}^{m}}R\left(f-g\right){{\Vert}}_{{L}^{2}\left(\mathbb{R}{\times}{S}^{1}\right)}$ to the least squares solution. Note that we are integrating over all directions in S¹ to enforce a full directional similarity in WF(f) and WF(g). Specifically in our case f = μ_E, g = n_e and we aim to minimize the quadratic functional

$$\begin{equation}\underset{{\mu }_{E},{n}_{\mathrm{e}}}{\text{arg min}}{{\Vert}\left(\begin{pmatrix}\hfill w{R}_{L}\hfill & \hfill 0\hfill \\ 0\hfill & \hfill \mathcal{T}\hfill \\ \hfill \alpha \left[\frac{{\mathrm{d}}^{m}}{\mathrm{d}{s}^{m}}R\right.\hfill & \hfill \left.-\nu \frac{{\mathrm{d}}^{m}}{\mathrm{d}{s}^{m}}R\right]\hfill \end{pmatrix}\right)\left(\begin{pmatrix}\hfill {\mu }_{E}\hfill \\ \hfill {n}_{\mathrm{e}}\hfill \end{pmatrix}\right)-\left(\begin{pmatrix}\hfill w{b}_{1}\hfill \\ \hfill {b}_{2}\hfill \\ 0\hfill \end{pmatrix}\right){\Vert}}_{2}^{2},\end{equation} \tag{ 5.2 }$$

where R_L denotes a discrete, limited data Radon operator, R is the discrete form of the full data Radon operator, $\mathcal{T}$ is the discrete form of the toric section transform, b₁ is known transmission data and b₂ is the Compton scattered data. Here α is a regularization parameter which controls the level of similarity in the image wavefront sets. The lambda regularizers enforce the soft constraint that μ_E = νn_e (since $\frac{{\mathrm{d}}^{m}}{\mathrm{d}{s}^{m}}Rf=0\;\Longleftrightarrow\;f=0$ for $f\in {L}_{c}^{2}\left(X\right)$), but with emphasis on the location, direction and magnitude of the image singularities in the comparison. Further we expect the lambda regularizers to have a smoothing effect given the nature of $\frac{{\mathrm{d}}^{m}}{\mathrm{d}{s}^{m}}R$ as a differential operator (i.e. the inverse is a smoothing operation). Hence we expect α to also act as a smoothing parameter. The weighting $w={\Vert}\mathcal{T}{{\Vert}}_{2}/{\Vert}{R}_{L}{{\Vert}}_{2}$ is included so as to give equal weighting to the transmission and scattering datasets in the inversion. We denote the joint reconstruction method using lambda tomography regularizers as 'JLAM'. A common choice for m in lambda tomography applications is m = 2 [7, 43] (hence the name 'lambda regularizers'). With complete x-ray data, the application of a Lambda term yields this ${R}^{{\ast}}\frac{{d}^{2}}{\mathrm{d}{s}^{2}}Rf=-4\pi \sqrt{-{\Delta}}f$ [31, example 9], so the singularities of f are preserved and emphasized by order 1 in Sobolev scale, so they will dominate the lambda reconstruction. Hence choosing m = 2 is sufficient for a full recovery of the image singularities. Since the singularities are dominant in the lambda term, they are matched accurately in (5.2). Indeed we have already seen the effectiveness of such a filtering approach in recovering the image contours earlier in the right-hand of figure 9. We find that setting m = 2 here works well as a regularizer on synthetic image phantoms and simulated data with added pseudo random noise, as we shall now demonstrate. We note that the derivative filters for m ≠ 2 are also worth exploration but we leave such analysis for future work.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

5.3. Proposed testing and comparison to the literature

To test our reconstruction method, we first consider two test phantoms, one simple and one complex (as in [55]). See figure 11. The phantoms considered are supported on Γ, the region in figure 7 in which there is full wavefront coverage from joint x-ray and Compton scattered data. The simple density phantom consists of a polyvinyl chloride (PVC) cuboid and an aluminium sphere with an approximate density ratio of 1:2 (PVC:Al). The complex density phantom consists of a water ellipsoid, a sulphur ellipsoid, a calcium sulphate (CaSO₄) right-angled-triangle and a thin film of titanium dioxide (TiO₂) in the shape of a cross. The density ratio of the materials which compose the complex phantom is approximately 1:2:3:4 (H₂O:S:CaSO₄:TiO₂). The density values used are those of figure 8 taken from the NIST database [28], and the background densities (≈0) were set to the density of dry air (near sea level). The corresponding attenuation coefficient phantoms are simulated similarly. The materials considered are widely used in practice. For example CaSO₄ is used in the production of plaster of Paris and stucco (a common construction material) [57], and TiO₂ is used in the making of decorative thin films (e.g. topaz) and in pigmentation [56, p 15].

To simulate data we set

$$\begin{equation}b=\left(\begin{pmatrix}\hfill {b}_{1}\hfill \\ \hfill {b}_{2}\hfill \end{pmatrix}\right)=\left(\begin{pmatrix}\hfill {R}_{L}{\mu }_{E}\hfill \\ \hfill \mathcal{T}{n}_{\mathrm{e}},\hfill \end{pmatrix}\right)\end{equation} \tag{ 5.3 }$$

and add a Gaussian noise

$$\begin{equation}{b}_{\eta }=b+\eta {\Vert}b{{\Vert}}_{2}\frac{{v}_{G}}{\sqrt{l}},\end{equation} \tag{ 5.4 }$$

for some noise level η, where l is the length of b and v_G is a vector of length l of draws from $\mathcal{N}\left(0,1\right)$. For comparison we present separate reconstructions of μ_E and n_e using total variation (TV regularizers). That is we will find

$$\begin{equation}\underset{{\mu }_{E}}{\text{arg min}}{\Vert}{R}_{L}{\mu }_{E}-{b}_{1}{\parallel }_{2}^{2}+\alpha \text{TV}\left({\mu }_{E}\right)\end{equation} \tag{ 5.5 }$$

to reconstruct μ_E and

$$\begin{equation}\underset{{n}_{\mathrm{e}}}{\text{arg min}}{\Vert}\mathcal{T}{n}_{\mathrm{e}}-{b}_{2}{\parallel }_{2}^{2}+\alpha \text{TV}\left({n}_{\mathrm{e}}\right)\end{equation} \tag{ 5.6 }$$

for n_e, where TV(f) = ||∇f||₁ and α > 0 is a regularization parameter. We will denote this method as 'TV'. In addition we present reconstructions using the state-of-the-art joint reconstruction and regularization techniques from the literature, namely the joint total variation (JTV) methods of [20] and the linear parallel level sets (LPLS) methods of [12]. To implement JTV we minimize

$$\begin{equation}\underset{{\mu }_{E},{n}_{\mathrm{e}}}{\text{arg min}}{{\Vert}\left(\begin{pmatrix}\hfill w{R}_{L}\hfill & \hfill 0\hfill \\ 0\hfill & \hfill \mathcal{T}\hfill \end{pmatrix}\right)\left(\begin{pmatrix}\hfill {\mu }_{E}\hfill \\ \hfill {n}_{\mathrm{e}}\hfill \end{pmatrix}\right)-\left(\begin{pmatrix}\hfill w{b}_{1}\hfill \\ \hfill {b}_{2}\hfill \end{pmatrix}\right){\Vert}}_{2}^{2}+\alpha {\text{JTV}}_{\beta }\left({\mu }_{E},{n}_{\mathrm{e}}\right),\end{equation} \tag{ 5.7 }$$

where $w={\Vert}\mathcal{T}{{\Vert}}_{2}/{\Vert}{R}_{L}{{\Vert}}_{2}$ as before, and

$$\begin{equation}{\text{JTV}}_{\beta }\left({\mu }_{E},{n}_{\mathrm{e}}\right)={\int }_{\left[-2,2\right]{\times}\left[-3,1\right]}{\left({\Vert}\nabla {\mu }_{E}\left(\mathbf{\text{x}}\right){{\Vert}}_{2}^{2}+{\Vert}\nabla {n}_{\mathrm{e}}\left(\mathbf{\text{x}}\right){{\Vert}}_{2}^{2}+{\beta }^{2}\right)}^{\frac{1}{2}}\mathrm{d}\mathbf{\text{x}},\end{equation} \tag{ 5.8 }$$

where β > 0 is an additional hyperparameter included so that the gradient of JTV_β is defined at zero. This allows one to apply techniques in smooth optimization to solve (5.7).

To implement LPLS we minimize

$$\begin{equation}\underset{{\mu }_{E},{n}_{\mathrm{e}}}{\text{arg min}}{{\Vert}\left(\begin{pmatrix}\hfill w{R}_{L}\hfill & \hfill 0\hfill \\ 0\hfill & \hfill \mathcal{T}\hfill \end{pmatrix}\right)\left(\begin{pmatrix}\hfill {\mu }_{E}\hfill \\ \hfill {n}_{\mathrm{e}}\hfill \end{pmatrix}\right)-\left(\begin{pmatrix}\hfill w{b}_{1}\hfill \\ \hfill {b}_{2}\hfill \end{pmatrix}\right){\Vert}}_{2}^{2}+\alpha {\text{LPLS}}_{\beta }\left({\mu }_{E},{n}_{\mathrm{e}}\right),\end{equation} \tag{ 5.9 }$$

where

$$\begin{equation}{\text{LPLS}}_{\beta }\left({\mu }_{E},{n}_{\mathrm{e}}\right)={\int }_{\left[-2,2\right]{\times}\left[-3,1\right]}{\Vert}\nabla {\mu }_{E}\left(\mathbf{\text{x}}\right){{\Vert}}_{\beta }{\Vert}\nabla {n}_{\mathrm{e}}\left(\mathbf{\text{x}}\right){{\Vert}}_{\beta }-\vert \nabla {\mu }_{E}\left(\mathbf{\text{x}}\right)\cdot \nabla {n}_{\mathrm{e}}\left(\mathbf{\text{x}}\right){\vert }_{{\beta }^{2}}\;\mathrm{d}\mathbf{\text{x}},\end{equation} \tag{ 5.10 }$$

where $\vert \mathbf{\text{x}}{\vert }_{\beta }=\sqrt{\vert \mathbf{\text{x}}{\vert }^{2}+{\beta }^{2}}$ and ${\Vert}\mathbf{\text{x}}{{\Vert}}_{\beta }=\sqrt{{\Vert}\mathbf{\text{x}}{{\Vert}}_{2}^{2}+{\beta }^{2}}$ for β > 0. The JTV and LPLS penalties seek to impose soft constraints on the equality of the image wavefront sets of μ_E and n_e. For example setting β = 0 in the calculation of LPLS_β yields

$$\begin{align}\hfill {\text{LPLS}}_{\beta }\left({\mu }_{E},{n}_{\mathrm{e}}\right)& ={\Vert}\nabla {\mu }_{E}{{\Vert}}_{2}{\Vert}\nabla {n}_{\mathrm{e}}{{\Vert}}_{2}-\vert \nabla {\mu }_{E}\cdot \nabla {n}_{\mathrm{e}}\vert ={\Vert}\nabla {\mu }_{E}{{\Vert}}_{2}{\Vert}\nabla {n}_{\mathrm{e}}{{\Vert}}_{2}\left(1-\vert \mathrm{cos} \theta \vert \right),\hfill \end{align} \tag{ 5.11 }$$

where θ is the angle between ∇n_e and ∇μ_E. Hence (5.11) is minimized for the gradients which are parallel (i.e. when θ = 0, π), and thus using LPLS_β as regularization serves to enforce equality in the image gradient direction and location (i.e. when LPLS_β is small, the gradient directions are approximately equal).

We wish to stress that the comparison with JTV and LPLS is included purely to illustrate the potential advantages (and disadvantages) of the lambda regularizers when compared to the state-of-the-art regularization techniques. Namely is the improvement in image quality due to joint data, lambda regularizers or are they both beneficial? We are not claiming a state-of-the-art performance using JLAM, but our results show JLAM has good performance, and it is numerically easier to implement, requiring only least squares solvers. There are two hyperparameters (α and β) to be tuned in order to implement JTV and LPLS, which is more numerically intensive (e.g. using cross validation) in contrast to JLAM with only one hyperparameter (α). Moreover the LPLS objective is non-convex [12, appendix A], and hence there are potential local minima to contend with, which is not an issue with JLAM, being a simple quadratic objective.

To minimize (5.2), we store the discrete forms of R_L, R and $\mathcal{T}$ as sparse matrices and apply the conjugate gradient least squares (CGLS) solvers of [22, 23] (specifically the 'IRnnfcgls' code) with non-negativity constraints (since the physical quantities n_e and μ_E are known a priori to be nonnegative). To solve equations (5.5) and (5.6) we apply the heuristic least squares solvers of [22, 23] (specifically the 'IRhtv' code) with TV penalties and non-negativity constraints. To solve (5.7) and (5.9) we apply the codes of [12], modified so as to suit a Gaussian noise model (a Poisson model is used in [12, equation 3]). The relative reconstruction error is calculated as = ||x − y||₂/||x||₂, where x is the ground truth image and y is the reconstruction. For all methods compared against we simulate data and added noise as in equations (5.3) and (5.4), and the noise level added for each simulation is η = 0.1 (10% noise). We choose α for each method such that is minimized for a noise level of η = 0.1. That is we are comparing the best possible performance of each method. We set β for JTV and LPLS to the values used on the 'lines2' data set of [12]. We do not tune β to the best performance (as with α) so as to give fair comparison between TV, JLAM, JTV and LPLS. After the optimal hyperparameters were selected, we performed 100 runs of TV, JLAM, JTV and LPLS on both phantoms for 100 randomly selected sets of materials. That is, for 100 randomly chosen sets of values from figure 8 and the NIST database, with the NIST values corresponding to the nonzero parts of the phantoms. We present the mean (μ) and standard deviation (σ) relative errors over 100 runs in the left-hand of tables 2 and 4 for the simple and complex phantom respectively. The results are given in the form _± = μ ± σ for each method. In addition to the relative error , we also provide metrics to measure the structural accuracy of the results. Specifically we will compare F-scores on the image gradient and support, as is done in [2, 49]. The gradient F-score [2] measures the wavefront set reconstruction accuracy, and the support F-score [49, p 5] (see DICE score) is a measure of the geometric accuracy. That is, the support F-score checks whether the reconstructed phantoms are the correct shape and size. The F-score takes values on [0, 1]. For this metric, values close to one indicate higher performance, and conversely for values close to zero. Similarly to , we present the F-scores of the randomized tests in the form F_± = μ_F ± σ_F, where μ_F and σ_F are the mean and standard deviation F-scores respectively. In all tables, the support F-scores are labelled by supp(n_e) and supp(μ_E), and by ∇n_e and ∇μ_E for the gradient F-scores.

5.4. Results and discussion

See figure 11 for image reconstructions of the simple phantom using TV, JLAM, JTV and LPLS, and see table 1 for the corresponding and F-score values. See table 2 for the _± and F_± values calculated over 100 randomized simple phantom reconstructions. For the complex phantom, see figure 12 for image reconstructions, and table 3 for the and F-score values. See table 4 for _± and F_±. In the separate reconstruction of n_e (using method TV) we see a blurring of the ground truth image edges (wavefront directions) in the horizontal direction and there are artefacts in the reconstruction due to limited data, as predicted by our microlocal theory. In the TV reconstruction of μ_E we see a similar effect, but in this case we fail to resolve the wavefront directions in the vertical direction due to limited line integral data. This is as predicted by the theory of section 4 and [5]. In section 3 we discovered the existence also of nonlocal artefacts in the n_e reconstruction, which were induced by the mappings λ_ij. However these were found to lie largely outside the imaging space unless the singularity in question (x, ξ) ∈ WF(n_e) were such that x is close to the detector array (see figures 3 and 4). Hence why we do not see the effects of the λ_ij in the phantom reconstructions, as the phantoms are bounded sufficiently away from the detector array. The added regularization may smooth out such artefacts also, which was found to be the case in [55].

Table 1. Simple phantom and F-score comparison using TV, JLAM, JTV and LPLS.

	TV	JLAM	JTV	LPLS
ne	0.26	0.14	0.05	0.03
μE	0.40	0.15	0.05	0.03
F-score	TV	JLAM	JTV	LPLS
supp(ne)	0.81	0.99	0.99	∼1
∇ne	0.76	0.87	0.82	0.87
supp(μE)	0.78	∼1	∼1	∼1
∇μE	0.49	0.87	0.82	0.86

Table 2. Randomized simple phantom _± and F_± comparison over 100 runs using JLAM, JTV and LPLS.

±	JLAM	JTV	LPLS
ne	0.15 ± 0.01	0.05 ± 0.005	0.06 ± 0.002
μE	0.16 ± 0.02	0.05 ± 0.005	0.06 ± 0.03
F±	JLAM	JTV	LPLS
supp(ne)	0.99 ± 0.01	0.99 ± 0.004	∼1 ± 0.005
∇ne	0.86 ± 0.01	0.84 ± 0.02	0.86 ± 0.04
supp(μE)	0.98 ± 0.04	0.99 ± 0.005	∼1 ± 0.005
∇μE	0.86 ± 0.01	0.85 ± 0.02	0.87 ± 0.04

Zoom In Zoom Out Reset image size

Table 3. Complex phantom and F-score comparison using TV, JLAM, JTV and LPLS.

	TV	JLAM	JTV	LPLS
ne	0.36	0.24	0.16	0.09
μE	0.63	0.30	0.19	0.13
F-score	TV	JLAM	JTV	LPLS
supp(ne)	0.78	0.98	0.97	0.99
∇ne	0.73	0.83	0.84	0.84
supp(μE)	0.65	0.94	0.98	0.99
∇μE	0.39	0.83	0.84	0.85

Table 4. Randomized complex phantom _± and F_± comparison over 100 runs using JLAM, JTV and LPLS.

±	JLAM	JTV	LPLS
ne	0.25 ± 0.02	0.13 ± 0.03	0.13 ± 0.03
μE	0.31 ± 0.03	0.16 ± 0.02	0.17 ± 0.04
F±	JLAM	JTV	LPLS
supp(ne)	0.98 ± 0.01	0.98 ± 0.005	0.99 ± 0.004
∇ne	0.76 ± 0.06	0.78 ± 0.05	0.75 ± 0.05
supp(μE)	0.90 ± 0.06	0.96 ± 0.03	0.98 ± 0.007
∇μE	0.73 ± 0.07	0.77 ± 0.05	0.74 ± 0.05

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 5. Al bar phantom and F-score comparison using TV, JLAM, JTV and LPLS.

	TV	JLAM	JTV	LPLS
ne	0.28	0.09	0.04	0.02
μE	0.68	0.11	0.09	0.07
F-score	TV	JLAM	JTV	LPLS
supp(ne)	0.71	0.99	∼1	∼1
∇ne	0.64	0.82	0.79	0.83
supp(μE)	0.54	∼1	∼1	∼1
∇μE	0.53	0.82	0.80	0.90

Table 6. Randomized bar phantom _± and F_± comparison over all NIST materials considered (153 runs) using JLAM, JTV and LPLS.

±	JLAM	JTV	LPLS
ne	0.10 ± 0.02	0.04 ± 0.004	0.03 ± 0.02
μE	0.12 ± 0.03	0.12 ± 0.03	0.08 ± 0.03
F±	JLAM	JTV	LPLS
supp(ne)	0.99 ± 0.02	∼1 ± 0.003	∼1 ± 0.001
∇ne	0.83 ± 0.02	0.79 ± 0.02	0.84 ± 0.02
supp(μE)	0.98 ± 0.06	0.99 ± 0.02	∼1 ± 0.008
∇μE	0.81 ± 0.03	0.77 ± 0.02	0.85 ± 0.02

Using the joint reconstruction methods (i.e. JLAM, JTV and LPLS) we see a large reduction in the image artefacts in n_e and μ_E, since with joint data we are able to stably resolve the image singularities in all directions. The improvement in and the F-score is also significant, particularly in the μ_E reconstruction. Thus it seems that the joint data is the greater contributor (over the regularization) to the improvement in the image quality, as the approaches with joint data each perform well. Upon comparison of JLAM, JTV and LPLS, the metrics are significantly improved when using JTV and LPLS over JLAM, but the image quality and F-scores are highly comparable. This indicates that, while the noise in the reconstruction is higher using JLAM, the recovery of the image edges and support is similar using JLAM, JTV and JLAM. As theorized, the lambda regularizers were successful in preserving the wavefront sets of μ_E and n_e. However there is a distortion present in the nonzero parts of the JLAM reconstruction. This is the most notable difference in JLAM and JTV/LPMS, and is likely the cause of the discrepancy. So while the edge preservation and geometric accuracy of JLAM is of a high quality (and this was our goal), the smoothing properties of JLAM are not up to par with the state-of-the-art currently. We note however that the JTV and LPLS objectives are nonlinear (with LPLS also non-convex) and require significant additional machinery (e.g. in the implementation of the code of [12] used here) in the inversion when compared to JLAM, which is a straight forward implementation of linear least squares solvers.

5.5. Reconstructions with limited data

Here we have introduced a new joint reconstruction method 'JLAM' for low effective Z imaging (Z < 20), based on ideas in lambda tomography. We considered primarily the 'parallel line segment' geometry of [54], which is motivated by system architectures for airport security screening applications. In section 3 we gave a microlocal analysis of the toric section transform $\mathcal{T}$, which was first proposed in [54] for a CST problem. Explicit expressions were provided for the microlocal artefacts and verified through simulation. Section 4 explained the x-ray CT artefacts using the theory of [5]. Following the theory of sections 3 and 4, we detailed the JLAM algorithm in section 5. Here we conducted simulation testing and compared JLAM to separate reconstructions using TV, and to the nonlinear joint inversion methods, JTV [20] and LPLS [12] from the literature. The joint inversion methods considered (i.e. JLAM, JTV and LPLS) were successful in preserving the image contours in the reconstruction, as predicted. However the smoothing applied by JLAM was not as effective as JTV and LPLS, and we saw a distortion in the JLAM reconstruction (see figures 11 and 12). JTV and LPLS were thus shown to offer better performance than JLAM, with LPLS producing the best results overall. The advantages of JLAM over JTV and LPLS are in the fast, linear inversion, and the reduction in tuning parameters (one for JLAM, two for JTV/LPLS). Given the linearity of JLAM, the ideas of JTV and LPLS can be easily integrated with lambda regularization to modify the objectives of the literature and improve further the edge resolution of the reconstruction. To preserve the linearity of JLAM we could also combine JLAM with a Tikhonov regularizer. This may help smooth out the distortion observed in the JLAM reconstruction. We leave such ideas for future work.

We would like to thank the journal reviewers for their helpful comments and insight towards this article, particularly in regards to the simulation study. This material is based upon work supported by the U.S. Department of Homeland Security, Science and Technology Directorate, Office of University Programs, under Grant Award 2013-ST-061-ED0001. The views and conclusions contained in this document are those of the authors and should not be interpreted as necessarily representing the official policies, either expressed or implied, of the U.S. Department of Homeland Security. The work of the second author was partially supported by US National Science Foundation Grant DMS 1712207. The authors thank the referees for thorough reviews and thoughtful comments that improved the article.

The generation of the plots of figure 8 is explained in more detail here. We will explain the generation of the plot for E = 100 keV. Refer to figure A1. We first plotted μ_E for E = 100 keV against n_e for all materials in the NIST database [28] with effective Z less than 20. This is the left-hand plot of figure A1. The set of materials with effective Z < 20 was

$$\begin{equation*}{Z}_{\text{eff}}=\left\{Z:{\sigma }_{E}\left(Z\right){< }{\sigma }_{E}\left(20\right),E=100\;\text{keV}\right\},\end{equation*}$$

where σ_E is the electron cross section. We noticed a large outlier (coal, or amorphous carbon) which corrupts the correlation in our favour, and hence we chose to remove the material from consideration in simulation. The outlier is highlighted in the left hand plot. After the outlier was removed we noticed a number of materials located at the origin (with negligible attenuation coefficient and density, such as air) in the middle scatter plot of figure A1. As such materials again bias the correlation and plot standard deviation in our favour, these were removed to produce the left hand plot of figure 8 in the right-hand of figure A1. The same points were removed in the generation of the right-hand plot of figure 8 also, for E = 1 MeV.

Zoom In Zoom Out Reset image size