# Material Identification From Radiographs
Without Energy Resolution

Michael T. McCann, , Elena Guardincerri, Samuel M. Gonzales,
Lauren A. Misurek, Jennifer L. Schei, and Marc L. Klasky All authors are with
Los Alamos National Laboratory, Los Alamos, NM, 87545 USA. M. T. McCann and M.
L. Klasky and are with the Theoretical Division; S. M Gonzales and L. A.
Misurek are with the Nuclear Engineering and Nonproliferation Division; E.
Guardincerri is with the Physics Division; and J. L. Schei is with the
Integrated Weapons Experiments Division.

###### Abstract

We propose a method for performing material identification from radiographs
without energy-resolved measurements. Material identification has a wide
variety of applications, including in biomedical imaging, nondestructive
testing, and security. While existing techniques for radiographic material
identification make use of dual energy sources, energy-resolving detectors, or
additional (e.g., neutron) measurements, such setups are not always
practical—requiring additional hardware and complicating imaging. We tackle
material identification without energy resolution, allowing standard X-ray
systems to provide material identification information without requiring
additional hardware. Assuming a setting where the geometry of each object in
the scene is known and the materials come from a known set of possible
materials, we pose the problem as a combinatorial optimization with a loss
function that accounts for the presence of scatter and an unknown gain and
propose a branch and bound algorithm to efficiently solve it. We present
experiments on both synthetic data and real, experimental data with relevance
to security applications—thick, dense objects imaged with MeV X-rays. We show
that material identification can be efficient and accurate, for example, in a
scene with three shells (two copper, one aluminum), our algorithm ran in six
minutes on a consumer-level laptop and identified the correct materials as
being among the top 10 best matches out of 8,000 possibilities.

###### Index Terms:

X-ray radiography, material identification, combinatorial optimization, branch
and bound algorithm.

##  I Introduction and Related Work

X-ray radiography is an important tool in many security applications,
including airport security [1, 2, 3, 4, 5, 6], cargo inspection [7, 8, 9, 10],
and in the emergency response to a suspicious item that may contain explosives
or other hazards, i.e., for bomb squads [11, 12]. Automating aspects of X-ray-
based screening is attractive in terms of speed, cost, and safety and is an
active area of study; see [13] for a recent review. X-ray radiography provides
both information about the internal geometry of what is being imaged and about
its material composition, because different materials absorb X-rays to
different degrees. Both types of information can contribute to the security
application.

Here, we are interested in using X-ray radiography for material
identification. There are several related yet distinct problems that may be
classified as radiographic material identification. In one setting, the
possible materials come from a finite, discrete set, e.g., [8] aims to
identify whether a slab of material is plywood, polyethylene, aluminum, iron,
copper, or lead. In another setting, the goal is to determine a continuous-
valued material property for each object (or each pixel in the radiograph)
that can be used to infer what materials are present, e.g. [14] estimates
atomic number (as a continuous quantity) in order to distinguish bone,
aluminum, soft tissue, and acrylic. Along similar lines, some approaches
attempt to determine the contribution of multiple materials to each
pixel/voxel of the reconstruction; this is sometimes called _material
decomposition_ [15]. Some works, e.g., [16], determine only _whether_ each
material is present, while others, e.g. [4], determine _where_ each material
is, i.e., segment the radiograph.

Many approaches to material identification from X-ray measurements exist in
the literature. One class of approaches relies on hardware such as dual energy
scanners [17, 18, 19, 20, 21] or photon-counting detectors [22, 23] to collect
energy-resolved measurements; these methods date back to at least the 1970s
[24]. In these cases, material identification can be achieved by solving a
system of equations that relate the energy-dependent measurements to material
properties, e.g., Compton and photoelectric coefficients [25]. Work exists
[26] that aims to elucidate the fundamental limits of dual-energy approaches.
Another, more recent approach, has been to use phase contrast X-ray imaging to
decompose materials [27, 28]. Finally, another class of approaches combines
X-ray measurements with another modality, e.g., neutron imaging [29, 8, 30].
Here, quantities such as the ratio of neutron to photo absorption can be
indicative of material identity.

In this work, we focus on the problem of identifying the materials present in
a radiograph without energy-resolved measurements or additional imaging
modalities, thus allowing a standard MeV X-ray system to provide material
identification information without requiring additional hardware. Solving this
problem could provide an important new capability for security applications,
but, to our knowledge, it has not been addressed in the literature.
Specifically, it could reduce the cost of security imaging systems (because
the source would only need to generate a single X-ray spectrum and the
detector would not need to provide energy resolution). It could also
accelerate and simplify the imaging procedure (because it requires only a
single radiograph) which is of particular importance in a suspicious package
scenario where radiography must be performed quickly, accurately, and
potentially in uncontrolled settings, i.e., “in the field”. Finally, it could
help to reduce the incidental radiation exposure received by people near a
threat device, e.g., in the unfortunate scenario described in [11].

Material identification without energy resolution is a challenging,
underdetermined problem which is impossible to solve in general (e.g., because
different materials of different thicknesses can absorb exactly the same
amount of X-ray radiation); in order to progress, we make a number of
simplifying assumptions that are nonetheless well-justified for security
applications. First, we treat each voxel as containing only one material from
a known list, rather than a mixture; this is reasonable because we expect to
image man-made objects and domain knowledge can inform what materials to look
for. Second, we assume that the geometry of the scene is already known; this
is also reasonable because the geometry can be determined e.g., by inference
from the measured radiograph based on domain knowledge, previous interrogation
with soft X-rays, conventional photography, or additional physical
measurements. Finally, our experiments use MeV X-rays and thick, highly
absorbing objects (see Figure 1), again this is justified because, in security
applications, we are looking for man-made, usually metallic objects.

Our specific contributions are

  1. 1.

the formulation of a novel and important material identification problem;

  2. 2.

an algorithm to efficiently solve the resulting combinatorial problem,
allowing scenes with up to eight objects to be processed in minutes on a
laptop computer; and

  3. 3.

experimental validation of the algorithm on real and synthetic data.

The remainder of the manuscript is organized as follows. In Section II we
formulate the material identification problem and present our proposed
algorithm for solving it. In Section III, we describe validation experiments
on real and synthetic data. We discuss these results in Section IV and
conclude in Section VI.

##  II Proposed Approach

We now formulate the material identification problem and present our
radiograph forward model, branch and bound algorithm, and preprocessing and
calibration procedures.

###  II-A Problem Formulation

Given a properly preprocessed radiograph ${\bm{t}}\in\mathbb{R}^{V\times U}$,
where $V$ and $U$ denote the height and width in pixels, of a scene with
${N_{\text{obj}}}$ single-material objects111Extending our approach to
multiple radiographs is straightforward; we present the single radiograph
version here for brevity and because our experiments use single radiographs.
and a set of possible materials ${\Omega}$ with
$|{\Omega}|={N_{\text{mats}}}$, our goal is to determine the material of each
object: ${\bm{x}}\in{\Omega}^{N_{\text{obj}}}$. For example radiographs, see
Figure 1. Note by our definition an object is a connected region in space
consisting of a single material: an empty coffee cup is one object; a full
coffee cup is two: mug and coffee. We assume that we have the following
ingredients of a radiographic forward model: a raytrace for each object at
each pixel, $\\{\bm{\ell}_{n}\\}_{n=1}^{N_{\text{obj}}}$ with
$\bm{\ell}_{n}\in\mathbb{R}^{V\times U}$, the attenuation coefficients for
each material at each of ${N_{\text{e}}}$ energies
$\\{{\bm{\mu}}_{m}\\}_{m\in{\Omega}}$ with ${\bm{\mu}}_{m}\in\mathbb{R}_{\geq
0}^{N_{\text{e}}}$, the spectrum of the source
${\bm{I}_{s}}\in\mathbb{R}_{\geq 0}^{N_{\text{e}}}$, and the (energy-
dependent) detector response ${\bm{g}}\in\mathbb{R}_{\geq 0}^{N_{\text{e}}}$.

