1. Introduction

In this paper we present a novel reconstruction technique for two-dimensional BST, based on the generalized Radon transform models of [1]. Our method employs a combination of ideas from multibang control, Compressed Sensing (CS), and microlocal analysis to derive new regularization penalties, which prove to be effective in combatting the high level of noise and systematic error in BST data (e.g. error due to beam attenuation and Compton scatter, the two primary physical mechanisms not accounted for under the BST model). Such ideas have been applied in conventional X-ray CT, for example, CS has been applied to reduce the radiation dose and the number of projections in medical imaging [2–5], and microlocal analysis is used to reduce image artifacts with limited X-ray CT data [6]. We apply similar ideas here to BST. The BST model arises from a scanning geometry (first introduced in [1]) which uses translating sources to inspect what is well approximated as a line in image space. See Fig. 1. A more detailed description of the sensing geometry is given later in section 2.1. Here BST refers to the imaging of the Bragg differential cross section function (denoted by $f$ in this paper) from Bragg scatter data, and is not exclusive to the sensing geometry of [1].

 

Fig. 1. The X-ray scanner geometry. The scanning target (i.e., the spherical object pictured) is contained in the scanning tunnel shown in Fig. 1(a). The detectors are collimated to planes, and the scattering events occur along lines $L=\{x_2=x',x_3=0\}$, for some $-w_{x_2}<x'<w_{x_2}$. The scatter from $L$ is measured exclusively by detectors with offset $\epsilon$, as pictured.

Download Full Size | PDF

The literature considers a variety of reconstruction techniques and experimental methods in BST [7–12]. In [8] the authors consider snap shot, pencil beam, coded aperture imaging of crystalline powders, such as NaCl and Al. The technique, referred to as Coded Aperture X-ray Scatter Imaging (CAXSI), uses a pencil beam source to illuminate a small sample of crystalline powder. The scattered rays are then passed through a coded aperture mesh, and the resulting intensity is recorded by a linear (1-D) array of energy-resolved detectors. The use of coded aperture offers information about the scattering directions of the incoming photons, and thus improves the problem stability. Mathematically, the coded aperture modelling is represented by a kernel weighting ($t$ in [8, Eq. (3)]). The physical modeling then leads the authors to a linear inverse problem to solve for the Bragg differential cross section function $f$. To obtain a solution the authors apply Total Variation (TV) regularization and minimize the Poisson log-likelihood function. An iterative, Generalized Expectation-Maximization (GEM) algorithm is then implemented to minimize the objective, with good results on experimental data.

In [9] CAXSI is considered, using an experimental setup much like that of [8] with a planar (2-D) array of detectors (the detectors are not energy-resolved). A Generalized Maximum Likelihood (GML) estimator is then applied to estimate $f$. The GML algorithm is a multiplicative, iterative update intended to match the BST model to the data with a mean background and Poisson statistics applied. GML is similar to GEM, except GEM applies an additional maximization step [8, Eq. (10)] after each iteration, which amounts to Poisson denoising with TV.

In [7,10–12] a number of CAXSI variations are considered, for example, using fan beam sources in [7,10]. In [10], a set of “reference" differential cross section curves is used to detemine the object location, whereby each pixel is assigned to a material from the reference library based on the normalized correlation between the reconstructed and reference cross section values at the pixel. In this paper, we do not assume knowledge of a cross section reference library, and wish to keep the material properties general. We assume only that the form factor curves are $L^{2}$ functions.

There are a number of papers which consider Coherent Scatter Computed Tomography (CSCT) [13–15]. In those papers, low angle ($\omega <10^{\circ }$) coherent scatter is used to recover the coherent scatter form factor function. For example, in [14], a Filtered Back Projection (FBP) type formula, where the back-projection is performed along curved lines in 3-D space, is used to obtain a coherent form factor reconstruction. In the simulations conducted in this paper, we consider the scatter angles in the range $5.2^{\circ }<\omega <114.6^{\circ }$, and photon energies $E\in [1,29]$keV. We consider such an energy range so that our forward model satisfies the uniqueness guarantees of [1]. Practically speaking, this makes the inversion more stable, as there can be no additional artifacts in the image due to null space. For $E\in [1,29]$keV, the small angle scatter does not correspond to a meaningful range of momentum transfer values. Hence, we consider a larger range of scattering angles than what is typically considered in CSCT to compensate for the low energy range. In CSCT, a higher energy range than $E\in [1,29]$keV is typically used (e.g., a mean source energy of $43$keV in [13]), which allows for larger, more dense material imaging, as higher energy X-rays have higher penetration and are less likely to be absorbed. For $E\in [1,29]$keV, we see little signal when imaging large, dense objects, due to photoelectric absorption. Hence, the proposed methods are applicable to small sample, low $Z$ (i.e. $Z<20$) material at this stage. This is a limitation of our work, when compared to the CSCT literature. Small sample, low $Z$, imaging is considered in other CSCT work, however, e.g., in [8,9] (in those papers the scanned materials have $\approx 1$cm dimension), and hence the imaging problem considered in this paper is of interest to the CSCT community. In [13–15], the authors assume some knowledge of transmission data to correct for self-attenuation. We assume no knowledge of transmission data in our work, and obtain reconstructions purely from diffraction data. In [16,17], similar machine geometries to that of Fig. 1 are considered for threat detection. In those papers, low-angle scatter is also used (similarly to, e.g., [13]), and it is suggested to use transmission data for attenuation correction.

We introduce a new regularization scheme for small sample ($<3\textrm {cm}$ dimension), low effective atomic number $Z$ (i.e., $Z<20$) BST. A motivation for this work is the detection of small-sample crystallites which are concealed by clutter (e.g., this could be a small crystalline sample in a mail parcel at customs). The Monte Carlo simulations conducted later in section 3.6 consider this scenario, where the clutter is modeled by cellulose (a main component of common types of clutter, e.g., paper and clothing) and the samples are spherical crystallites with radius $r<3$cm. The technique we propose uses a collection of ideas from CS, multibang control and microlocal analysis. The application of such ideas has not yet been investigated in the BST literature. We also make no use of coded apertures in our experimental setup to restrict the scattering directions and increase problem stability. The linear collimation technology of the scanner of Fig. 1 is used to the same effect, to restrict the scattering sites to lines parallel to $x_1$ in the $(x_1,x_2)$ plane.

The reconstruction target $f$ (illustrated as a sphere in Fig. 1) is four-dimensional $f=f(q,x_1,x_2,x_3)$, where $q$ denotes the momentum transfer of the scattering interaction (see Eq. (2.7)). As in [1], we consider only 2-D reconstruction here (in $(q,x_1)$ space), with $x_2$ and $x_3=0$ fixed, in our simulations, and leave the piecing together of the 2-D slices (to form a full 4-D image) to future work. Note that, in the Monte Carlo data generation considered later in section 3.6, we include the attenuative and absorption effects from the full 3-D objects, so the densities at different $(x_1,x_2,x_3)$ are accounted for. Further, for a fixed $\epsilon$ (or detector row), the scattering contributions from multiple $x_2$ and $x_3$ are not considered in the model due to the linear collimation technology and fan-beam source (the collimation fixes $x_2$, and $x_3=0$ is fixed due to the fan-beam source), and the mapping $x_2\mapsto \epsilon$ is assumed to be diffeomorphic. In the next section, we will formalize the mapping between $x_2$ and $\epsilon$. A major part of the regularization idea we propose is to assume that $f$ is separable as a function of $q$ and $(x_1,x_2,x_3)$ (in the spatial domain). That is, we assume $f$ can be decomposed $f(q,x_1,x_2,x_3)=\sum _{j=1}^{l} f_j(q)\chi _j(x_1,x_2,x_3)$, where $l\in \mathbb {Z}_+$. Later, in section 2.3, we define the order of the model $l$, and the types of function $f_j$ and $\chi _j$ used. We model the spatial component of $f$, that is, the $\chi _j$ functions, using an overcomplete piecewise constant dictionary. This is a standard idea in CS [3,18,19], although typically in CS the library (basis) functions (or atoms [19]) are chosen to span the whole imaging domain. The basis functions we use only cover the $(x_1,x_2,x_3)$ domain, with $f$ having more general $L^{2}$ properties in the $q$ domain. The assumptions made here regarding the piecewise constant model for $f$ discussed above are consistent with the BST literature [7–12], and what is expected in practice. That is, we expect $f$ to be expressible as a finite sum of characteristic functions in the spatial domain. For example, $f$ could be a block of explosives in airport baggage or a sample of narcotic powder (e.g. fentanyl) in mail at customs.

Another major component of our regularization penalty is the use of ideas from multibang control [20–22]. The multibang penalty [20] is used to enforce solutions whereby the function outputs are constrained to a finite set. We aim to apply multibang ideas here to enforce the piecewise constant structure of $f$ in $(x_1,x_2,x_3)$ space, as discussed in the last paragraph. We do this by defining a set of binary switches $a_j$ to either activate or deactivate a characteristic function from our library. The multibang penalties are applied to enforce binary solutions for the $a_j$. So the finite set of solutions for the fitted $a_j$, in our case, is $\{0,1\}$. As our proposed objective function has smoothly defined gradients, we seek a relaxed alternative to the multibang penalty (the multibang penalty of [20] is not smooth), which we introduce later in section 2.3.

In addition to CS and multibang techniques, we also employ filtering ideas from microlocal analysis [6,23], to suppress the boundary artifacts typically observed in BST reconstruction, e.g., as observed in [24] in reconstructions from Bragg integral data. The filtering techniques from the literature are shown in section 2.5 to offer significantly improved image quality and artifact reduction in reconstructions from Monte Carlo data.

The remainder of this paper is organized as follows. In section 2 we explain our methodology. This includes a review of the physical model and Bragg transform from [1] in section 2.2, before moving on to explain our new reconstruction method in section 2.3. The reconstruction technique we propose is formalized as an algorithm in section 2.4, and the pre-processing using microlocal filters is explained in section 2.5. In section 3 we present our results on a wide variety of Monte Carlo and analytic simulations of interest in threat detection, and give comparison to a TV regularized solution.

2. Methodology

In this section we introduce our reconstruction technique and explain the filtering techniques used as data pre-processing from microlocal analysis. First we review the sensing geometry and Bragg transform of [1], and introduce some notation.

2.1 Sensing geometry

The scanner of Fig. 1 is equipped with linear detector collimation technology, which we will refer to as “Venetian blind" type collimation. The scanned object $f$ travels through the scanning tunnel (the red rectangle in Fig. 1(a)) in the $x_3$ direction on a conveyor belt, and is illuminated by a line of X-ray sources, located opposite a plane of detectors. The scanner sources (with coordinate $\mathbf {s}$) are fixed and switched along $\{x_2=-w_{x_2},x_3=0\}$, and are assumed to be polychromatic 2-D fan-beam (in the $(x_1,x_2)$ plane) with opening angle $\beta$. The detectors (with coordinate $\mathbf {d}$) are assumed to be energy-resolved and lie on the $\{x_2=w_{x_2}\}$ plane, with small (relative to the scanning tunnel size) offset $\epsilon$ in the $x_3$ direction. The detectors are collimated to record photons which scatter on planes in $\mathbb {R}^{3}$, and the planes of collimation are orientated to intersect the source $(x_1,x_2)$ plane along horizontal lines (parallel to $x_1$). Hence the photon arrivals measured by the scanner detectors are scattered from horizontal lines embedded in the $(x_1,x_2)$ plane. An example $\epsilon$ is shown in Fig. 1(b), which maps to the line $L=\{x_2=0,x_3=0\}$ at the half way point.

2.2 Bragg transform

Let $\mathcal {I}=[-w_{x_2},w_{x_2}]$, let $\mathfrak {E}=[E_m,E_M]$ be the energy range, and let $\Phi : \mathcal {I}\to \Phi (\mathcal {I})$ be a diffeomorphic map from the scanned line profile $x_2$ to the detector array position $\epsilon$. Let $L^{2}_0(\Omega )$ denote the set of square integrable functions with compact support on $\Omega \subset \mathbb {R}^{3}$. Then the Bragg transform $\mathfrak {B}_{a}$ defines a mapping from the target $f\in L^{2}_0(\mathfrak {E}\times [-w_{x_1},w_{x_1}]\times \mathcal {I})$ to the Bragg scatter measured by the scanner of Fig. 1 [1, page 7]

The solid angle is

The polarisation factor $P(\theta )$ is given by

We consider the recovery of the 2-D functions $f(\cdot ,\cdot ,x_2,0)$ from $\mathfrak {B}_af(\cdot ,\cdot ,\cdot ,\Phi (x_2))$, for each $x_2\in \mathcal {I}$. That is we consider the slice-by-slice reconstruction of $f$ from the 4-D data $\mathfrak {B}_af$. We focus exclusively on 2-D reconstruction here. This is to say that we do not consider the piecing together of the 2-D slices to form a 4-D image. This is left to future work. With this in mind we adopt the short-hand notation $f(q,x_1)=f(\cdot ,\cdot ,x_2,0)$, for some fixed $x_2\in \mathcal {I}$.

The operator $\mathfrak {B}_a$ is the same as considered in [1, Eq. (3.13)], but with the attenuation terms $A_1$ and $A_2$ removed from the modeling (using the notation of [1]). We remove $A_1$ and $A_2$ from the model of [1] as, in this paper, we assume no prior knowledge of the attenuation map or transmission data, and want to investigate what can be done with diffraction data alone. Further, as discussed in [1], the neglection of attenuation effects from the modeling is needed to prove linear invertiblity. Later, in the simulations conducted in section 3.6, the attenuation of the incoming and scattered rays is accounted for in the Monte Carlo data generation. Thus, there is a systematic error in the model due to the neglection of attenuation effects in (2.1).

In our case, $\mathfrak {B}_a$ reduces to the offset Bragg transform $\mathfrak {B}_{\epsilon }$ of [1, Eq. (6.2)] when $s_1=d_1$, with weight $W$ (using the notation of [1, Eq. (6.2)]) consisting of the physical modeling terms of Eq. (2.1). In [1, Theorem 6.1] $\mathfrak {B}_{\epsilon }$ is proven to be injective, and hence, trivially, $\mathfrak {B}_a$ is injective (i.e., $\mathfrak {B}_af=0\implies \mathfrak {B}_{\epsilon }f=0\implies f=0$), and we can solve for $f$ uniquely. Note that, in [1, Corollary 4.3], the physical modeling terms of (2.1) are shown to satisfy the conditions on $W$ in [1, Theorem 6.1]. Also, the set of Bragg curve integrals used in our simulations is sufficient to satisfy the saddle point conditions illustrated in [1, Fig. 7]. In all simulations conducted in this paper the conditions of [1, Theorem 6.1] are satisfied and uniqueness is guaranteed.

2.3 Reconstruction method

Throughout this paper, $f(q,x_1)$ (when discretized) will be represented as an $n\times m$ image, with $n$ the number of $q$ samples, and $m$ the number of $x_1$ samples. Let us fix $x_2\in \mathcal {I}$, and let $A\in \mathbb {R}^{p\times (mn)}$ denote the discretized form of the linear operator $\mathfrak {B}_a$, where $p$ is total number of data points. See Fig. 2 for a visualization of the discretized Bragg operator in $(q,x_1)$ space.

 

Fig. 2. Bragg curve examples for varying $E$, $s_1$ and $d_1$. $x_2=0$ is fixed. The figures were formed by reshaping the rows of $A$ into logical 2-D images. The curve in the left-hand figure is chosen so that $s_1=d_1$, as is considered in [1]. We see that the left-hand curve has the same shape as those shown in [1, Fig. 7].

Download Full Size | PDF

The regularization method we propose is derived from a set of a-priori assumptions regarding the target function $f$. We assume that $f$ is of the form

We introduce a library (finite set) of characteristic functions $\chi _1,\ldots ,\chi _l$ from which $f$ can be formed. The $a_j\in \{0,1\}$ act as binary switches to either activate or deactivate characteristic $j$ from the library. The $\chi _j$ library is chosen to comprehensively cover the image domain, but to also be restrictive enough to sufficiently regularize the solution. That is, we choose the smallest $l$ such that $\left [-\frac {2}{3}w_{x_1},\frac {2}{3}w_{x_1}\right ]$ (i.e., the width of the scanning tunnel in Fig. 1(a)) is adequately covered by the $\chi _j$. The set of $\chi _j$ used later in section 3., in our simulations, is shown to be small enough to offer high levels of regularization, while maintaining sufficient coverage of $\left [-\frac {2}{3}w_{x_1},\frac {2}{3}w_{x_1}\right ]$, and we see good results, for the most part, on phantoms which are comprised of characteristic functions lying outside of the chosen library.

Let the $\chi _j$ be sampled uniformly on $\left [-\frac {2}{3}w_{x_1},\frac {2}{3}w_{x_1}\right ]$, with $m$ samples, to form the vector $\mathbf {z}_j$. let $C_j=\left ({\mathbf {z}_{j1}I_n,\ldots ,\mathbf {z}_{jm}I_n}\right )^{T}$, where $\mathbf {z}_{ji}$ is the $i^{\textrm {th}}$ entry of $\mathbf {z}_j$ and $I_n$ is the $n\times n$ identity matrix. Then we define $A_j=AC_j$ as the restriction of $A$ to characteristic $\chi _j$. We define $G$ as the matrix with entries

Then $G_{ij}\neq 0$ if $\chi _i$ and $\chi _j$ intersect and $G_{ij}= 0$ otherwise. We are now ready to define our objective function.

We propose to minimize the functional

The term

To minimize $\mathcal {C}$ we use the L-BFGS-B method of [27] with box contraints, setting non-negativity constraints on the $\mathbf {y}_j$ and constraining $a_j\in [0,1]$. A solution is obtained by first fixing $\mathbf {a}$, and solving (2.12) for $\mathbf {y}$ for a number of inner iterates $n_2$, then fixing the $\mathbf {y}$ output and solving for $\mathbf {a}$ over $n_2$ iterates. The entire process is then repeated for $n_1$ outer iterations until convergence is reached. A more detailed explanation of the reconstruction algorithm is given below.

2.4 Algorithm: 2-D Bragg scatter reconstruction (2DBSR)

Initialize $b$, $\alpha$, $\lambda$, $\gamma$, $n_1$ and $n_2$. Initialize $\mathbf {a}^{0}$, $\mathbf {y}^{0}$, and characteristic library $(\mathbf {z}_1,\ldots ,\mathbf {z}_l)$. Define $G$ as in (2.11). Then for $n_{\textrm {outiter}}=1$ to $n_1$ do:

2.5 Artifact reduction and data pre-processing using microlocal filters

Here we discuss the idea to apply filtering techniques from microlocal analysis to the Bragg data, with the aim to suppress boundary type artifacts in the reconstruction, e.g., as are discovered in [24] in reconstructions from Bragg curve integral data.

In the literature on microlocal analysis in 2-D limited data X-ray CT [6,23], the authors make use of smoothing filters to reduce streaking artifacts which appear along straight lines at the boundary of the data set. Let $R$ denote the classical Radon transform and let $S\subset S^{1}\times \mathbb {R}$ be the subset of sinogram space for which $Rf\in L^{2}(S^{1}\times \mathbb {R})$ is known, where $f$ is the reconstruction target. Then the reconstruction in limited data CT is typically obtained by direct filtered backprojection

In our case the Bragg data $\mathfrak {B}_af$ has domain $D=[E_m,E_M]\times \Omega \times \Phi (\mathcal {I})$, where $\Omega =[-w_{x_1},w_{x_1}]^{2}$ is the space of source and detector $x_1$ coordinates, and $\Phi (\mathcal {I})$ is the set of detector offsets $\epsilon$ (refer to Fig. 1). Thus, the boundary of the sinogram in our case is a 3-D cube in energy, source position and detector position, for each $x_2$ considered. The full data consists of all $(E,s_1,d_1)\in \mathbb {R}_+\times \mathbb {R}^{2}$. Due to the finite scanner width and limitations on the energy range however, we consider the limited data on $[E_m,E_M]\times \Omega$, for each fixed $x_2$. To deal with the cutoff at the boundary of $[E_m,E_M]\times \Omega$ we construct a filtering function $\psi \in C^{\infty }_p\left ({D}\right )$ which goes to zero smoothly in the energy variable $E$, as $E\to E_M$. In the remaining variables $s_1$ and $d_1$, we see a natural smoothing effect in the data towards zero for $s_1, d_1$ near the $x_1$ limit of the portal scanner (i.e. at $x_1=\pm w_{x_1}$) due to solid angle effects and the reduction in source intensity with increasing $x_1$ (see (2.6)). That is, with a wide enough scanner (or large enough $w_{x_1}$), the detectors with $d_1\to \pm w_{x_1}$ have a more restrictive view of the scatter than those $d_1\approx 0$, and thus measure less scatter due to lower solid angle. The source projections with $s_1\to \pm w_{x_1}$ produce less counts due to reduced illumination of the object (i.e. at the scanner edge much of the source fan-beam does not intersect the scanning region), and decreased input intensity at further distances from the source. Also, as $E\to E_m$, we see a natural smoothing of the signal towards zero due to increased effects due to attenuation and self-absorption at low $E$. The self-absorbed photons are those which are absorbed by a material (due to the photoelectric effect) after being scattered initially from the same material. While the natural smoothing does not reach zero exactly at the boundary (e.g. at $s_1, d_1= \pm w_{x_1}$), we found the smoothing sufficient to reduce significantly the presence of boundary artifacts in the reconstruction. Therefore, we require no additional smoothing in the data in the $s_1$ and $d_1$ variables, and in $E$ as $E\to E_m$.

With this in mind, we define $\psi$ as the third order polynomial in $E$

 

Fig. 3. Plot of $\psi$, as in Eq. (2.16), with $E_M=29$keV and $E_m=1$keV.

Download Full Size | PDF

To show the effectiveness of this idea, and how it works, refer to Fig. 4. Figure 4(a) shows a plot of Monte Carlo data $\mathfrak {B}_af$ of the scatter measured by the portal scanner from NaCl and Carbon (diamond structure) spheres. In this case $x_2=0$ and the spheres are centered along the line $\{x_3=0,x_2=0\}$. Figure 4(b) shows the same data filtered by $\psi$. To clarify, the plots in Figs. 4(a) and 4(b) are the vectorized forms of the full 3-D data cube (in $(E,s_1,d_1)$), where each energy bin on the plot contains the corresponding vectorized 2-D datasets in $(s_1,d_1)$. The $s_1$ and $d_1$ positions used are specified later in section 3.1. We present the plots in this way purely to show the effects of the $\psi$ filtering. We can see that the plot of Fig. 4(b) is that of Fig. 4(a) point multiplied by $\psi$, as in Fig. 3.

 

Fig. 4. Microlocal filtering examples. The unfiltered data in (a) falls to zero sharply at $E=29$keV, which causes boundary artifacts in the corresponding reconstruction in (d). To alleviate the sharp cutoff, and the artifacts, we point multiply the data by $\psi$ as a pre-processing step before reconstruction. See (b) for the filtered data, and (e) for the corresponding reconstruction. We see a significant reduction in artifacts using the microlocal filtering approach. The ground truth image is shown in (c), for reference. The key take-away from this result is that we can use microlocal filtering methods to reduce artifacts in the reconstruction, by smoothing points near the boundary of the data where we see sharp jumps in value.

Download Full Size | PDF

The reconstructions of $f(q,x_1,0)$ corresponding to Figs. 4(a) and 4(b) are shown in Figs. 4(d) and 4(e), using $\mathfrak {B}_af$ and $\psi \mathfrak {B}_af$ as input to Algorithm 2.4, respectively, with $n_{\textrm {outiter}}=1$ outer iterations performed. Later in section 3 we will show the results with higher $n_{\textrm {outiter}}$, and explain the imaging parameters in more detail. This example is included purely to illustrate the effectiveness of microlocal filters for our problem. In Fig. 4(d) we see significant image artifacts, which roughly follow the shape of the Bragg integration curves of Fig. 2. The left-hand curve in Fig. 2 corresponds to the case when $E=E_M$, that is, at the sinogram boundary in $E$. Thus we see artifacts in the reconstruction due to the sharp cutoff in sinogram at $E=E_M$. When the cutoff is smoothed by $\psi$, using $\psi \mathfrak {B}_af$ as input to Algorithm 2.4 in Fig. 4(e), we notice a significant reduction in the artifacts, and a marked improvement in image quality.

3. Testing and results

In this section we perform the analytic and Monte Carlo data testing of Algorithm 2.4 in a variety of imaging scenarios of interest in airport baggage screening and threat detection. First we discuss the imaging specifics (i.e. source/detector positions, source energies etc.) and then give simulated reconstructions from analytic and Monte Carlo data.

3.1 Imaging parameters

Let us consider the scanning geometry depicted in Fig. 1. Throughout the simulations conducted in this section we convert to a more practical measurement system, using units of millimeters (mm), and set $w_{x_1}=300$mm and $w_{x_2}=410$mm. $\Phi$ is given explicitly by the linear map

As discussed in section 2.3, we aim to recover $f(q,x_1)$ on the range [0,2] Å^-1$\times [-300,300]\textrm {mm}$. The maximum energy $E_M=29$keV is chosen to correspond with the maximum $q$ in the reconstruction space, namely $q_M=2$. That is,

The input distribution $I_0(E)\propto \sum _{j=1}^{29}\delta (E-j)$ is chosen to be uniform in this study. This choice is well founded since, given the separability of $I_0(E,x_1)$ (Eq. (2.6)), we can divide through by $I_0(E)>0$ in Eq. (2.1) to obtain the Bragg signal with uniform input spectra. We use 600 detector positions $d_1\in \{-300+(j-1) : 1\leq j\leq 600\}$mm spaced at 1mm intervals. A 1mm detector pitch is realistic given the pixel sizes of the energy sensitive detectors on the market [29–31] (e.g., the detectors in [29] have a 250$\mu$m pitch). Silicon semiconductor detectors [30,31] can offer high energy resolution (up to $\approx 0.1$keV) at a low energy range (0-30keV), which is sufficient for our problem. The specific choice of detector is the subject for future work. The main focus of this paper is the methodology, and showing how the proposed regularizers can deal with noise effects such as attenuation.

The $\epsilon$ coordinate of the detectors varies with the chosen $x_2\in [-410,410]$mm, and is determined by $\epsilon =\Phi (x_2)$ (Eq. (3.1)). The detectors are assumed to be energy-resolved with energy resolution 1keV, so we are able to distinguish between all 29 energies in our range $E\in [1,29]$keV. The number of source positions (31) is low here as we anticipate longer scanning times per source projection so as to allow for sufficient photon counts. The increased number of detector positions (600) and energies bins (29) are at no detriment to the scan time. Note that due to the delta-comb nature of $I(E)$, there is no overlap in the energy bins. In total, the number of data points (or the number of rows of $A$) is $p=31\times 29\times 600 =539400$. We sample $f$ as a $750\times 600$ (high resolution) image, with $n=750$ $q$ samples in $\{q=\frac {2(j-1)}{750} : 1\leq j\leq 750\}$Å^-1 and $m=600$ $x_1$ samples in $x_1\in \{-300+(j-1) : 1\leq j\leq 600\}$mm. The length of $\mathbf {y}$ (i.e. the number of reconstructed pixels) is therefore $n\times m=750\times 600=45000$ (also the number of columns of $A$).

We use the characteristic library consisting of $l=405$ elements, with widths $w_j\in [1,3]$cm at $5$mm intervals, and centers $x^{c}_j\in [-200,200]$mm at 5mm intervals. We are assuming $f(q,x_1)=0$ for $|x_1|>200$mm, i.e., $f$ is zero outside of the scanning tunnel depicted in Fig. 1. We also assume that $f$ is composed of crystallite samples with width not exceeding 3cm, e.g., this could be a small sample of narcotics hidden amongst clutter in a mail parcel.

3.2 Quantitative analysis and hyperparameter selection

As a quantitative measure of the accuracy of our results we use the edge $F_1$ score [32–34]. Let $f$ and $\hat {f}$ denote ground truth and reconstructed images respectively. We employ the code of [35] here, which detects large jumps in the gradient images (i.e. in $\nabla f$ and $\nabla \hat {f}$), as the basis for computing $F_1$. The edge $F_1$ score is a measure of how well we have recovered $\nabla f$ and the locations of the Bragg peaks (i.e. the large spikes in the Bragg spectra of Fig. 5). The edge $F_1$ score is chosen as a quality metric as the locations of the Bragg peaks are the most valuable feature for materials characterization (e.g., when compared to the peak magnitudes). In classical powder Bragg spectroscopy (e.g., in experiments such as Debye-Scherrer [36,37]), a similar analysis is performed, whereby the locations of the Bragg peaks are detemined by measuring the radii of rings scattered onto a plane (or by measuring the scattering angles at which peaks occur in the case of [36,37]). For example, in [36, Fig. 8], the powder photographs show the locations at which peaks in intensity occur. The peak locations are then used to detemine the crystal lattice parameters. The absolute magnitude of the peaks, in this case, is not of special importance with regards to materials characterization. See the discussions after the reference to Fig. 8 on [36, page 575]. Thus, we can interpret a high edge $F_1$ score value as yielding a high level of materials characterization, and vice versa.

 

Fig. 5. $\hat {F}(\cdot ,Z)$ curve plots for a variety of $Z$. The plots are normalized by $L^{\infty }$ norm (max value). We restrict the $q$ range to $[0,0.6]$Å^-1 ($1$Å^-1 corresponds to $\sim 12.4$keV) to better highlight the more significant (larger) Bragg peaks.

Download Full Size | PDF

Let TP, FN and FP denote the proportion of true positives, false negatives and false positives in a classification result, respectively. Then the $F_1$ score is defined

Equivalently, the $F_1$ score is the geometric mean of the recall and precision of a classification result. The edge $F_1$ score is the $F_1$ score of (3.2) calculated between the edge maps $\nabla f$ and $\nabla \hat {f}$. Specifically, we use the code of [35] (with high gradient threshold) to calculate $\nabla f$ and $\nabla \hat {f}$ and convert to binary images. Then, we vectorize said binary images and use as input to Eq. (3.2). The gradient threshold used is kept fixed throughout the simulations conducted in this paper.

In addition to the edge $F_1$ score, we also quantify our results using the normalized correlation [7,8] (in those papers the normalized correlation is used to quantify the results, and for classification). Similarly to [7,8], we consider the normalized correlation of 1-D line profiles of the image. Let $f(\cdot ,x_1)$ and $\hat {f}(\cdot ,x_1)$ denote vertical line profiles, at $x_1$, of the ground truth and reconstructed images respectively. Then, the normalized correlation is defined [8, Eq. (11)]

If $f(\cdot ,x_1)=0$, or $\hat {f}(\cdot ,x_1)=0$, then we set $s(x_1)=0$. The normalized correlation provides a measure of how well the magnitudes of the Bragg peaks (in particular the larger Bragg peaks) match in $f(\cdot ,x_1)$ and $\hat {f}(\cdot ,x_1)$ (up to scaling). Let $[-w+x_c,w+x_c]\subset [-300,300]$mm be a subinterval of the imaging domain, where $x_c$ is the center of the interval, and $2w$ is the width. Then we define

Refer back to Eq. (2.12) and the definitions of the hyperparameters $\alpha ,\gamma ,\lambda$. Through experimentation we found that the hyperparameter values $\alpha =10^{6}$ (multibang parameter) and $\gamma =10^{10}$ (no overlap constraint) worked well for the examples considered here. We choose the smoothing parameter $\lambda$ ($L^{1}$ smoothing) for the best results in terms of edge $F_1$ score.

3.3 Materials considered

Let $Z$ denote the atomic number of the material. Then we consider the $F(q,Z)$ curves for three materials here, namely Carbon with a graphite structure (denoted C-graphite), salt and Carbon with a diamond structure (denoted C-diamond). We choose these three materials as they fall into our $Z<20$ range, and exhibit a variety of crystalline structures (e.g. hexagonal for graphite, and face-center cubic for salt), and are thus suitable for testing our algorithm. Here $Z$ replaces $\mathbf {x}\in \mathbb {R}^{2}$ from section 2.2 to represent the material. Furthermore, the $F(q,Z)$ curves for C-graphite, salt and C-diamond are well known and readily available in the literature (e.g. see [38] for the crystal structure of salt). Calculating $F(q,Z)$ for general compounds is a difficult task and there is work to be done in the spectroscopy literature to compile a wider database of $F$ curves.

The Bragg differential cross section $F(q,Z)$ is defined [1,39,40]

3.4 Comparison to other methods

For comparison, we consider a TV regularized solution, and use the microlocal filtering techniques explained in section 2.5 as data pre-processing to remove boundary artifacts. Specifically we aim to find

3.5 Analytic data testing

Here we consider the analytic data testing of Algorithm 2.4 by means of a Poisson noise model. That is, we use a scaled, matched model as mean to a multivariate Poisson distribution to generate data, where the scaling factor controls the level of noise. We consider two imaging phantoms for reconstruction. See Fig. 6. The left-hand phantom is comprised of an NaCl and C-graphite object with centers $x_c=-75$mm and $x_c=75$mm respectively, both with width $w=20$mm. The widths and centers used are chosen to be part of the imaging basis used (as detailed in section 3.1). The right-hand phantom is comprised of two NaCl objects with centers $x_c=-72.5$mm and $x_c=77.5$mm, both with width $w=1.25$cm. The widths and centers in this case are chosen to lie outside of the imaging basis. Let $\mathbf {y}$ denote the vetorized form of $f(\cdot ,\cdot ,x_2)$ for a given $x_2\in [-410,410]$mm. Then the phantoms are normalized to have max value $\|\mathbf {y}\|_{\infty }=1$.

 

Fig. 6. Imaging phantoms.

Download Full Size | PDF

The analytic data is generated by

The noisy data $\mathbf {b}_{\eta }$ is generated as a Poisson draw. The scaling parameter $\eta _c$ is the average counts per detector and controls the level of noise, i.e., larger $\eta _c$ implies less noise and vise-versa. We consider two average count levels $\eta _c=10$ and $\eta _c =1$, which approximately equate to relative least square errors of $\eta _{\textrm {ls}}=15\%$ and $\eta _{\textrm {ls}}=50\%$ respectively, where $\eta _{\textrm {ls}}=\frac {\|\mathbf {b}-\mathbf {b}_{\eta }\|_2}{\|\mathbf {b}\|_2}$. We consider three imaging line profiles $x_2\in [-410,410]$mm, namely $x_2=0$mm, $x_2=-205$mm, and $x_2=205$mm. Note that the operator $A$ varies with $x_2$, so we are considering three 2-D linear inverse problems with three different $A$ operators in total.

See Tables 1 and 2 for our results using 2DBSR (Algorithm 2.4) and FTV in terms of edge $F_1$ score. In Figs. 7 and 8 we present 2-D image reconstructions along the central line profile ($x_2=0$), for $\eta _c=1$. The image quality and $F_1$ scores are good and comparible using both methods. The average normalized correlation scores using 2DBSR and FTV are presented in Tables 3 and 4. The normalized correlation results show that FTV has higher performance in terms of resolving the magnitudes of the Bragg peaks (in particular the larger Bragg peaks) in the reconstruction. We notice a significantly increased $F_1$ score using FTV on phantom 2. This is because the phantom 2 materials were chosen to lie outside imaging basis used for 2DBSR. In this case the 2DBSR algorithm chooses the centers and widths closest to the ground truth, namely $x_c=\pm 75$mm and $w=1.5$cm. Thus in the 2DBSR reconstruction the centers and widths are slightly shifted from the ground truth causing a reduction in $F_1$ score. We could increase the size of the characteristic library to combat this, although this would be at a cost to the level of regularization and the implementation speed. When using FTV, we have no such restrictions on imaging basis. In general, with a matched model, and low level of noise, FTV offers an improved performance over 2DBSR. This makes sense, as with a low level of noise, the more powerful regularizers of 2DBSR are not necessary, and cause a higher level of regularization error than FTV.

 

Fig. 7. Phantom 1 reconstructions using 2DBSR (left) and FTV (right) along the central line profile ($x_2=0$), with count level $\eta _c=1$ ($50\%$ noise).

Download Full Size | PDF

 

Fig. 8. Phantom 2 reconstructions using 2DBSR (left) and FTV (right) along the central line profile ($x_2=0$), with count level $\eta _c=1$ ($50\%$ noise).

Download Full Size | PDF

Table 1. Phantom 1 edge $F_1$ score results. The measurements in the left column are given in mm.

View Table | View all tables in this article

Table 2. Phantom 2 edge $F_1$ score results. The measurements in the left column are given in mm.

View Table | View all tables in this article

Table 3. Phantom 1 average normalized correlation results. The measurements in the left column, and the subscripts in the second and third columns are in millimeters. We give $\textrm {R}_{x_c,w}$ for each pair $(x_c,w)$ where materials are located, for each experiment conducted. For example, $\textrm {R}_{-75,10}$ is the average normalized correlation corresponding to the NaCl phantom.

View Table | View all tables in this article

Table 4. Phantom 2 average normalized correlation results. The measurements in the left column, and the subscripts in the second and third columns are in millimeters. We give $\textrm {R}_{x_c,w}$ for each pair $(x_c,w)$ where materials are located, for each experiment conducted. For example, $\textrm {R}_{-72.5,6.25}$ is the average normalized correlation corresponding to the left-hand NaCl phantom.

View Table | View all tables in this article

3.6 Monte Carlo experiments

In the Monte Carlo tests conducted in this section, we aim to recover spherical crystalline objects which are centered on the central ($x_2=0$) line profile. Consider the schematic shown in Fig. 9. Here we have displayed the attenuation image ($\mu _E$ for $E=100$keV) of two spherical crystallites centered on the $x_2=0$ line profile, which are embedded in clutter (e.g. paper, clothing). The experimental setup of Fig. 9 is motivated by security screening problems. For example, Fig. 9 could simulate a parcel in commerce mail or a carry-on bag at an airport, whereby a small sample of crystalline material is concealed in amongst a large volume of clutter.

 

Fig. 9. Test parcel.

Download Full Size | PDF

We use the imaging parameters (e.g. source/detector positions, source energies) specified in section 3.1. The $\epsilon$ coordinate of the detector array is determined by $\Phi$ (Eq. (3.1)), and is fixed at $\epsilon _0=\Phi (0)=37.5$mm throughout this section. $\epsilon _0$ is the $\epsilon$ coordinate of the detector array which corresponds to $x_2=0$. As discussed in section 3.1, we use 29 source energies $E\in [1,29]$keV at 1keV intervals. The source intensity used here is $I_0(E)=6\times 10^{10}$ counts per projection, per energy, which amounts to $1.74\times 10^{12}$ counts per projection over all 29 energy bins. Equivalently the source intensity is $2.94\times 10^{9}$ counts, per cm$^{2}$, per energy bin, at 1m from the source. To achieve such intensity, we estimate a source projection time in the order of seconds (e.g., using an X-ray tube with low filtering, so as to allow for low energies). See Fig. 10. Here we have shown an example X-ray tube spectra generated using the Spekcalc software [43–45]. In this case, the source projection time would be $\approx 2$ seconds to achieve $2.94\times 10^{9}$ counts per energy bin at 1m, yielding a full scan time, over the 31 source projections used, in the order of minutes (in this example the total scan time would be $\approx 1$ minute). A more precise estimate of the scan time is the subject of future work, once the source tube specifics (e.g. voltages, ampage) are decided upon. Further, in this paper we wish to investigate what happens when we allow for long scanning times, before moving onto the more challenging problem of high-speed scanning in future work.

 

Fig. 10. Example unfiltered 50kV Tungsten target X-ray tube spectra, with 10mA current. Source spectra generated using Spekcalc software [43–45].

Download Full Size | PDF

The clutter is simulated as a 30cm$^{3}$ box centered at the origin. The clutter attenuation coefficient is simulated as one-tenth that of cellulose. Cellulose is chosen as a material as it is a primary component of paper and clothing, which are commonly found in airport luggage and mail parcels. The attenuation coefficient of the cellulose is scaled by $1/10$ following the assumption that the clutter will largely be comprised of air, thus reducing the attenuative effects. The simulated clutter is an amorphous material with significantly lower attenuation coefficient than that of the crystalline materials considered in this paper (i.e. C-graphite, NaCl, C-diamond). Due to the linear collimation technology of the detectors, and fan-beam source, the scatter from the clutter occurs along 1-D lines in the image space. The scatter from a 1-D line of material with 1/10 the attenuation coefficient of Cellulose is negligible, especially when we consider the relatively low range of energies used (i.e., $E\in [1,29]$keV) and self-absorption effects. Hence, the scattering contribution from the clutter is modeled as zero. The role of the Cellulose is to simulate the attenuation, and signal reduction effects of clutter. The attenuation of the incoming and scattered rays though the full 3-D volume of clutter is included in the Monte Carlo data generation (this is a non-negligible effect).

Let $x_c\in [-200,200]$mm denote the sphere center $x_1$ coordinate. The sphere center $x_2$ and $x_3$ coordinates are fixed at $x_2=x_3=0$mm. Let $r\in [5,15]$mm denote the sphere radius. For example, the spheres shown in Fig. 9 have center $x_1$ coordinates $x_c=\pm 75$mm and radius $r=15$mm. Then we consider nine Monte Carlo experiments in total which are designed to cover a range of $x_c$ and $r$, for materials contained inside and outside the imaging basis used (as detailed in section 3.1).

To generate the data, we use a novel, single scatter Monte Carlo (MC) code developed in Matlab, which we introduce in this paper. A step-by-step description of the MC code and a pseudo code is provided in the appendix. The full code is available from the authors upon request. While it is ideal to test with multiple scatter, here, since the samples are small ($<3$cm in dimension), and given the type of detector collimation used (shown in Fig. 1), multiple scatter will be negligible. In [14], the collimator blades are used in a similar way to suppress the multiple scatter (see the discussions at the end of page 2467). Further, the experiments conducted here are initial tests of the algorithm to show how our regularizers can deal with the more major mismatch between the physics used to generate the data and the physics employed by the inversion model; namely attenuation effects and Compton scatter.

Our code can provide an easier but less accurate simulation when compared to other scatter simulation packages on the market, e.g., such as GEANT4 [46], for researchers who are looking to go beyond a Poisson or additive Gaussian noise model. Some researchers may also prefer the Matlab implementation of our code, compared to GEANT4, which is written in C++. GEANT4 has a steep learning curve, and requires a high level of training to use, but is more accurate than our code.

Our code fully incorporates unmodeled effects in the data generation, including attenuation and Compton scatter. We do not assume knowledge of the attenuation coefficient in our model (i.e., the model of (2.1)), and thus there is systematic error in the model due to the neglection of attenuation modeling. In the energy and scattering angle range considered ($E<30$keV and $\omega <114.6^{\circ }$) the Compton scatter intensity is significantly lower than that of the Bragg signal. Hence the Compton scatter will act as an additional background noise in our simulations.

In Fig. 11 we have plotted the systematic attenuation error and Compton signal for a single sphere of Carbon-graphite with center $x_c=0$, for varying $r$ and $E$. We see in Fig. 11(a) that as $r$ increases and $E$ decreases (towards the bottom-right corner of the figures) the systematic error due to attenuation increases due to higher self-absorption effects at low energy (fewer photons penetrate) and higher r (those that suffer greater attenuation due to Beer’s Law). We see the converse effect for the Compton noise in Fig. 11(b) as much of the Compton scatter is self-absorbed for low $E$ and high $r$. This is as we’d expect and Fig. 11(a) appears as the negative image of Fig. 11(b) (after scaling and translation). So there is a trade-off to consider here, and no size of object or photon energy is necessarily preferred over the other. The MC experiments considered here are designed to reflect a variety of levels of Compton noise and attenuation error in the data (i.e. we consider a variety of $r$ and $E$).

 

Fig. 11. Plots of relative, systematic attenuation error (left) and absolute Compton signal (right) with $E$ and sphere radius ($r$).

Download Full Size | PDF

To obtain the 2DBSR reconstructions presented in this section we implement the following two-stage reconstruction procedure. First, we apply Algorithm 2.4, with the filtered Bragg data $\mathbf {b}$ (using the filters specified in section 2.5), characteristic library and hyperparameters (as specified in section 3.1) as input, choosing $\lambda$ for the best results in terms of edge $F_1$ score. Then, we perform a second full pass of Algorithm 2.4, replacing $\lambda \to \lambda \times 10$. We found that increasing $\lambda$ by an order of magnitude on the second pass of Algorithm 2.4 helped to reduce background noise in the reconstruction, once the more significant artifacts (e.g. boundary artifacts) had been removed in the first pass of Algorithm 2.4. To show how this works, see Fig. 12 for image reconstructions of an NaCl sphere, center $x_c=0$mm, radius $r=15$mm, at each stage of the two-stage process. We notice an improvement in the image quality and $F_1$ score as we advance from stage 1 to stage 2.

 

Fig. 12. Illustration of two-stage reconstruction process used for 2DBSR. We show reconstructions of a C-diamond sphere, center $x_c=-6.25$mm and radius $r=8.75$mm. at each stage of the reconstruction process. Figures (a)-(c) show the ground truth and the 2-D reconstructions, and figures (d) and (e) show the central line profiles of the 2-D reconstructions corresponding to each stage. The $F_1$ scores and normalized correlations corresponding to each stage are given in the subfigure captions.

Download Full Size | PDF

Upon applying the same two-stage process to the FTV reconstructions, we did not see an improvement in the results, and hence why we did not apply the two-stage process to FTV. This is as expected since the objective of (3.8) is convex, and hence, as we chose $\lambda$ for FTV for the best results, increasing $\lambda \to \lambda \times 10$ gave the same result (i.e. the best result for FTV over all $\lambda$).

In Table 5, we have presented our results in terms of edge $F_1$ score for the nine MC experiments conducted. In Table 5 we use the shorthand notation “St" for NaCl (salt), “C" for C-graphite and “Di" for C-diamond. Experiment $i$. for $1\leq i\leq 9$, is denoted by $\textrm {Exp}_i$. In brackets we specify the centers and radii of the spheres in millimeters. For example, “St ($x_c=0$, $r=10$)" means a single NaCl sphere with center $x_c=0$mm and radius $r=10$mm, and “St ($x_c=-75, r=5$), Di ($x_c=75, r=5$)" means an NaCl and C-diamond sphere with center $\pm 75$mm and radius $r=5$mm. Stages 1 and 2 of 2DBSR detailed above are denoted by S1 and S2 in Table 5. We see an improved performance at all stages of 2DBSR over FTV, in all experiments conducted. We found that S2 offered better image quality and $F_1$ score overall, when compared to S1.

Table 5. $F_1$ scores for each experiment conducted, comparing stage 1 (S1), and stage 2 (S2) of 2DBSR (ours) with FTV (literature). Experiment $i$, for $1\leq i\leq 9$, is labeled by $\textrm {Exp}_i$. The measurements in the brackets in the left-hand column are in millimeters, and describe the centers ($x_c$) and radii ($r$) of the spheres scanned. Here “St" denotes NaCl (salt), “C" denotes C-graphite and “Di" is C-diamond. The asterisked entries denote materials which are outside of the imaging basis. The bold rows correspond to the image reconstructions shown in Figs. 13–17.

View Table | View all tables in this article

In Table 6, we present our results in terms of the average normalized correlation ($\textrm {R}_{x_c,r}$), for all experiments conducted. We report $\textrm {R}_{x_c,r}$ for each pair $(x_c,r)$ where materials are located. For example, in experiment 6, $\textrm {R}_{-75,10}$ is the average normalized correlation corresponding to the NaCl phantom. In the majority of cases, we see a significantly improved performance, in terms of $\textrm {R}_{x_c,r}$, at all stages of 2DBSR over FTV. The normalized correlation scores at S1 and S2 are more comparable overall, when compared to the $F_1$ scores.

Table 6. Average normalized correlations ($\textrm {R}_{x_c,r}$) for each experiment conducted, comparing stage 1 (S1), and stage 2 (S2) of 2DBSR (ours) with FTV (literature). Experiment $i$, for $1\leq i\leq 9$, is labeled by $\textrm {Exp}_i$, as in Table 5. The measurements in the subscripts in the second, third, and fourth columns are in millimeters. We give $\textrm {R}_{x_c,r}$ for each pair $(x_c,r)$ where materials are located, for each experiment conducted. For example, in $\textrm {Exp}_6$, $\textrm {R}_{-75,10}$ is the average normalized correlation corresponding to the salt phantom. The bold rows correspond to the image reconstructions shown in Figs. 13–17, as in Table 5.

View Table | View all tables in this article

In Tables 5 and 6 we have highlighted five rows in bold. These correspond to the image reconstructions in Figs. 13–17, which are chosen to highlight the best and worst results using our method (at stage 2), for materials within and outside of the imaging basis. In Figs. 13–17, “GT" denotes the “Ground Truth", and “LP" denotes “Line Profile". We specify the line profile equation in brackets in each case. In Fig. 13 we have presented image reconstructions of a C-graphite sphere, center $x_c=0$mm and radius $r=15$mm. In this case the attenuation error is high and the Compton noise is low (see Fig. 11). We see structured distortions (systematic error) of the GT Bragg peaks in the FTV reconstruction. That is, the horizontal line segments of the GT are stretched in the positive $q$ direction, and we see a shifting of the Bragg peaks in the $x_1=0$mm LP. Given the assumed structure of $f$ when using 2DBSR (as in Eq. (2.10)), the 2DBSR algorithm does not allow for such structured distortions in the reconstruction. Hence we see a complete removal of the artifacts in the 2DBSR reconstruction, with an increase in $F_1$ and normalized correlation score as a result. In Fig. 14 we see a similar effect, although the $F_1$ score and normalized correlation using 2DBSR is decreased as, in this case, the C-diamond sphere is out of basis. The $F_1$ score and normalized correlation using FTV is improved however, and the shifting of the Bragg peaks is reduced in the reconstruction. This is because the attenuation error is decreased due to the smaller sphere radius ($r=8.75$mm).

 

Fig. 13. C-graphite sphere reconstruction, $x_c=0$mm, $r=15$mm.

Download Full Size | PDF

 

Fig. 14. C-diamond sphere reconstruction, $x_c=-6.25$mm, $r=8.75$mm.

Download Full Size | PDF

 

Fig. 15. NaCL and C-diamond sphere reconstruction, centers $x_c=-75$mm and $x_c=75$mm respectively, both with $r=10$mm.

Download Full Size | PDF

 

Fig. 16. NaCL and C-graphite sphere reconstruction, centers $x_c=-75$mm and $x_c=75$mm respectively, both with $r=10$mm.

Download Full Size | PDF

 

Fig. 17. Two NaCl spheres reconstruction, centers $x_c=-72.5$mm and $x_c=77.5$mm, both with $r=6.25$mm.

Download Full Size | PDF

In Fig. 15, we have presented reconstructions of NaCL and C-diamond spheres within the imaging basis, with centers $x_c=-75$mm and $x_c=75$mm respectively, both with $r=10$mm. Here we see some significant noise effects in-between the Bragg peaks in the 2DBSR reconstruction, at a detriment to the $F_1$ score ($F_1=0.65$), and the normalized correlation of the right-hand phantom ($\textrm {R}_{75,10}=0.58$), as some of the noise is mistaken for additional (unwanted) Bragg peaks. We may be able to combat this by continuing to increase $\lambda$ (the $L_1$ regularization parameter) and running further iterations of 2DBSR as in the second stage of the two-stage reconstruction process outlined above, with the idea that the sparse regularizers would remove the noise and focus on the more significant (true) Bragg peaks in the reconstruction. The centers and widths of the spheres in the 2DBSR reconstruction are recovered correctly, and the $F_1$ score is acceptable. Furthermore, the results using 2DBSR are much improved over FTV, where we see severe artifacts in the reconstruction, with low $F_1$ score ($F_1=0.32$), and normalized correlation scores ($\textrm {R}_{-75,10}=0.19$ and $\textrm {R}_{75,10}=0.29$).

In Fig. 16 we have presented $f(q,x_1)$ reconstructions of a NaCl and C-graphite sphere with centers $x_c=-75$mm and $x_c=75$mm, repsectively, both with radius $r=10$mm. We see that the NaCl sphere is effectively zeroed out in the 2DBSR and FTV reconstruction, leading to a low $F_1$ score and normalized correlation using both methods. The reason for the poor image quality is that the scattered signal contribution from the Carbon-graphite sphere is more than twice that of the salt sphere. The sum of counts (taken over all $p=539400$ data points) scattering from the salt sphere is $7.1\times 10^{5}$, and the sum of the counts from the Carbon-graphite sphere is $1.7\times 10^{6}$. The NaCl scatter therefore accounts for $29\%$ of the total signal. As the signal from the NaCl sphere is largely hidden by that of the Carbon, this makes the NaCl $f(q,x_1)$ function more difficult to recover, and it is set close to zero by 2DBSR and FTV. Similar effects are seen in conventional X-ray CT, for example, in head scans, where the bone and skull shows up strongly in the image (due to higher absorption) and the skull interior and soft tissue are harder to see [26, Fig. 1].

In Fig. 17 we show $f(q,x_1)$ reconstructions of two NaCl spheres with centers $x_c=-72.5$mm and $x_c=77.5$mm, respectively, both with radius $r=6.25$mm. Thus, in this example, the NaCl spheres are outside the 2DBSR imaging basis. In this case we notice that the Bragg peaks of the right-hand NaCl sphere are blurred in the 2DBSR reconstruction. The blurring effect is mild however and the image quality is improved when compared to FTV. In the FTV reconstruction the shifting of the Bragg peaks, although still prominent, is less significant than in Figs. 13–16. We notice, in general, a reduction in the Bragg peak shift in the FTV reconstructions with decreasing $r$ (i.e. as the attenuation error decreases). This observation is supported by the results of Tables 5 and 6, as the $F_1$ and normalized correlation scores using FTV increase (and become more competitive with 2DBSR) with decreasing $r$. Thus it seems that the main cause of error in the FTV reconstructions is due to attenuation modeling, based on the error analysis of Fig. 11.

4. Conclusions and further work

Here we have presented the 2DBSR algorithm - a novel reconstruction technique for two-dimensional Bragg scatter imaging. The regularization penalties we applied are based on ideas in multibang control and compressive sensing. We also incorporated filtering ideas from microlocal analysis for data pre-processing, with the aim to suppress boundary type artifacts (e.g. such as those observed in [24]) in the reconstruction. In section 2.3 we formalized our approach and detailed the 2DBSR algorithm in section 2.4. The microlocal filters used for data pre-processing were defined as polynomials in $E$ in section 2.5. In section 3. we designed and conducted a variety of BST experiments, using analytic and Monte Carlo data, which are of interest in airport baggage screening and threat detection. Here we compared the performance of 2DBSR to a Filtered Total Variation (FTV) approach. In the analytic data experiments both methods performed well, with FTV slightly outperforming 2DBSR in terms of $F_1$ score, in most cases. In the Monte Carlo experiments, 2DBSR was shown to offer significantly higher performance than that of FTV in all cases considered (nine MC experiments in total). In the FTV reconstructions we saw severe image artifacts due to errors in the attenuation modeling. The artifacts appeared in the image as a shifting/distorting of the ground truth Bragg spectra in the positive $q$ direction, which caused a reduction in $F_1$ score and normalized correlation. In the majority of cases, the 2DBSR algorithm was successful in combatting the image artifacts due to attenuation modeling (among other errors, such as those due to Compton scatter), and the $F_1$ score and normalized correlation was much improved over FTV.

In this paper, the hyperparameter $\lambda$ for FTV and 2DBSR was tuned to the best results in terms of $F_1$ score. We aim to consider automatic hyperparameter selection methods (e.g., cross validation) in future work to improve the robustness of our method. The parameters, $\alpha$, $\gamma$, and the gradient threshold used for calculating the $F_1$ score are kept fixed throughout the simulations conducted in this paper, which suggests some level of robustness with regards to these parameters in terms of the performance of our algorithm. In future work we aim to also consider the automatic tuning of $\alpha$, $\gamma$, and the gradient threshold to further improve the robustness and performance of our algorithm

In the simulations conducted here, the energy range $E\in [1,29]$keV was used, which is relatively low in terms of X-ray CT applications. This energy range was chosen so that our forward model satisfied the uniqueness conditions of [1]. Given the low energy range, we were constrained to imaging small sample (i.e., $<3$cm in dimension), low $Z$ (i.e., $Z<20$) material in this paper, as large, more dense objects would absorb many of the low energy photons. In further work, we aim to consider a higher energy range, whilst also retaining the uniqueness conditions of [1], so that the restrictions of our method to small, low $Z$ samples can be lifted. This will likely require a reduction in the detector offsets ($\epsilon$), so that the scatter angles ($\omega$) are more restricted to a lower range, e.g., as is done in CSCT [13,14,47], where $\omega < 10^{\circ }$ and higher energy photons are used to image large objects. To determine if such a reduction in detector offset is possible, more discussions with industry collaborators are required.

We have applied our method in this paper to benign, crystalline materials, e.g., NaCl. Further testing on a wider variety of materials of interest in threat detection (e.g., narcotics), and amorphous materials (e.g., water) is a core goal we aim to pursue in future work.

We aim to consider full 4-D Bragg imaging using low-rank tensor ideas [48]. This follows the idea that the neighboring 2-D slices (i.e., the $f(\cdot ,\cdot ,x_2,x_3)$ for neighboring $(x_2,x_3)$) are highly correlated. Specifically, we aim to minimize the nuclear norm (as used in [48] in multispectral X-ray CT) to enforce linear dependence among the 2-D image slices, and regularize the 4-D inversion.

We also aim to consider joint electron density/Bragg spectra reconstruction using a combination of Bragg and Compton scatter data. Here we aim to use the models from the Compton scatter tomography literature [49–55] to improve the forward model (2.1), and the recovery of $f(q,x_1)$ in $x_1$ space. That is, the Compton data will be primarily used for density recovery and locating the crystallites (i.e. recovering the centers and widths of the materials), with the recovery in $q$ space coming mainly from the Bragg data. We expect that such ideas will be effective in combatting the more significant errors observed in Figs. 16 and 17 in the 2DBSR reconstructions, where we saw the zeroing out of the NaCl sphere (in Fig. 16) and blurring of the ground truth Bragg spectra (in Fig. 17), which yielded poor $F_1$ scores.

Appendix: description of Monte Carlo code

Consider the scattering event pictured in Fig. 18. Let $L_{\mathbf {x}_1\mathbf {x}_2}=\{\mathbf {x}_1+t(\mathbf {x}_2-\mathbf {x}_1) : t\in \mathbb {R}\}$ be the line through $\mathbf {x}_1$ and $\mathbf {x}_2$, and let us consider the scattering interaction space $\Omega =\{\textrm {NI},\textrm {PE},\textrm {INC},\textrm {COH}\}$, where NI denotes “no interaction" (transmitted photons), PE denotes “photoelectric absorption", INC denotes “incoherent scatter", and COH denotes “coherent scatter". Let $\mu (E,Z)=\sum _{\mathcal {E}\in \Omega \backslash \{\textrm {NI}\}}\mu _{\mathcal {E}}(E,Z)$ denote the total attenuation coefficient, where $\mu _{\mathcal {E}}$ is the attenuation coefficient for the interaction type $\mathcal {E}\in \Omega \backslash \{\textrm {NI}\}$, and $Z$ is the atomic number. Then the Monte Carlo algorithm, for a single photon trial, scattering site ($\mathbf {x}$) and source position ($\mathbf {s}$), reads as follows:

• 1. Sample the initial photon energy $E$ from the source spectrum. The photon is emitted from $\mathbf {s}$ and travels in the direction $\mathbf {x}-\mathbf {s}$ towards $\mathbf {x}$.
• 2. Sample from a Bernoulli distribution, where $e^{-\int _{L_{\mathbf {s}\mathbf {x}}}\mu (E,Z)}$ is the probablity of success, to decide if the photon reaches $\mathbf {x}$.
• 3. Determine the photon interaction by sampling from a finite distribution on $\Omega$, with $P(\textrm {NI})=e^{-\mu (E,Z)}$ and

  (A.1)$$P(\mathcal{E})=(1-P(\textrm{NI}))\frac{\mu_{\mathcal{E}}(E,Z)}{\mu(E,Z)}$$

  for $\mathcal {E}\in \Omega \backslash \{\textrm {NI}\}$.
• 4. Sample $\omega \in [0,\pi ]$ (polar angle of the scatter direction) from the differential cross section distribution for the given interaction type, keeping a uniform spread over each scattering cone, namely the cone with central axis direction $\mathbf {x}-\mathbf {s}$ (the dashed line in Fig. 18), vertex $\mathbf {x}$ and opening angle $\omega$. That is, we sample the azimuth angle of the scatter direction (in $[0,2\pi ]$) from a uniform distribution.
• 5. If the scattered photon travels in the direction $\mathbf {d}-\mathbf {x}$ towards $\mathbf {d}$, sample from a second Bernoulli distribution, with $e^{-\int _{L_{\mathbf {x}\mathbf {d}}}\mu (E_s,Z)}$ as the probablity of success, to decide whether the photon reaches $\mathbf {d}$.
• 6. If the scattered photon reaches $\mathbf {d}$, count one photon with energy $E_s$, else continue. Repeat for all detectors in the array.

The psuedo code detailed above is applied to the geometry of Fig. 1, with the differential cross section distribution set to one of the Bragg spectra of Fig. 5, depending on the material. For Compton scatter, the differential cross section is set to the product of the Klein-Nishina distribution [56,57], and the incoherent scattering function [41]. To generate the data used in section 3.6, we expand the single trial, single scatter site code to multiple scattering events along lines parallel to $x_1$ (e.g. $L$ as pictured in Fig. 1).

 

Fig. 18. A scattering event occurs at a scattering site $\mathbf {x}$, for photons emitted from a source $\mathbf {s}$ and recorded at a detector $\mathbf {d}$ (displayed as a small square in the $x_1x_3$ plane). The initial photon energy is $E$ and the scattered energy is $E_s$. Here $\mathbf {v}$ is the direction normal to the detector surface. The scattering angle is $\omega$.

Download Full Size | PDF

Funding

U.S. Department of Homeland Security (18STEXP00001-03-02).

Acknowledgments

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 Awards 18STEXP00001-03-02, formerly 2013-ST-061-ED0001 and 70RSAT19FR0000155. 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.

Disclosures

The authors declare no conflicts of interest.

Data availability

Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.

References

1. J. W. Webber and E. L. Miller, “Bragg scattering tomography,” Inverse Probl. Imaging 15(4), 683–721 (2021). [CrossRef]  

2. S. Hashemi, S. Beheshti, P. R. Gill, N. S. Paul, and R. S. Cobbold, “Accelerated compressed sensing based CT image reconstruction,” Comput. Mathematical Methods Medicine 2015, 1–16 (2015). [CrossRef]  

3. X. Li and S. Luo, “A compressed sensing-based iterative algorithm for CT reconstruction and its possible application to phase contrast imaging,” Biomedical engineering online 10(1), 73 (2011). [CrossRef]  

4. Z. Zhu, K. Wahid, P. Babyn, D. Cooper, I. Pratt, and Y. Carter, “Improved compressed sensing-based algorithm for sparse-view CT image reconstruction,” Comput. Mathematical Methods Medicine 2013, 1–15 (2013). [CrossRef]  

5. G.-H. Chen, J. Tang, and S. Leng, “Prior image constrained compressed sensing (PICCS): a method to accurately reconstruct dynamic CT images from highly undersampled projection data sets,” Med. Phys. 35(2), 660–663 (2008). [CrossRef]  

6. L. Borg, J. Frikel, J. S. Jørgensen, and E. T. Quinto, “Analyzing reconstruction artifacts from arbitrary incomplete X-ray CT data,” SIAM J. Imaging Sci. 11(4), 2786–2814 (2018). [CrossRef]  

7. M. Hassan, J. A. Greenberg, I. Odinaka, and D. J. Brady, “Snapshot fan beam coded aperture coherent scatter tomography,” Opt. Express 24(16), 18277–18289 (2016). [CrossRef]  

8. J. A. Greenberg, K. Krishnamurthy, and D. Brady, “Snapshot molecular imaging using coded energy-sensitive detection,” Opt. Express 21(21), 25480–25491 (2013). [CrossRef]  

9. K. MacCabe, K. Krishnamurthy, A. Chawla, D. Marks, E. Samei, and D. Brady, “Pencil beam coded aperture X-ray scatter imaging,” Opt. Express 20(15), 16310–16320 (2012). [CrossRef]  

10. J. A. Greenberg, M. Hassan, K. Krishnamurthy, and D. Brady, “Structured illumination for tomographic X-ray diffraction imaging,” Analyst 139(4), 709–713 (2014). [CrossRef]  

11. J. A. Greenberg, K. Krishnamurthy, M. Lakshmanan, K. MacCabe, S. Wolter, A. Kapadia, and D. Brady, “Coding and sampling for compressive X-ray diffraction tomography,” in Wavelets and Sparsity XV, vol. 8858 (International Society for Optics and Photonics, 2013), p. 885813.

12. J. A. Greenberg, M. N. Lakshmanan, D. J. Brady, and A. J. Kapadia, “Optimization of a coded aperture coherent scatter spectral imaging system for medical imaging,” in Medical Imaging 2015: Physics of Medical Imaging, vol. 9412 (International Society for Optics and Photonics, 2015), p. 94125E.

13. D. L. Batchelar, M. T. Davidson, W. Dabrowski, and I. A. Cunningham, “Bone-composition imaging using coherent-scatter computed tomography: Assessing bone health beyond bone mineral density,” Med. Phys. 33(4), 904–915 (2006). [CrossRef]  

14. U. Van Stevendaal, J.-P. Schlomka, A. Harding, and M. Grass, “A reconstruction algorithm for coherent scatter computed tomography based on filtered back-projection,” Med. Phys. 30(9), 2465–2474 (2003). [CrossRef]  

15. G. Harding and J. Kosanetzky, “Elastic scatter computed tomography,” Phys. Med. Biol. 30(2), 183–186 (1985). [CrossRef]  

16. G. Harding, H. Fleckenstein, D. Kosciesza, S. Olesinski, H. Strecker, T. Theedt, and G. Zienert, “X-ray diffraction imaging with the multiple inverse fan beam topology: Principles, performance and potential for security screening,” Appl. Radiat. Isot. 70(7), 1228–1237 (2012). [CrossRef]  

17. G. Harding, “X-ray diffraction imaging—a multi-generational perspective,” Appl. Radiat. Isot. 67(2), 287–295 (2009). [CrossRef]  

18. L. Rencker, F. Bach, W. Wang, and M. D. Plumbley, “Sparse recovery and dictionary learning from nonlinear compressive measurements,” IEEE Trans. Signal Process. 67(21), 5659–5670 (2019). [CrossRef]  

19. E. Herrholz, D. Lorenz, G. Teschke, and D. Trede, Sparsity and Compressed Sensing in Inverse Problems (Springer International Publishing, Cham, 2014), pp. 365–379.

20. C. Clason and k. Kunisch, “Multi-bang control of elliptic systems,” Inverse Probl. Imaging 31(6), 1109–1130 (2014). [CrossRef]  

21. C. Clason and K. Kunisch, “A convex analysis approach to multi-material topology optimization,” ESAIM: Math. Modell. Numer. Anal. 50(6), 1917–1936 (2016). [CrossRef]  

22. C. Clason, F. Kruse, and K. Kunisch, “Total variation regularization of multi-material topology optimization,” ESAIM: Math. Modell. Numer. Anal. 52(1), 275–303 (2018). [CrossRef]  

23. J. Frikel and E. T. Quinto, “Characterization and reduction of artifacts in limited angle tomography,” Inverse Probl. 29(12), 125007 (2013). [CrossRef]  

24. J. W. Webber and E. T. Quinto, “Microlocal analysis of generalized Radon transforms from scattering tomography,” arXiv preprint, arXiv:2007.00208 (2020)

25. W. H. Bragg and W. L. Bragg, “The reflection of X-rays by crystals,” Proc. Royal Soc. London. Ser. A, Containing Pap. a Math. Phys. Character 88, 428–438 (1913). [CrossRef]  

26. M. J. Ehrhardt, K. Thielemans, L. Pizarro, D. Atkinson, S. Ourselin, B. F. Hutton, and S. R. Arridge, “Joint reconstruction of PET-MRI by exploiting structural similarity,” Inverse Probl. 31(1), 015001 (2015). [CrossRef]  

27. R. H. Byrd, P. Lu, J. Nocedal, and C. Zhu, “A limited memory algorithm for bound constrained optimization,” SIAM J. Sci. Comput. 16(5), 1190–1208 (1995). [CrossRef]  

28. A. Kak and M. Slaney, “Principles of computerized tomographic imaging IEEE press,” New York (1988).

29. C. K. Egan, S. D. Jacques, T. Connolley, M. D. Wilson, M. C. Veale, P. Seller, and R. J. Cernik, “Dark-field hyperspectral X-ray imaging,” Proc. R. Soc. A 470(2165), 20130629 (2014). [CrossRef]  

30. D. Pennicard, B. Pirard, O. Tolbanov, and K. Iniewski, “Semiconductor materials for X-ray detectors,” MRS Bull. 42(06), 445–450 (2017). [CrossRef]  

31. “Hamamatsu handbook of X-ray detectors,” https://www.hamamatsu.com/resources/pdf/ssd/e09_handbook_xray_detectors.pdf.

32. J. Webber, E. T. Quinto, and E. L. Miller, “A joint reconstruction and lambda tomography regularization technique for energy-resolved X-ray imaging,” Inverse Problems (2020).

33. H. Andrade-Loarca, G. Kutyniok, and O. Öktem, “Shearlets as feature extractor for semantic edge detection: The model-based and data-driven realm,” arXiv preprint arXiv:1911.12159 (2019).

34. A. A. Taha and A. Hanbury, “Metrics for evaluating 3D medical image segmentation: analysis, selection, and tool,” BMC Med. Imaging 15(1), 29 (2015). [CrossRef]  

35. J. Canny, “A computational approach to edge detection,” IEEE Trans. Pattern Anal. Mach. Intell. PAMI-8(6), 679–698 (1986). [CrossRef]  

36. A. Bradley and A. Jay, “A method for deducing accurate values of the lattice spacing from X-ray powder photographs taken by the Debye-Scherrer method,” Proc. Phys. Soc. 44(5), 563–579 (1932). [CrossRef]  

37. A. Taylor and H. Sinclair, “On the determination of lattice parameters by the Debye-Scherrer method,” Proc. Phys. Soc. 57(2), 126–135 (1945). [CrossRef]  

38. C. Zhu, C. Arson, and A. Pouya, “Theoretical study of damage accommodation in salt subject to viscous fatigue,” in Mechanical Behaviour of Salt VIII, (CRC Press, 2015), pp. 331–342.

39. R. Bryan, “International tables for crystallography. Vol. C. Mathematical, physical and chemical tables edited by A. J. C. Wilson,” (1993).

40. J. J. DeMarco and P. Suortti, “Effect of scattering on the attenuation of X-rays,” Phys. Rev. B 4(4), 1028–1033 (1971). [CrossRef]  

41. J. H. Hubbell, W. J. Veigele, E. Briggs, R. Brown, D. Cromer, and d. R. Howerton, “Atomic form factors, incoherent scattering functions, and photon scattering cross sections,” J. Phys. Chem. Ref. Data 4(3), 471–538 (1975). [CrossRef]  

42. J. H. Hubbell and I. Overbo, “Relativistic atomic form factors and photon coherent scattering cross sections,” J. Phys. Chem. Ref. Data 8(1), 69–106 (1979). [CrossRef]  

43. G. G. Poludniowski and P. M. Evans, “Calculation of X-ray spectra emerging from an X-ray tube. part I. electron penetration characteristics in X-ray targets,” Med. Phys. 34(6Part1), 2164–2174 (2007). [CrossRef]  

44. G. G. Poludniowski, “Calculation of X-ray spectra emerging from an X-ray tube. part II. X-ray production and filtration in X-ray targets,” Med. Phys. 34(6Part1), 2175–2186 (2007). [CrossRef]  

45. G. Poludniowski, G. Landry, F. Deblois, P. Evans, and F. Verhaegen, “Spekcalc: a program to calculate photon spectra from tungsten anode X-ray tubes,” Phys. Med. Biol. 54(19), N433–N438 (2009). [CrossRef]  

46. G. Collaboration and S. Agostinelli, “GEANT4–a simulation toolkit,” Nucl. Instrum. Meth. A 506(3), 250–303 (2003). [CrossRef]  

47. G. Harding, J. Kosanetzky, and U. Neitzel, “X-ray diffraction computed tomography,” Med. Phys. 14(4), 515–525 (1987). [CrossRef]  

48. O. Semerci, N. Hao, M. E. Kilmer, and E. L. Miller, “Tensor-based formulation and nuclear norm regularization for multienergy computed tomography,” IEEE Trans. on Image Process. 23(4), 1678–1693 (2014). [CrossRef]  

49. V. Palamodov, “An analytic reconstruction for the Compton scattering tomography in a plane,” Inverse Probl. 27(12), 125004 (2011). [CrossRef]  

50. M. Nguyen and T. Truong, “Inversion of a new circular-arc Radon transform for Compton scattering tomography,” Inverse Probl. 26(6), 065005 (2010). [CrossRef]  

51. G. Rigaud, M. K. Nguyen, and A. K. Louis, “Novel numerical inversions of two circular-arc Radon transforms in Compton scattering tomography,” Inverse Probl. Sci. Eng. 20(6), 809–839 (2012). [CrossRef]  

52. T. T. Truong and M. K. Nguyen, “New properties of the V-line Radon transform and their imaging applications,” J. Phys. A: Math. Theor. 48(40), 405204 (2015). [CrossRef]  

53. T. T. Truong, M. K. Nguyen, and H. Zaidi, “The mathematical foundations of 3D Compton scatter emission imaging,” Int. J. Biomed. Imaging 2007, 1–11 (2007). [CrossRef]  

54. G. Rigaud and B. N. Hahn, “3D Compton scattering imaging and contour reconstruction for a class of Radon transforms,” Inverse Probl. 34(7), 075004 (2018). [CrossRef]  

55. C.-Y. Jung and S. Moon, “Inversion formulas for cone transforms arising in application of Compton cameras,” Inverse Probl. 31(1), 015006 (2015). [CrossRef]  

56. O. Klein and Y. Nishina, “The scattering of light by free electrons according to Dirac’s new relativistic dynamics,” Nature 122(3072), 398–399 (1928). [CrossRef]  

57. O. Klein and Y. Nishina, “Über die Streuung von Strahlung durch freie Elektronen nach der neuen relativistischen Quantendynamik von Dirac,” Z. Physik 52(11-12), 853–868 (1929). [CrossRef]