A joint reconstruction and lambda tomography regularization technique for energy-resolved X-ray imaging
\ddmmyyyydate \currenttime

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 $\mu_{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, 51, 41, 43, 55], and the sources are fan-beam (in the plane) with opening angle $\pi$ (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, $\log\left(\frac{I_{0}}{I_{A}}\right)=\int_{L}\mu_{E}\mathrm{d}l$ [31, page 2] where $I_{A}$ is the photon intensity measured at the detector, $I_{0}$ is the initial source intensity and $\mu_{E}$ is the attenuation coefficient at energy $E$. Here $L$ is a line through $S$ and $D_{A}$, with arc measure $\mathrm{d}l$. Hence the transmission data determines a set of integrals of $\mu_{E}$ over lines, and the problem of reconstructing $\mu_{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 $\mu_{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]

(1.1)

$$I_{C}=\int_{T}n_{e}\mathrm{d}t,$$

where $I_{C}$ is the Compton scattered intensity measured at a point on $D_{C}$. A toric section $T=C_{1}\cup C_{2}$ is the union of two intersecting circles of the same radii (as displayed in figure 2), and $\mathrm{d}t$ 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.

The line and circular arc Radon transforms with full data, are known [35, 30] 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 $\|\frac{\mathrm{d}^{m}}{\mathrm{d}s^{m}}R(\mu_{E}-n_{e})\|_{L^{2}(\mathbb{R}\times S^{1})}$ for some $m\geq 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 $\mu_{E}$ and $n_{e}$. Further $\frac{\mathrm{d}^{m}}{\mathrm{d}s^{m}}R$ for $m\geq 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, page 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, 19, 40, 46, 47, 6, 50, 12, 20, 11, 48, 10, 45, 4]. 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 $\mu_{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 $\mu_{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}^{*}\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 artifacts in the X-ray CT reconstruction. In section 5, we detail our joint reconstruction method for the simultaneous reconstruction of $\mu_{E}$ and $n_{e}$. Later in section 5.4 we present simulated reconstructions of $\mu_{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}(X)$ be the space of smooth functions compactly supported on $X$ with the standard topology and let $\mathcal{D}^{\prime}(X)$ denote its dual space, the vector space of distributions on $X$. Let $\mathcal{E}(X)$ be the space of all smooth functions on $X$ with the standard topology and let $\mathcal{E}^{\prime}(X)$ denote its dual space, the vector space of distributions with compact support contained in $X$. Finally, let $\mathcal{S}({{\mathbb{R}}^{n}})$ be the space of Schwartz functions, that are rapidly decreasing at $\infty$ along with all derivatives. See [44] for more information.

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

$$\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.$$

If $f$ is a function of $({\mathbf{y}},{\mathbf{x}},{\boldsymbol{\sigma}})$ then $\partial^{\alpha}_{\mathbf{y}}f$ and $\partial^{\alpha}_{\boldsymbol{\sigma}}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 $\phi$ is a function of $({\mathbf{y}},{\mathbf{x}},{\boldsymbol{\sigma}})\in Y\times X\times{{\mathbb{R}}}^{N}$ then we define $\mathrm{d}_{{\mathbf{y}}}\phi=\left(\frac{\partial\phi}{\partial y_{1}},\frac{\partial\phi}{\partial y_{2}},\cdots,\frac{\partial\phi}{\partial y_{n}}\right)$, and $\mathrm{d}_{\mathbf{x}}\phi$ and $\mathrm{d}_{\boldsymbol{\sigma}}\phi$ are defined similarly. We let $\mathrm{d}\phi=\left(\mathrm{d}_{{\mathbf{y}}}\phi,\mathrm{d}_{{\mathbf{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, page 16]:

This definition follows the intuitive idea that the elements of $\mathrm{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 $\mathrm{WF}(f)=\{({\mathbf{x}},t{\mathbf{x}}):{\mathbf{x}}\in S^{2},t\neq 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^{*}(X)=X\times{{\mathbb{R}}^{n}}$ and consider $\mathrm{WF}(f)$ as a subset of $X\times{{\mathbb{R}}^{n}}\setminus\left\{\mathbf{0}\right\}$.

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

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, 27, 26].

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

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 $\delta$ functions, before proving our main microlocal theory.

Let $r_{m}>1$ and define the set of points to be scanned as

$$X:=\{(x_{1},x_{2})\in{{\mathbb{R}}^{2}}\hskip 0.85358pt:\hskip 0.85358pt2-r_{m}<x_{2}<1\}.$$

Note that $r_{m}$ controls the depth of the scanning tunnel as in figures 1 and 2. Let

$$Y:=(0,\infty)\times\mathbb{R}$$

then for $j=1,2$, and $(s,x_{0})\in Y$, we define the circles $C_{j}$ and their centers $\mathbf{c}_{j}$

(3.1)

$$\begin{gathered}r=\sqrt{s^{2}+1},\quad\mathbf{c}_{j}(s,x_{0})=((-1)^{j}s+x_{0},2)\\ C_{j}(s,x_{0})=\{{\mathbf{x}}\in\mathbb{R}^{2}:|{\mathbf{x}}-\mathbf{c}_{j}(s,x_{0})|^{2}-s^{2}-1=0\}.\end{gathered}$$

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

(3.2)

$$\mathcal{T}_{1}f(s,x_{0})=\int_{C_{1}}f\mathrm{d}s,\ \ \ \ \ \ \ \ \mathcal{T}_{2}f(s,x_{0})=\int_{C_{2}}f\mathrm{d}s.$$

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

(3.3)

$$\mathcal{T}f(s,x_{0})=\int_{C_{1}\cup C_{2}}f\mathrm{d}s=\mathcal{T}_{1}(f)(s,x_{0})+\mathcal{T}_{2}(f)(s,x_{0})$$

where $\mathrm{d}s$ denotes the arc element on a circle and $(s,x_{0})\in Y$.

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

(3.4)

$$\begin{split}\mathcal{T}f(s,x_{0})&=\mathcal{T}_{1}f(s,x_{0})+\mathcal{T}_{2}f(s,x_{0})\\ &=\frac{1}{2r}\sum_{j=1}^{2}\int_{\mathbb{R}^{2}}\delta(|{\mathbf{x}}-\mathbf{c}_{j}(s,x_{0})|^{2}-s^{2}-1)f({\mathbf{x}})\mathrm{d}{\mathbf{x}}\\ &=\frac{1}{2r}\sum_{j=1}^{2}\int_{-\infty}^{\infty}\int_{\mathbb{R}^{2}}e^{-i\sigma(|{\mathbf{x}}-((-1)^{j}s+x_{0},2)|^{2}-s^{2}-1)}f({\mathbf{x}})\mathrm{d}{\mathbf{x}}\mathrm{d}\sigma.\end{split}$$

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(|{\mathbf{x}}-\mathbf{c}_{j}(s,x_{0})|-\sqrt{s^{2}+1}\right)f(x)\,\mathrm{d}{\mathbf{x}}$. So the toric section transform is the sum of two FIO’s with phase functions

$$\phi_{j}(s,x_{0},{\mathbf{x}},\sigma)=\sigma(|{\mathbf{x}}-((-1)^{j}s+x_{0},2)|^{2}-s^{2}-1)$$

for $j=1,2$. Our distributions $f$ are supported away from the intersection points of $C_{1}$ and $C_{2}$, 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}$.

Because $\mathcal{C}_{j}$ satisfies the Bolker Assumption, the composition $\mathcal{C}^{*}_{j}\circ\mathcal{C}_{j}\subset\Delta$, where $\Delta$ 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_{1}(s,x_{0})\cap C_{2}(s,x_{0}))\cap\text{supp}(f)=\emptyset$ for any $(s,x_{0})\in Y$.

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

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}(X)$ and $\mathcal{T}_{j}^{t}$ is not defined for such distributions. Therefore, we introduce a smooth cutoff function. Choose $r_{M}>2$ and let $\psi(s)$ be a smooth compactly supported function equal to one for $s\in\left[1,\sqrt{1-r^{2}_{M}}\right]$ and define

(3.14)

$$\mathcal{T}_{j}^{*}g=\mathcal{T}^{t}_{j}(\psi g)$$

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