5 cm (a) Al-Cu Shells. Left: Cu stepwedge; outer shell: Al; inner shell: Cu;
right: Al stepwedge.

5 cm (b) Eight Cubes. Top row: Cu, Bi, Fe, Ta; bottom row: Ti, W, Mo, Al;
center: resolution target; right: Al stepwedge.

5 cm (c) Al-Cu Shells Brem. Outer shell: Al; inner shell: Cu; right: Cu
stepwedge.

5 cm (d) Al-Cu-Cu Shells. Outer shell: Al; middle shell: Cu; inner shell: Cu;
right: Cu stepwedge.

5 cm (e) Eight Cubes Brem. Top row: Mo, W, Cu, Fe; bottom row: Al, Ta, Bi, Ti;
right: Cu step wedge.

Figure 1: Radiographs (after preprocessing, see Section II-D) of the scenes
used for material identification. Images are $2352\times 2878$.

Practically speaking, it is useful for algorithms to return a ranked list of
guesses for ${\bm{x}}$, thereby reducing the potentially large number of
possible material assignments to a tractable number. For example, for a bomb
squad investigating a suspicious package, this would allow the team to focus
their efforts on a handful of possible weapon designs. Thus, for performance
evaluation, we consider top-$N$ accuracy, i.e., an algorithm is top-$N$
correct if the true material assignment ${\bm{x}_{\text{gt}}}$ is among the
$N$ ${\bm{x}}$s returned. We pose the material identification problem
described above in terms of minimizing the functional

$J({\bm{x}})=\min_{\alpha,\theta}\|{\bm{t}}-(\alpha{\bm{d}}({\bm{x}})+{\bm{s}}(\theta))\|_{2}^{2},$
(1)

where ${\bm{t}}$ is the measured radiograph, $\alpha$ is a scalar, ${\bm{d}}$
is our model of the direct signal, and ${\bm{s}}(\theta)$ is our parametric
scatter model with parameters $\theta$ (which we detail later in (4)). Seeking
a set of top-$N$ material assignments amounts to finding

$\begin{split}{\bm{x}}^{*}_{1}&\in\operatorname*{argmin}_{{\bm{x}}\in{\Omega}^{N_{\text{obj}}}}\
J({\bm{x}})\\\
{\bm{x}}^{*}_{2}&\in\operatorname*{argmin}_{{\bm{x}}\in{\Omega}^{N_{\text{obj}}}\setminus{\bm{x}}^{*}_{1}}J({\bm{x}})\\\
&\vdots\\\
{\bm{x}}^{*}_{N}&\in\operatorname*{argmin}_{{\bm{x}}\in{\Omega}^{N_{\text{obj}}}\setminus{\bm{x}}^{*}_{1},{\bm{x}}^{*}_{2},\dots,{\bm{x}}^{*}_{N-1}}J({\bm{x}}),\end{split}$
(2)

where $\setminus$ denotes set difference. We note that, in general, the
problem (2) is underdetermined in that several different material assignments
may explain the measured radiograph equally well; we expect such cases to be
rare in practice and, even if they do occur, it can still be valuable to find
these assignments and rule out any others.

In the following sections, we describe each element of the functional
$\eqref{eq:functional}$ and how we solve the minimization problem (2).

###  II-B Radiographic Forward Model

We model the signal entering the detector as the sum of the direct and
scattered signal. The direct signal accounts for photons that travel directly
from the source to the detector; objects create contrast in the direct signal
because they remove photons from that straight path (either by absorption or
scattering) according to Beer’s law. The scattered signal represents photons
that scatter from the objects of interest as well as the remainder of the
scene (ground, walls, detector, etc.); this signal is difficult to model
accurately. We now describe each in detail.

We model the direct signal entering the detector at pixel ${\bm{r}}$ for a
scene with materials ${\bm{x}}$ as

${\bm{d}}({\bm{x}})[{\bm{r}}]=\sum_{e}{\bm{g}}[e]{\bm{I}_{s}}[e]\exp\left(-\sum_{n=1}^{{N_{\text{obj}}}}{\bm{\mu}}_{{\bm{x}}[n]}[e]\bm{\ell}_{n}[{\bm{r}}]\right),$
(3)

where ${\bm{g}}$, ${\bm{I}_{s}}$, $\\{{\bm{\mu}}_{m}\\}_{m\in{\Omega}}$, and
$\\{\bm{\ell}_{n}\\}_{n=1}^{N_{\text{obj}}}$ are known physical constants and
aspects of the scene as described in II-A. This model is a standard
discretized version of the polyenergetic Beer’s law [31]; for more detail on
the discretization, see [32]. The direct signal model (3) has no term to
account for spatial variations in the intensity of the source, this is handled
during preprocessing of the radiograph, see Section II-D

We do not directly model scatter; instead, our cost functional
$\eqref{eq:functional}$ is designed so that we fit a parametric function,
${\bm{s}}$ to the difference between the prediction of the forward model
${\bm{d}}({\bm{x}})$ and the measured radiograph ${\bm{t}}$. The effect of
this fitting is that when the difference between ${\bm{d}}({\bm{x}})$ and
${\bm{t}}$ is in the (nonlinear) span of ${\bm{s}}$, the discrepancy between
them is zero, i.e., $J({\bm{x}})=0$. We emphasize that, unlike other scatter
modeling approaches [33, 34], our scatter model need not be predictive.
Instead, it needs to strike a balance between representing enough—so that it
can bridge the gap between our model for the direct signal and what is
actually measured—and too much—so that it does not make incorrect material
assignments appear to match the measured radiograph.

We use low-order polynomial fields to model the scatter, ${\bm{s}}$. These
allow us to account for smoothly-varying scatter, enable rapid fitting during
the evaluation of the cost function (1), and provide control over their
complexity via changing the polynomial order. Specifically,

${\bm{s}}(\theta)[{\bm{r}}]=\sum_{m\geq 0,n\geq 0,m+n\leq
P}\theta_{m,n}u({\bm{r}})^{m}v({\bm{r}})^{n},$ (4)

where $P$ is the degree of the polynomial, $\theta_{m,n}\in\mathbb{R}$ is one
element of the scatter model parameters, $\theta$, and $u({\bm{r}})$ and
$v({\bm{r}})$ give the $u$ and $v$ coordinate of the pixel ${\bm{r}}$ (where
the coordinate system is arbitrary; we place $(0,0)$ at the upper-left corner
of the radiograph with positive $u$ to the right and positive $v$ down). For a
demonstration that polynomial fields are a reasonable way to approximate
scatter, see Section III-E.

###  II-C Branch and Bound Algorithm

The optimization problem (2) is a combinatorial problem: the space of
solutions is discrete and finite, so the problem may be solved via exhaustive
search, however the size of the space grows exponentially with the number of
objects (${N_{\text{mats}}}^{N_{\text{obj}}}$). For example, performing
material identification among twenty materials on five objects requires
evaluating the forward model $20^{5}\approx 3.2$ million times. For
comparison, solving a convex problem would typically involve evaluating the
forward model and its adjoint only hundreds of times. While this search may be
parallelized, with even a few more materials it quickly becomes infeasible.

In order to efficiently solve the optimization problem (2), we propose an
algorithm in the branch and bound family. Branch and bound algorithms rely on
finding lower bounds on the objective value of partial solutions to a
combinatorial problem; see [35], Chapter 5 for an introduction. To design a
branch and bound algorithm for our problem, we need a way to express partial
solutions and bound their objective values. We also need to extend the
standard branch and bound framework to return the top $N$ solutions rather
than simply the best. (Such extensions have been used before in other
contexts, e.g., [36, 37].)

Algorithm 1 Branch and bound for top-$N$ optimization

1:$topN\leftarrow\\{\\}$, $J^{*}\leftarrow\infty$,
$cands\leftarrow\\{[\varnothing,\dots,\varnothing]\\}$

2:while $cands$ is not empty do

3: Remove a candidate ${\bm{x}}$ from $cands$

4: if ${\bm{x}}$ is a full solution and $J({\bm{x}})<J^{*}$ then

5: Add ${\bm{x}}$ to $topN$

6: if $|topN|>N$, remove $\operatorname*{argmax}_{{\bm{x}}\in
topN}J({\bm{x}})$

7: $J^{*}\leftarrow\max_{{\bm{x}}\in topN}J({\bm{x}})$

8: if ${\bm{x}}$ is a partial solution and $bound({\bm{x}})<J^{*}$ then

9: $cands\leftarrow cands\cup branch({\bm{x}})$

We outline our branch and bound algorithm in Algorithm 1. The approach follows
a standard branch and bound structure with the exception that a list of
solutions, $topN$, must be maintained instead of a single best solution. The
$branch$ and $bound$ subroutines are novel and adapted to the problem at hand.
Our bounding approach makes use of the fact that, for typical radiographs,
only a subset of objects are present at each pixel; thereby allowing us to
compute a lower bound on the error at that pixel for any partial solution that
assigns materials to each of the objects that are present. We now describe the
algorithm in detail.

To represent partial solutions, we augment our list of possible materials with
an additional symbol, $\varnothing$, which we take to mean “not yet assigned”;
partial solutions are then elements of
$({\Omega}\cup\varnothing)^{N_{\text{obj}}}$ and full solutions are the
elements of the same set that do not involve any $\varnothing$s. We will say
that a partial or full solution $\bm{y}$ is a descendent of a partial solution
${\bm{x}}$, written $\bm{y}\in{\bm{x}}$, when $\bm{y}$ matches ${\bm{x}}$ in
every non-$\varnothing$ position of ${\bm{x}}$. The list of candidate
solutions, $cands$, is initialized with a single partial solution that is full
of $\varnothing$s, this represents all possible material assignments.

The subroutine $branch({\bm{x}})$ creates new partial (or possibly full)
solutions from a partial solution. To do this, we replace the first
$\varnothing$ in ${\bm{x}}$ with each material in the set of materials
${\Omega}$. For example, if ${\Omega}=\\{\text{Cu},\text{Ni},\text{Sn}\\}$,
then
$branch([\text{Cu},\varnothing,\varnothing])=\\{[\text{Cu},\text{Cu},\varnothing],[\text{Cu},\text{Ni},\varnothing],[\text{Cu},\text{Sn},\varnothing]\\}$.

The subroutine $bound({\bm{x}})$ produces a lower bound on the value of
$J(\bm{y})$ for any full solution $\bm{y}\in{\bm{x}}$. Good bounds are
critical to the efficiency of the algorithm because they are what allow bad
partial solutions (and all of their descendant full solutions) to be skipped.
To design $bound$, we make use of the fact that, in many real-world scenarios,
only a subset of materials are present at each pixel in the radiograph.
Because of this, when a partial solution fixes all the materials present at a
set of pixels, the error incurred at those pixels cannot be reduced by making
further assignments and therefore provides a lower bound on the total error
for the partial solution. Specifically, to compute $bound({\bm{x}})$ we note
that

$\begin{split}\min_{\bm{y}\in{\bm{x}}}\min_{\alpha,\theta}\|{\bm{t}}-(\alpha{\bm{d}}(\bm{y})&+{\bm{s}}(\theta))\|_{2}^{2}\\\
&\geq\min_{\bm{y}\in{\bm{x}}}\min_{\alpha,\theta}\|\bm{M}_{\bm{x}}({\bm{t}}-(\alpha{\bm{d}}(\bm{y})+{\bm{s}}(\theta)))\|_{2}^{2}\\\
&=\min_{\alpha,\theta}\|\bm{M}_{\bm{x}}({\bm{t}}-(\alpha{\bm{d}}(\bm{y}_{0})+{\bm{s}}(\theta)))\|_{2}^{2},\end{split}$
(5)

where $\bm{M}_{\bm{x}}$ is a binary masking operator chosen to mask away any
pixels that depend on $\varnothing$ values in ${\bm{x}}$ and $\bm{y}_{0}$ is
any full solution descendent of ${\bm{x}}$. The inequality holds because the
fact that $\|\cdot\|_{2}^{2}$ is a sum of squares means that, for all
($\bm{y},\alpha,\theta$), we have
$\|{\bm{t}}-(\alpha{\bm{d}}(\bm{y})+{\bm{s}}(\theta))\|_{2}^{2}\geq\|\bm{M}_{\bm{x}}({\bm{t}}-(\alpha{\bm{d}}(\bm{y})+{\bm{s}}(\theta)))\|_{2}^{2}$;
thus the same relationship holds for the minimum values. The equality holds
because the mask $\bm{M}_{\bm{x}}$ is defined to zero out exactly the pixels
that would make two different full solution descents of $\bm{x}$ have
different functional values. The mask is recomputed for each partial solution
${\bm{x}}$, but this operation is fast because it involves only pixelwise,
boolean operations on binary masks.

The final piece of our branch and bound algorithm is to determine the order in
which candidate solutions are evaluated, i.e., the method for selecting a
candidate from $cands$ in line 3 of Algorithm 1. Common options are depth-
first (first in, last out), breadth-first (first in, first out), or best-first
(least bound, first out). In our context, we need to avoid breadth-first
search because the $cands$ list would become huge (containing nearly every
possible partial solution) before $J^{*}$ was ever updated. Depth-first search
has a different shortcoming, in that it can spend a lot of time exploring
solutions with poor initial material assignments. Finally, best-first search
ends up being similar to breadth-first search because we tend to give overly
optimistic bounds to partial solutions with few assignments. Our solution is
use first in, last out ordering (depth-first search) but to also sort new
candidates according to their $bound$ value when adding them to $cands$. The
result is a depth-first search that chooses the most promising candidate at
each depth first.

###  II-D Preprocessing

We have so far assumed that we have access to a properly preprocessed
radiograph, ${\bm{t}}$. Here, we detail our preprocessing steps, which include
detector pixel calibration and flat-field correction.

To account for the fact that each detector pixel may have its own gain and
offset value (also called _pedestal_), we model a (uncorrected) radiograph,
${\bm{w}}$ as

${\bm{w}}[{\bm{r}}]=t\alpha[{\bm{r}}]{\bm{t}}[{\bm{r}}]+t\beta[{\bm{r}}]+\gamma[{\bm{r}}]$
(6)

where $\alpha$, $\beta$, and $\gamma$ represent pixelwise gain, dark count
rate, and offset values, ${\bm{t}}$ is the unnormalized total transmission,
and $t$ is the exposure time. The goal of our preprocessing is to use
calibration data to estimate $\alpha$, $\beta$, and $\gamma$ in order to undo
their effects and to simultaneously remove the effect of spatial variation in
the X-ray source intensity.

As calibration data, we use a set of flat-field radiographs (radiographs of an
empty scene) taken using a high activity ${}^{60}Co$ radioactive source and a
set of dark-field radiographs (radiographs without an X-ray source). The
reason for using a radioactive source as opposed to a betatron is that for the
former the dose rate at the detector is known and constant in time, whereas
for the latter this is not true. We can assume that there is no transmission
present in the dark-field radiographs, so ${\bm{t}}_{\text{dark}}[r]=0$. For
the flat-field radiographs, we can write

${\bm{t}}_{\text{flat}}[{\bm{r}}]={\bm{I}_{p}}[{\bm{r}}]\sum_{e}{\bm{g}}[e]{\bm{I}_{s}}[e]+{\bm{s}}_{\text{flat}}[{\bm{r}}]$
(7)

where ${\bm{I}_{p}}$ is the spatial source profile, ${\bm{g}}$ is the energy-
dependent detector response, ${\bm{I}_{s}}$ is the beam spectrum, and
${\bm{s}}_{\text{flat}}$ is an unknown scatter term. Note that (7) assumes
that the beam spectrum and detector response are multiplicatively separable
with respect to position, ${\bm{r}}$, and energy, $e$; if the relative
contribution of different energies change significantly between pixels, it
would not be accounted for by this calibration. With the approximation that
${\bm{s}}_{\text{flat}}[{\bm{r}}]=0$ and the convention to scale ${\bm{g}}$
and ${\bm{I}_{s}}$ so that their inner product is one, we have
${\bm{t}}_{\text{flat}}={\bm{I}_{p}}$. Substituting into (6), we have

$\displaystyle{\bm{w}}_{\text{dark}}[{\bm{r}}]$
$\displaystyle=t\beta[{\bm{r}}]+\gamma[{\bm{r}}]$ (8)
$\displaystyle{\bm{w}}_{\text{flat}}[{\bm{r}}]$
$\displaystyle=t\alpha[{\bm{r}}]{\bm{I}_{p}}[{\bm{r}}]+t\beta[{\bm{r}}]+\gamma[{\bm{r}}].$
(9)

Therefore, using several dark- and flat-field radiographs, we can estimate
$\alpha[{\bm{r}}]{\bm{I}_{p}}[{\bm{r}}]$, $\beta[{\bm{r}}]$, and
$\gamma[{\bm{r}}]$ by solving the linear system

$\begin{bmatrix}{\bm{w}}_{\text{dark},1}[{\bm{r}}]\\\
{\bm{w}}_{\text{dark},2}[{\bm{r}}]\\\ \vdots\\\
{\bm{w}}_{\text{dark},N_{\text{dark}}}[{\bm{r}}]\\\
{\bm{w}}_{\text{flat},1}[{\bm{r}}]\\\ {\bm{w}}_{\text{flat},2}[{\bm{r}}]\\\
\vdots\\\ {\bm{w}}_{\text{flat},N_{\text{flat}}}[{\bm{r}}]\\\
\end{bmatrix}=\begin{bmatrix}0&t_{\text{dark},1}&1\\\ 0&t_{\text{dark},2}&1\\\
\vdots&\vdots&\vdots\\\ 0&t_{\text{dark},N_{\text{dark}}}&1\\\
t_{\text{flat},1}&t_{\text{flat},1}&1\\\
t_{\text{flat},2}&t_{\text{flat},2}&1\\\ \vdots&\vdots&\vdots\\\
t_{\text{flat},N_{\text{flat}}}&t_{\text{flat},N_{\text{flat}}}&1\end{bmatrix}\begin{bmatrix}\alpha[{\bm{r}}]{\bm{I}_{p}}[{\bm{r}}]\\\
\beta[{\bm{r}}]\\\ \gamma[{\bm{r}}]\end{bmatrix},$ (10)

To account for stuck or dead pixels on the detector, we set a threshold on the
R-squared value of the fit and remove each pixel that falls below it.

Finally, to preprocess a radiograph with exposure time $t$, we use

$\text{prep}({\bm{w}})[{\bm{r}}]=\frac{{\bm{w}}[{\bm{r}}]-t\beta[{\bm{r}}]-\gamma[{\bm{r}}]}{t\alpha[{\bm{r}}]{\bm{I}_{p}}[{\bm{r}}]}.$
(11)

Substituting (11) into (6) confirms that

$\text{prep}({\bm{w}})[{\bm{r}}]=\frac{{\bm{t}}[{\bm{r}}]}{{\bm{I}_{p}}[{\bm{r}}]}.$
(12)

In short, our preprocessing removes detector gain and pedestal and also
normalizes the transmission with respect to the flat-field image.

###  II-E Source Spectrum and Detector Response Calibration

Our model depends on knowing the product of the source spectrum and the
detector response over the energies where they are nonzero. In practice, there
are challenges to knowing each of these: The spectrum for a line source is
known but affected by shielding around the source, e.g., the container holding
it. A betatron generates, ideally, a Bremsstrahlung spectrum, but in practice
involves random fluctuations as well as systematic changes as the source
warms. Detector responses may be provided by the manufacturer; if they are
not, they may be estimated from the materials used in the detector, but both
of these responses come with some degree of uncertainty and are difficult to
validate.

To overcome these challenges, we propose to calibrate the quantity
${\bm{I}_{s}}{\bm{g}}$ from a calibration radiograph by solving

$\operatorname*{argmin}_{{\bm{I}_{s}}{\bm{g}}}J({\bm{x}_{\text{gt}}}),$ (13)

where the data fidelity term $J$ is given in (1) and ${\bm{x}_{\text{gt}}}$ is
the ground truth material assignment. Note that computing
$J({\bm{x}_{\text{gt}}})$ itself involves computing a linear fit and that the
problem is generally ill-posed and nonconvex. To minimize (13) we used
standard gradient descent, where $\nabla J$ was computed using automatic
differentiation software [38].

##  III Experiments and Results

We now describe our experiments on real and synthetic data. All material
identifications were performed on a consumer-level laptop (MacBook Pro, 2.4
GHz processor, 32 GB of RAM). Objects were identified as being one of the 20
materials listed in Table I. Attenuation coefficients for these materials were
obtained from the MCNP particle transport software [39], specifically the
mcplib04 library, which derives from ENDF/B-VI.8.

TABLE I: Materials list air | aluminum | beryllium | bismuth
---|---|---|---
boron | carbon | copper | gold
high explosive | iron | lead | lithium
molybdenum | plutonium | polycarbonate | tantalum
tin | titanium | tungsten | uranium

###  III-A Material Identification From Experimental Data—Cobalt-60 Source

As a first validation with experimental data, we performed material
identification on radiographs taken with a cobalt-60 source. The source had an
activity of 3.501 Ci at the time of the experiment. All experiments were
performed in the basement of at the Nonproliferation International Security
Center at the Los Alamos National Laboratory; all the images were collected
using a Pixium 3543EZ digital readout panel. For the source spectrum, we used
the two emission lines of cobalt-60: 1.17 MeV and 1.33 MeV. We assumed the
detector response was the same for both lines and did not perform further
calibration. The source-to-detector distance was 100 cm and the object-to-
detector distance was 17.5 cm.

Al-Cu Shells. The first scene consisted of concentric shells of aluminum (1
inch thick, 5 inch outer diameter) and copper (0.5 inch thick, 4 inch outer
diameter) and aluminum and copper step wedges. The scene was imaged with an
exposure time of 100 seconds; the radiograph is shown in Figure 1. The task
was to identify the materials of the two shells from the list of possible
materials in Table I. We preprocessed the image (Section II-D) and cropped out
the table, step wedges, and stand.

In order to compute path lengths through the shells, we assumed that the
center of the sphere was on the optical axis and used the known outer diameter
of the aluminum shell, source-to-shells distance, source-to-detector distance,
and diameter of the sphere in pixels to fix the rest of the scene geometry. We
then raytraced the shells analytically. The result of this process was an
image of the path length at each pixel for each shell, $\bm{\ell}_{0}$ and
$\bm{\ell}_{1}$; Figure 2 shows $\bm{\ell}_{0}$.

05001000150020000500100015002000250001234567path length (cm) Figure 2: Path
length of the aluminum shell in the Al-Cu Shells scene, $\bm{\ell}_{0}$.

We ran the branch and bound algorithm (Section II-C) with scatter model order
$P=2$ to identify the top 20 material configurations for the scene. The
correct material assignment, [aluminum, copper], was the one with the lowest
fitting error (RMSE=0.0128), i.e., it was the value ${\bm{x}}$ that minimized
the loss functional (1) out of all possible assignments. meaning that our
algorithm was able to correcting identify these materials; the full ranking
and RMSEs for the top 20 assignments is given in Figure 3.

The total number of possible configurations was $20^{2}=400$ and the branch
and bound algorithm evaluated all of them plus an additional 20 partial
solutions, i.e., $[\text{air},\varnothing,]$, $[\text{B},\varnothing,]$,
$[\text{high explosive},\varnothing,]$, …. This process took under three
minutes on a consumer-level laptop.

Figure 3: Results of material identification on the aluminum and copper shell
scene with a cobalt-60 source.

To further explore these results, we plot the direct radiograph that our model
predicted for the material assignment
${\bm{x}}=[\text{aluminum},\text{copper}]$ (Figure 4) and the corresponding
modeled scatter field (Figure 5). The fact that that the modeled direct only
reaches a value around 0.9 in the background region (instead of approximately
1.0, as direct theoretically should) suggests that the gain parameter and the
scatter field parameters ($\alpha$ and $\theta$ in (1)) work together to model
scatter and any other aspect of the real imaging system that our forward model
does not capture. Finally, we plot a lineout across the diameter of the
shells, Figures 6 and 7. These show good agreement with the measured
radiograph, but also show evidence of systematic error in that the predicted
transmission consistently undershoots the measured transmission.

0500100015002000050010001500200025000.10.20.30.40.50.60.70.80.9

Figure 4: Modeled direct radiograph. The isolated white pixels are those that
were identified as dead or stuck pixels during preprocessing (Section II-D).
0500100015002000050010001500200025000.060.070.080.090.100.11 Figure 5: Modeled
scatter radiograph. 600800100012001400160018000.20.40.60.81.0 Figure 6:
Lineout across the diameter of the sphere comparing the measured radiograph to
the first prediction (blue, dashed) and the second prediction (orange,
dotted). 8008509009501000105011000.150.200.250.300.350.400.450.500.55 Figure
7: Zoomed in lineout across the diameter of the sphere comparing the measured
radiograph to the first prediction (blue, dashed) and the second prediction
(orange, dotted).

Eight Cubes. As a more challenging validation, we radiographed a scene of
eight, one-inch cubes made of copper, bismuth, iron, tantalum, titanium,
tungsten, molybdenum, and aluminum (Figure 1). The scene was imaged with an
exposure time of 100 seconds. We set the path length through each cube at one
inch, neglecting the conebeam geometry; this is justified because all the
cubes are within a few inches of the optical axis, which is small compared to
the 82.5 cm source-to-object distance.

For this configuration, we used a warm start technique: we first ran a branch
and bound search with only lithium, tin, and uranium (low, medium, and high
absorbing materials) and used the top 20 assignments as an initialization for
the $topN$ list and the bound $J^{*}$ in Algorithm 1. This approach ensured a
good initial value for the bound, allowing the full search to finish more
quickly.

The initial and full search took approximately 31 minutes on the same
consumer-level laptop. The correct material assignment had an RMSE of
$6.79\times 10^{-3}$ and did not rank among the top twenty assignments, which
had RMSEs between $6.24\times 10^{-3}$ and $6.41\times 10^{-3}$. However all
of the top ten assignments were partially correct (Figure 8) and each cube was
identified corrected in at least one of them.

Figure 8: Material identification accuracy for the top ten assignments in for
the Eight Cubes data. A check mark in row $n$, column $m$ indicates that the
$n$th assignment correctly identified material $m$.

Here, the benefits of the branch and bound algorithm were much more apparent.
There were $20^{8}=2.56\times 10^{10}$ possible configurations, and the full
search used only 130,720 bounding operations and evaluated 12,640 complete
configurations. This means that the vast majority (99.99995%) of the cases
where excluded by the algorithm, representing an approximately $1.8\times
10^{5}$ times speedup over exhaustive search.

###  III-B Material Identification From Experimental Data—Betatron Source

In a separate experiment, we collected radiographs using a SEA-7 betatron
accelerator by Instauro where the energy of the primary electron beam was set
to 7 MeV. We modeled the source spectrum as a 7 MeV Bremsstrahlung spectrum
and used a simulation-based model for the detector response. The response
function of the digital readout panel was calculated using MCNP. The panel was
modeled as a sequence of material layers as indicated by the manufacturer and
shown in table II. Multiple simulations were run to calculate the detector
response. In each of them, a monochromatic X-ray pencil beam was shot on the
panel, perpendicular to it, and the total energy deposited in the sensitive
layers of the panel was obtained. The energy of the incident X-rays was
different in each simulation run and varied in the range (0 MeV, 7 MeV]. The
resulting response function is shown in figure 9.

TABLE II: Pixium 3543EZ digital readout panel layer | thickness (mm) | modeled as
---|---|---
front side housing | 3.6 | plastic
gadox scintillator layer | 0.4 | gadox
photodiode plate | 0.7 | silicon
mechanical assembly | 2.3 | air
lead | 1.0 | lead
electronics | 7.25 | plastic
back side housing | 1.5 | plastic
$\displaystyle{10^{-2}}$$\displaystyle{10^{-1}}$$\displaystyle{10^{0}}$X-ray
energy (MeV)2461e$-$5energy deposited (MeV/g) Figure 9: The response function
of the Pixium 3543EZ digital readout panel calculated using MCNP.
$\displaystyle{10^{-2}}$$\displaystyle{10^{-1}}$$\displaystyle{10^{0}}$X-ray
energy (MeV)0.000.020.040.06normalized responseuncalibratedcalibrated Figure
10: Product of the detector response and source spectrum before and after
calibration.

We subsequently calibrated the product of source spectrum and detector
response on the Al-Cu Shells Brem data using the method described in Section
II-E, see Figure 10 for the results of the calibration. Calibration emphasized
energies between about 0.1 MeV and 1.0 MeV. The source-to-detector distance
was 300 cm and the object-to-detector distance was 100 cm.

Al-Cu Shells Brem. We repeated the Al-Cu Shells scene with the betatron
source, with an exposure time of 8 seconds. We again used the known shell
diameters, source-to-object distance, and source-to-detector distance to
analytically compute the path lengths. Branch and bound again reduced to
exhaustive search and took 7.2 minutes. The correct assignment was at rank 2
(RMSE=0.014); the error for the top 15 results and lineouts for the top 2
results are given in Figure 11.

02004006008001000120014000.30.40.50.60.70.80.91.0

1502002503003504000.250.300.350.400.45

Figure 11: Left: results of material identification for Al-Cu Shells scene
with a betatron source. Center: Lineout across the diameter of the shells
comparing the measured radiograph to the first prediction (blue, dashed) and
the second prediction (orange, dotted). Right: Zoomed lineout.

Al-Cu-Cu Shells Brem. As a more challenging problem, we performed material
identification on a scene with three shells: aluminum, copper, copper (Figure
1). For the outer shells, we analytically computed the path lengths as in the
previous experiments. The center shell in this configuration is not on the
optical axis, which complicates the process of computing its path lengths. We
manually searched for a position (in 3D) for the shell such that when we
computed path lengths for a shell of the correct dimensions, they aligned well
with the measured radiograph. The exposure time was 8 seconds. The search for
the top 20 assignments took under six minutes on the same consumer-level
laptop. The branch and bound algorithm computed error for 140 full assignments
and excluded the remaining 7860 assignments by computing error for 180 partial
assignments; this represents an approximately $25$ times speedup over
exhaustive search. The ground truth was at position six (RMSE=$9.95\times
10^{-3}$); errors for the top 15 assignments and lineouts are shown in Figure
12.

40060080010001200140016000.20.30.40.50.60.70.80.9

8009001000110012000.2000.2250.2500.2750.3000.3250.3500.3750.400

Figure 12: Left: results of material identification for the Al Cu Cu Shells
with a betatron source. Center: Lineout across the diameter of the sphere
comparing the measured radiograph to the first prediction (blue, dashed) and
the second prediction (orange, dotted). Right: Zoomed lineout. The
disagreement between measurement and model can be reduced by increasing the
polynomial scatter model order, but at the cost of decreased identification
performance, see Figure 17.

Eight Cubes Brem We repeated the eight cubes experiment with the betatron
source with an exposure time of 8 seconds. We again set the path length
through each cube to one inch, ignoring conebeam effects. Using the same warm
start approach as described for the Eight Cubes data, branch and bound took
5.3 hours. Again, the correct assignment did not appear in the top 20
assignments. The RMSE of the ground truth was $7.91\times 10^{-3}$, with the
top 20 assignments having RMSEs between $5.36\times 10^{-3}$ and $5.77\times
10^{-3}$. The top results still had some correct assignments, but several
materials (molybdenum, tantalum, and bismuth) were not identified correctly in
any of the top 10 assignments, shown in Figure 13.

Figure 13: Material identification accuracy for the top ten assignments in for
the Eight Cubes Brem data. A check mark in row $n$, column $m$ indicates that
the $n$th assignment correctly identified material $m$.

###  III-C Effect of Material Thickness

To evaluate the effect of material thickness on material identification
performance, we made a series of simulated radiographs with the MCNP Monte
Carlo particle transport code [39]. For all simulations, we used a 6 MeV
Bremsstrahlung source with an source-to-object distance of 4 m and a source-
to-detector distance of 6.5 m and a detector pixel size of 50 µm. In all
simulations, the scene consisted of a ball of lithium surrounded by shells of
polycarbonate, boron, iron, and plutonium. The simulations differed in that
the radius of the ball and the thickness of the shells varied from 0.1 cm to 2
cm.

Figure 14 shows the results of performing material identification on these
simulations. The ground truth ranked at the first position for the 1 cm
simulation and the second position for the rest, suggesting that the contrast
between materials is largest at intermediate lengths. This makes sense because
with large lengths, transmissions are always small no matter the material, and
vise versa for small lengths. At the smaller lengths, confusions were between
the boron shell and other materials with small attenuation coefficients: high
explosive and carbon.

Figure 14: Rank (bars) and error (crosses) of the ground truth material
assignment as a function of shell thickness, along with errors for the top
five incorrect assignments (points, horizontally spread for visibility) on
simulated data. The 2 cm shells are the most challenging in that more
assignments have an error similar to that of the ground truth.

###  III-D Effect of Noise and Path Length Error

In this section, we examine how noise and error in path length estimation may
impact material identification. In general, this is a complex question: our
algorithm returns a discrete answer, meaning that performance cannot smoothly
degrade as noise or path length errors increase. Further, the effect of noise
and errors will be scene-dependent, limiting the value of any particular
empirical study.

As a partial solution to these challenges, we analytically computed the value
of the forward model (3) using our betatron source spectrum and detector
response for several materials over a range of path lengths, see Figure 15.
This analysis provides an analytical way of determining what level of noise
and path length error is acceptable for a given application: the smallest
vertical cap between a pair of curves puts an upper limit on the acceptable
measurement noise before those two materials will be confused; likewise the
smallest horizontal gap limits the acceptable error in path length estimation.

01234567path length, $\displaystyle\bm{\ell}$ (cm)0.00.20.40.60.81.0direct
signal (normalized)aluminumironcopperlead Figure 15: Expected direct signal
for aluminum, iron, copper, and lead for path lengths up to 7 cm. At 15% path
length error for objects between 1 and 2 cm thick (left gray box) copper may
be confused with iron, because the green line falls within the grey box. At
15% measurement noise for objects between 4 and 5 cm thick (right gray box),
aluminum is unlikely to be confused with any of the other materials.

###  III-E Justification for Polynomial Field Scatter Model

To validate the appropriateness of our polynomial field scatter model (Section
II-B), we performed a photon transport simulation similar to the the Al-Cu
Shells experiment with the cobalt-60 source using MCNP [39]. The resulting
radiograph (including both direct and scattered photons) after flat-field
correction and the scatter radiograph (only scattered photons) are shown in
Figure 16. We observe that the scatter is mostly smooth, although the effect
of the left step-wedge is barely visible.

010020030040050060001002003004005000.20.40.60.81.0

010020030040050060001002003004005000.020.040.060.080.10

01002003004005006007000.000.020.040.060.080.10

Figure 16: Total signal (left), scatter signal (center), and scatter lineouts
with polynomial fits (right) of orders one, two, four, and eight for a
simulated scene. Scatter reaches over 10% of the flat-field signal (note the
different color scales) and can account for more than half of the total signal
in areas where the direct is small.

We fit the scatter field with polynomial fields of order one, two, four, and
eight; Figure 16 shows the resulting lineouts. By order two, the polynomial
fit is reasonably close to the scatter, and by order eight, the maximum
absolute error is less than 0.009 over the whole field of view. In the context
of material identification, the advantage of using lower orders is that it is
more computationally efficient (fewer coefficients lead to faster solutions
during the fitting in (1)) and is less prone to making incorrect material
assignments match the data well.

While this experiment shows that polynomial fields are reasonable for
representing scatter fields, it is not straightforward to generalize these
result to real data. The simulation does not (and can not) represent scene
scatter, i.e. photons leaving the source, interacting with material in the
scene but outside the field of view of the detector, and deflecting back into
the detector. Nor have we accounted for scattering within the detector itself.
It is unlikely we will ever be able to precisely measure the scatter component
of a real radiograph to further validate scatter models. However, we suspect
most of the effects that are difficult to model would tend to make scatter
more spatially smooth than our simulation, again justifying the use of a
smooth field to model scatter.

It is also important that the scatter model _not_ be able to account for too
much of the direct signal. As an extreme case, if we allowed the scatter to be
any image, then scatter would account for the entire measured signal and
material identification using our method would be impossible. This question
boils down to how closely the scatter model (in our case, a low-order
polynomial) can fit a direct signal, which is the path length as a function of
space passed through the polyenergetic Beer’s law (3), which is a pointwise
nonlinearity resembling a decaying exponential. For the types of objects that
we expect to encounter, the path length increases quickly at the edge of the
object, creating a sharp decrease in the direct signal (see, e.g., Figure 6)
that polynomials do not fit well.

###  III-F Effect of Polynomial Scatter Model Order

To explore the effect of the order of the polynomial scatter model, we ran the
Al-Cu-Cu Shells Brem material identification with model orders zero through
ten. The results of this experiment are summarized in Figure 17. The fitting
error of all assignments decreases as the polynomial model order increases,
which must be the case because higher order polynomials fit strictly more of
the residual than lower order ones. What is more important is the error of the
ground truth relative to other, incorrect assignments; this value is at its
minimum for orders two and three, which place the correct assignment at rank
six.

Figure 17: Rank (bars) and error (crosses) of the ground truth material
assignment as a function of polynomial scatter model order, along with errors
for the top ten incorrect assignments (points, horizontally spread for
visibility) on the Al-Cu-Cu Shells Brem data. Model orders two and three are
optimal for this data.

##  IV Discussion

Comparing experiments with the cobalt-60 source to those with the betatron
source, those with the betatron source resulted material identification that
was slower and less accurate. One reason for the difference in speed is that
the forward model is slower to evaluate due to the larger number of energy
bins (101 versus 2). Material identification in the betatron experiments was
worse than for the cobalt-60 experiments, (the Al-Cu Shells correct assignment
moved from rank one to two and the Eight Cubes top assignments changed from
five of eight to three of eight correct). We believe that this difference in
accuracy stems from the betatron source being relatively unstable, making
measurements less repeatable and flat-field correction difficult. This effect
is apparent in Figure 1, where the background for the Bremsstrahlung datasets
is notably nonuniform even after flat-field correction. While this effect can
be partially mitigated by dividing by the open-beam value for each radiograph
individually, doing so does not affect our algorithm’s material assignments
because the scalar $\alpha$ in (1) makes the approach invariant to rescaling
the radiograph ${\bm{t}}$. It also appears in Figure 12, where the
transmission appears nonsymmetric, i.e., the minimum transmission on the left
side of the inner copper shell is higher than the minimum transmission on its
right side (although this effect may also be due to scatter or blur).

The efficiency of branch and bound in our experiments depended on the
arrangement of the scene. In the Al-Cu Shells scene, branch and bound was
moderately less efficient than exhaustive search would have been because none
of the partial solution bounds resulted in exclusion of any full solutions.
This is probably because the shells overlap significantly in the radiograph
and therefore choosing the material of the outer shell only determines a small
amount of the total error. In the extreme case where all objects overlap
completely in the radiograph, e.g., a set of sheets perpendicular to the
optical axis, our branch and bound algorithm would never be faster than
exhaustive search because all partial bounds would be zero. However, we note
that the overhead is modest (roughly log of the total number of solutions), so
that in situations where exhaustive search is computationally feasible, branch
and bound will be, also. On the other hand, our experiments on the Eight Cubes
scene demonstrate that our algorithm can provide tremendous time savings over
exhaustive search. This happens because there are many partial solutions that
are so inaccurate that no material assignment for the remaining objects can
make them viable. We also note that both our branch and bound algorithm and
exhaustive search are simple to parallelize, which should allow decreasing the
runtimes reported here by several orders of magnitude when running on a
moderately-sized cluster (100s of nodes).

Our experiments also showed that the accuracy of material identification
depends on the arrangement of objects: the algorithm performed better on
scenes with concentric shells than on the cubes data. We believe that the
cubes data is especially challenging because it presents no diversity of path
lengths. This means that identification must be performed on the basis of one
number (the transmission through the cube) per object. Practically speaking,
such arrangements of objects are unlikely because they require the objects to
be perfectly aligned with the detector. We suspect that performance would have
been better if the cubes had been rotated so as to break this alignment and
provide a range of path lengths. We also note that, because of the exponential
growth of the number of possible solutions, placing the ground truth in the
top $N$ best fits is more difficult when there are more objects. For example,
the top 20 assignments represent $0.25\%$ of the total possible assignments
when there are 3 objects, but only $7.8\times 10^{-8}\%$ of the total when
there are 8 objects. This effect may also explain why our method did not place
the correct material assignment into the top 20 on the cubes data.

Finally, our experiments demonstrate how challenging it is to accurately model
real radiographic data. For example, even in the relatively simple Al-Cu
Shells scene, our best model showed deviations as large as ten percent of the
measured signal (Figure 7). This mismatch points to the need for more
sophisticated scatter models and/or further system modeling and calibration.
However, our experiments also show that the method tolerates some level of
approximation in modeling the data.

##  V Limitations

We now discuss several limitations of the current work.

First, we require significant knowledge of the scene: a list of possible
materials, their attenuation coefficients, the source spectrum, the detector
response, and the object geometry. Attenuation coefficients for many materials
are freely available, e.g., see the National Institute of Standards and
Technology Photon Cross Sections Database [40]. The challenge of producing a
list of possible materials depends on the application, but one can imagine
security scenarios wherein only a handful of materials are relevant and the
rest can be safely covered by a few generic material definitions. We have
demonstrated that the source spectrum and detector response can be accounted
for via calibration, which represents a relatively small time investment
amortized over the life of the system (or time between calibrations, if done
repeatedly). The most restrictive assumption, then, is the known object
geometry. Solving this challenge in general is beyond the scope of this paper,
but we argue that it can be approached by performing a kind of tomographic
reconstruction using multiple radiographs. This tomographic problem should be
less challenging than a full tomographic reconstruction because only object
boundaries must be recovered, rather than a map of internal densities. In
settings where the shape of the objects is known but their 3D pose is not,
techniques from the computer vision community for camera calibration and pose
detection/registration [41] will prove useful.

Second, our model does not include source or detector blur, and we have not
performed the experiments needed to characterize the blur in our setups. Our
model could be extended to include blur (if properly characterized), however
this would likely incur significant computational cost. Another approach would
be to mitigate the effect of blur by excluding pixels near boundaries from the
analysis.

Finally, our experiments covered either relatively simple (two or three
spherical shells) or very challenging (eight axis-aligned cubes) scenarios.
So, while our study both provides a proof of concept and identifies a regime
where performance degrades, more experimentation is needed to elucidate the
performance of the approach in intermediate settings.

##  VI Conclusions and Future Work

We have presented a method for identifying materials from radiographs without
energy-resolved measurements, thus allowing standard X-ray systems to provide
material identification information without requiring additional hardware. To
solve the underlying combinatorial optimization problem, we developed a branch
and bound algorithm with a bounding function specifically tailored to the
radiographic setting. Our experiments on metallic objects with MeV X-rays show
that the method can successfully identify materials from real radiographs,
although its accuracy and runtime are scene-dependent. This method could prove
to be a useful new capability in a variety of security applications.

This method can be extended in several ways: We could explore other scatter
models (e.g., piecewise polynomials, smooth images) to allow finer control of
the representation power or try physically-motivated kernel-based scatter
models such as [33, 34]. We could add a technique to automatically estimate
scene geometry and/or pose (which would presumably require multiple
radiographs). Finally, we could seek improved initialization, branching, and
bounding algorithms to speed up the optimization process, allowing even larger
problems to be tackled efficiently.

## References

  * [1] E. M. Hussein and E. J. Waller, “Review of one-side approaches to radiographic imaging for detection of explosives and narcotics,” _Radiation Measurements_ , vol. 29, no. 6, pp. 581–591, Dec. 1998.
  * [2] R. D. R. Macdonald, “Design and implementation of a dual-energy x-ray imaging system for organic material detection in an airport security application,” in _SPIE Proceedings_ , M. A. Hunt, Ed. SPIE, Apr. 2001.
  * [3] O. E. Wetter, “Imaging in airport security: Past, present, future, and the link to forensic and clinical radiology,” _Journal of Forensic Radiology and Imaging_ , vol. 1, no. 4, pp. 152–160, Oct. 2013. [Online]. Available: https://doi.org/10.1016/j.jofri.2013.07.002
  * [4] L. Martin, A. Tuysuzoglu, W. C. Karl, and P. Ishwar, “Learning-based object identification and segmentation using dual-energy CT images for security,” _IEEE Transactions on Image Processing_ , vol. 24, no. 11, pp. 4069–4081, Nov. 2015.
  * [5] A. Mouton and T. P. Breckon, “Materials-based 3D segmentation of unknown objects from dual-energy computed tomography imagery in baggage security screening,” _Pattern Recognition_ , vol. 48, no. 6, pp. 1961–1978, Jun. 2015\.
  * [6] L. Zhang and T. YangDai, “Determination of liquid’s molecular interference function based on x-ray diffraction and dual-energy CT in security screening,” _Applied Radiation and Isotopes_ , vol. 114, pp. 179–187, Aug. 2016.
  * [7] A. Chalmers, “Cargo identification algorithms facilitating unmanned/unattended inspection at high throughput portals,” in _SPIE Proceedings_ , E. M. Carapezza, Ed. SPIE, Oct. 2007. [Online]. Available: https://doi.org/10.1117/12.738354
  * [8] Y. Yang, B. Wu, T. Li, and Y. Li, “Materials identification by X-ray and photoneutron transmission,” in _2009 IEEE Nuclear Science Symposium Conference Record (NSS/MIC)_. IEEE, Oct. 2009.
  * [9] P. M. Shikhaliev, “Megavoltage cargo radiography with dual energy material decomposition,” _Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment_ , vol. 882, pp. 158–168, Feb. 2018.
  * [10] D. Lee, J. Lee, J. Min, B. Lee, B. Lee, K. Oh, J. Kim, and S. Cho, “Efficient material decomposition method for dual-energy X-ray cargo inspection system,” _Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment_ , vol. 884, pp. 105–112, Mar. 2018.
  * [11] C. Hill, I. Rutel, and F. Cain, “SU-e-i-105: Optimizing bomb squad x-ray systems for incidental human exposures,” _Medical Physics_ , vol. 39, no. 6Part5, pp. 3649–3649, Jun. 2012. [Online]. Available: https://doi.org/10.1118/1.4734822
  * [12] R. G. Beckett and G. J. Conlogue, _Advances in Paleoimaging_. CRC Press, Aug. 2020. [Online]. Available: https://doi.org/10.4324/9781315203089
  * [13] S. Akcay and T. Breckon, “Towards automatic threat detection: A survey of advances of deep learning within x-ray security imaging,” _Pattern Recognition_ , vol. 122, p. 108245, Feb. 2022. [Online]. Available: https://doi.org/10.1016/j.patcog.2021.108245
  * [14] N. Kimoto, H. Hayashi, T. Asahara, Y. Mihara, Y. Kanazawa, T. Yamakawa, S. Yamamoto, M. Yamasaki, and M. Okada, “Precise material identification method based on a photon counting technique with correction of the beam hardening effect in X-ray spectra,” _Applied Radiation and Isotopes_ , vol. 124, pp. 16–26, Jun. 2017.
  * [15] P. Sukovle and N. Clinthorne, “Basis material decomposition using triple-energy x-ray computed tomography,” in _IMTC/99. Proceedings of the 16th IEEE Instrumentation and Measurement Technology Conference (Cat. No.99CH36309)_. IEEE, 1999. [Online]. Available: https://doi.org/10.1109/imtc.1999.776097
  * [16] M. N. Lakshmanan, A. J. Kapadia, P. Sahbaee, S. D. Wolter, B. P. Harrawood, D. Brady, and E. Samei, “An X-ray scatter system for material identification in cluttered objects: A Monte Carlo simulation study,” _Nuclear Instruments and Methods in Physics Research Section B: Beam Interactions with Materials and Atoms_ , vol. 335, pp. 31–38, Sep. 2014.
  * [17] Z. Ying, R. Naidu, and C. R. Crawford, “Dual energy computed tomography forexplosive detection,” _Journal of X-Ray Science and Technology_ , vol. 14, no. 4, pp. 235–256, Sep. 2006.
  * [18] R. C. Bokun and S. M. Osadchii, “Generalization of dual-energy method for material identification with X-ray introscopy,” _Journal of Surface Investigation. X-ray, Synchrotron and Neutron Techniques_ , vol. 4, no. 4, pp. 591–593, Aug. 2010.
  * [19] Y. Long and J. A. Fessler, “Multi-material decomposition using statistical image reconstruction for spectral CT,” _IEEE Transactions on Medical Imaging_ , vol. 33, no. 8, pp. 1614–1626, Aug. 2014.
  * [20] D. Modgil, D. S. Rigie, Y. Wang, X. Xiao, P. A. Vargas, and P. J. L. Rivière, “Material identification in x-ray microscopy and micro CT using multi-layer, multi-color scintillation detectors,” _Physics in Medicine and Biology_ , vol. 60, no. 20, pp. 8025–8045, Sep. 2015.
  * [21] J. Duvillier, M. Dierick, J. Dhaene, D. V. Loo, B. Masschaele, R. Geurts, L. V. Hoorebeke, and M. N. Boone, “Inline multi-material identification via dual energy radiographic measurements,” _NDT & E International_, vol. 94, pp. 120–125, Mar. 2018.
  * [22] G. Beldjoudi, V. Rebuffel, L. Verger, V. Kaftandjian, and J. Rinkel, “An optimised method for material identification using a photon counting detector,” _Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment_ , vol. 663, no. 1, pp. 26–36, 2012. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S0168900211017451
  * [23] X. Wu, Q. Wang, J. Ma, W. Zhang, P. Li, and Z. Fang, “A hyperspectral X-ray computed tomography system for enhanced material identification,” _Review of Scientific Instruments_ , vol. 88, no. 8, p. 083111, Aug. 2017\.
  * [24] A. Macovski, R. Alvarez, J.-H. Chan, J. Stonestrom, and L. Zatz, “Energy dependent reconstruction in X-ray computerized tomography,” _Computers in Biology and Medicine_ , vol. 6, no. 4, pp. 325–336, Oct. 1976\.
  * [25] Y. Yuan, B. Tracey, and E. Miller, “Robust x-ray based material identification using multi-energy sinogram decomposition,” in _Anomaly Detection and Imaging with X-Rays (ADIX)_ , A. Ashok, M. A. Neifeld, and M. E. Gehm, Eds., vol. 9847, International Society for Optics and Photonics. SPIE, 2016, pp. 207 – 219.
  * [26] B. Tobias and M. Klasky, “Limits of special material detectability fundamental to idealized dual-energy radiographic systems,” _Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment_ , vol. 944, p. 162563, Nov. 2019.
  * [27] E. Braig, J. Böhm, M. Dierolf, C. Jud, B. Günther, K. Mechlem, S. Allner, T. Sellerer, K. Achterhold, B. Gleich, P. Noël, D. Pfeiffer, E. Rummeny, J. Herzen, and F. Pfeiffer, “Direct quantitative material decomposition employing grating-based x-ray phase-contrast CT,” _Scientific Reports_ , vol. 8, no. 1, Nov. 2018. [Online]. Available: https://doi.org/10.1038/s41598-018-34809-6
  * [28] F. Schaff, K. S. Morgan, J. A. Pollock, L. C. P. Croton, S. B. Hooper, and M. J. Kitchen, “Material decomposition using spectral propagation-based phase-contrast X-ray imaging,” _IEEE Transactions on Medical Imaging_ , vol. 39, no. 12, pp. 3891–3899, Dec. 2020.
  * [29] J. Jones, Y. Harker, W. Yoon, and L. Johnson, “Material identification technology (MIT) concept technical feasibility study,” Westinghouse Idaho Nuclear Co., Inc., Idaho Falls, ID, Tech. Rep., Sep. 1993.
  * [30] T. Cui, Y. Yu, Y. Yang, Z. Zhang, and X. Wang, “Material identification using dual particle interrogation,” _Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment_ , vol. 954, p. 161827, Feb. 2020.
  * [31] A. Kak and M. Slaney, _Principles of Computerized Tomographic Imaging_ , ser. Classics in Applied Mathematics. New York, NY: Society for Industrial and Applied Mathematics, 2001.
  * [32] M. T. McCann, M. L. Klasky, J. L. Schei, and S. Ravishankar, “Local models for scatter estimation and descattering in polyenergetic X-ray tomography,” _Opt. Express_ , vol. 29, no. 18, pp. 29 423–29 438, Aug 2021.
  * [33] V. N. Hansen, W. Swindell, and P. M. Evans, “Extraction of primary signal from EPIDs using only forward convolution,” _Medical Physics_ , vol. 24, no. 9, pp. 1477–1484, Sep. 1997. [Online]. Available: https://doi.org/10.1118/1.598036
  * [34] B. Ohnesorge, T. Flohr, and K. Klingenbeck-Regn, “Efficient object scatter correction algorithm for third and fourth generation CT scanners,” _European Radiology_ , vol. 9, no. 3, pp. 563–569, Mar. 1999.
  * [35] R. G. Parker and R. L. Rardin, _Discrete Optimization_. San Diego, CA: Academic Press, 1988.
  * [36] Y. Tao, V. Hristidis, D. Papadias, and Y. Papakonstantinou, “Branch-and-bound processing of ranked queries,” _Information Systems_ , vol. 32, no. 3, pp. 424–445, May 2007.
  * [37] P. Poensgen and R. Möller, “Branch-and-bound ranked search by minimizing parabolic polynomials,” _Open Journal of Databases (OJDB)_ , vol. 7, no. 1, pp. 12–20, 2020.
  * [38] J. Bradbury, R. Frostig, P. Hawkins, M. J. Johnson, C. Leary, D. Maclaurin, G. Necula, A. Paszke, J. VanderPlas, S. Wanderman-Milne _et al._ , “JAX: Autograd and XLA,” _Astrophysics Source Code Library_ , pp. ascl–2111, Nov. 2021.
  * [39] C. J. Werner, J. S. Bull, C. J. Solomon, F. Brown, G. McKinney, M. Rising, D. Dixon, R. Martz, H. Hughes, L. Cox, A. J. Zukaitis, J. Armstrong, R. A. Forster, and L. Casswell, “MCNP version 6.2 release notes,” Los Alamos National Laboratory, Tech. Rep. LA-UR-18-20808, 2018.
  * [40] M. Berger, J. Hubbell, S. Seltzer, J. Chang, J. Coursey, R. Sukumar, and K. Zucker, D.S. an Olsen. XCOM: Photon cross section database (version 1.5). [Online]. Available: http://physics.nist.gov/xcom
  * [41] D. A. Forsyth and J. Ponce, _Computer Vision_ , 2nd ed. Upper Saddle River, NJ: Pearson, Oct. 2011.