# Local Models for Scatter Estimation and Descattering in Polyenergetic X-Ray
Tomography

Michael T. McCann  1,2,* Marc L. Klasky  2 Jennifer L. Schei  3 and Saiprasad
Ravishankar1,4  1 Department of Computational Mathematics, Science and
Engineering, Michigan State University, East Lansing, MI, 48824, USA
(previously)
2 Theoretical Division, Los Alamos National Laboratory, Los Alamos, NM, 87545,
USA
3 Physics Division, Los Alamos National Laboratory, Los Alamos, NM, 87545, USA
4 Department of Biomedical Engineering, Michigan State University, East
Lansing, MI, 48824 USA<EMAIL_ADDRESS>
# Local Models for Scatter Estimation and Descattering in Polyenergetic X-Ray
Tomography: supplemental document

Michael T. McCann  1,2,* Marc L. Klasky  2 Jennifer L. Schei  3 and Saiprasad
Ravishankar1,4  1 Department of Computational Mathematics, Science and
Engineering, Michigan State University, East Lansing, MI, 48824, USA
(previously)
2 Theoretical Division, Los Alamos National Laboratory, Los Alamos, NM, 87545,
USA
3 Physics Division, Los Alamos National Laboratory, Los Alamos, NM, 87545, USA
4 Department of Biomedical Engineering, Michigan State University, East
Lansing, MI, 48824 USA<EMAIL_ADDRESS>
###### Abstract

We propose a new modeling approach for scatter estimation and descattering in
polyenergetic X-ray computed tomography (CT) based on fitting models to local
neighborhoods of a training set. X-ray CT is widely used in medical and
industrial applications. X-ray scatter, if not accounted for during
reconstruction, creates a loss of contrast in CT reconstructions and
introduces severe artifacts including cupping, shading, and streaks. Even when
these qualitative artifacts are not apparent, scatter can pose a major
obstacle in obtaining quantitatively accurate reconstructions. Our approach to
estimating scatter is, first, to generate a training set of 2D radiographs
with and without scatter using particle transport simulation software. To
estimate scatter for a new radiograph, we adaptively fit a scatter model to a
small subset of the training data containing the radiographs most similar to
it. We compared local and global (fit on full data sets) versions of several
X-ray scatter models, including two from the recent literature, as well as a
recent deep learning-based scatter model, in the context of descattering and
quantitative density reconstruction of simulated, spherically symmetrical,
single-material objects comprising shells of various densities. Our results
show that, when applied locally, even simple models provide state-of-the-art
descattering, reducing the error in density reconstruction due to scatter by
more than half.

###### Abstract

This document contains supplemental methods and figures for “Local Models for
Scatter Estimation and Descattering in Polyenergetic X-Ray Tomography.” For
the complete context, refer to the main text, Sections 3.5 and 4.3.

††journal: oe

##  1 Introduction

In the event of a possible nuclear threat, the United States government
employs an Emergency Response Team to determine the proper course of action. A
crucial step in assessing the threat is a determination of the object
characteristics, and radiographic examination of the object is an important
part of this process. While tomographic reconstruction techniques such as
filtered backprojection or model-based iterative reconstruction (MBIR) [1] are
available, great difficulty remains in determining the material composition of
the object due to scattered radiation, spectral effects, and issues associated
with the detector. Similar challenges exist in the context of baggage
screening and nondestructive testing.

This work focuses on addressing scattered radiation in service of performing
quantitative radiographic reconstructions. In general, scatter creates loss of
contrast impairing the ability to find sharp features as well as introducing
severe image artifacts, e.g., cupping, shading, streaks, etc. Consequently, in
applications in which a quantitative density reconstruction is desired,
scatter must be mitigated in hardware (e.g., by increasing the distance
between the object and the image plane, or through use of a Bucky-Potter grid)
and/or removed in software during reconstruction.

Scatter is a complex phenomenon with many possible types of individual
radiation-matter interactions (as X-ray photons pass through an object), e.g.,
Compton scatter, Rayliegh scatter, pair production, and the photoelectric
effect [2]. The total amount of scatter along a ray is simply proportional to
the probability of interaction [2], but the difficulty lies in computing the
spatial distribution of scatter reaching the detector. The transport of the
scattered radiation is in a direction away from the direction of the initial
scattering interaction. Consequently, the resulting scatter that contributes
to the image is a function of the magnitude of the scattered radiation
produced, its angular distribution and the density distribution of the object.
This motivates models in which a convolutional kernel is utilized to capture
the spatial aspects of scatter, but where the kernel somehow changes in
response to the density distribution of the object, e.g., [3].

There exists a vast literature on scatter correction, beginning in at least
1976 [4] and presented in the context of many different scanners and
applications, from keV CT for medical imaging to MeV CT for nondestructive
testing (where keV and MeV refer to the nominal energy of the X-ray photons).
Here, we aim to discuss the main trends in computational approaches to scatter
correction, see the pair of reviews [5, 6] for more comprehensive coverage.

Most scatter correction techniques rely on scatter estimation, and can
therefore be classified by what information they estimate the scatter from,
e.g., the total transmission radiograph, a radiograph without scatter (which
we refer to as a direct radiograph), or the reconstructed density [5]. Scatter
estimation from a reconstructed density volume can be highly accurate,
because, in principle, the density fully determines the scatter. In practice,
it is slow because it necessitates alternating between approximate
reconstruction of the 3D density and scatter estimation and removal. Scatter
estimation from the total transmission (i.e., the radiographic measurements
corrupted by scatter) has the advantage of providing single-shot descattering:
no iteration is required; however, these methods have limited accuracy because
some approximation is involved when predicting scatter without access to the
underlying density. Scatter estimation from a radiograph without scatter
strikes a balance between accuracy and speed, requiring only iteration in the
radiograph domain (the scatter is repeatedly estimated and removed based on
iteratively estimated direct radiographs) thereby avoiding the need to perform
multiple high-dimensional density reconstructions. Which scatter estimation
approach is best depends on the application: the amount and dominant mechanism
of scatter, the requirements on descattering and reconstruction time, and the
desired accuracy of the reconstructions.

In this paper, we propose a locally learned model of scatter, wherein the free
parameters of a given model are fit on training data in a local neighborhood
around the 2D radiograph to be descattered. We emphasize that the proposed
model is local in the sense that the model changes as a function of the
input—as distinct from spatially local models, where the model changes as a
function of location in the radiograph.

We use the proposed local scatter model as a component of an algorithm to
reconstruct the material density field of an object (with known material
composition) from radiographs obtained with a polyenergetic X-ray beam. Our
focus on _quantitative_ material density reconstruction is somewhat different
than in typical medical applications, where the material composition of the
object being imaged (and therefore the physical properties that would allow
mapping linear attenuation to density) is unknown and therefore typically only
the linear attenuation map is reconstructed (with some exceptions, e.g., [7]).
The use of a polyenergetic beam introduces another complication to our problem
as compared to medical imaging, where often (although not always, e.g., [8]) a
monoenergetic source is assumed. In our setting, we assume known materials,
which gives us the opportunity to reconstruct density without facing an
underdetermined problem. Future work will examine objects consisting of
multiple material compositions.

The rest of the manuscript is organized as follows. The remainder of this
section discusses related work. In Section 2, we discuss the descattering
problem, including the forward model for monoenergetic and polyenergetic
tomography. In Section 3, we discuss local modeling in general, then describe
local models for scatter estimation and our descattering algorithm. In Section
4, we present experiments validating our proposed local modeling approach. We
provide additional discussion in Section 5 and conclude in Section 6.

##  2 Monoenergetic and polyenergetic tomography with scatter

The overall goal of this work is to recover a discretized, space-varying
density ${\bm{\rho}}\in\mathbb{R}^{N_{1}\times N_{2}\times N_{3}}$, from a set
of $V$ X-ray projection measurements (radiographs), ${\bm{t}}_{1}$,
${\bm{t}}_{2}$, …, ${\bm{t}}_{V}\in\mathbb{R}^{M_{1}\times M_{2}}$, where $v$
indexes the viewing direction. (If the object being imaged is axially
symmetric, as is the case in some important security applications, a single
view suffices for reconstruction, since, ignoring noise, all views are
identical). In this section, we describe the forward model that relates
${\bm{\rho}}$ to $\\{{\bm{t}}_{v}\\}_{v=1}^{V}$; much of our presentation is
adapted from the standard textbook [9].

X-ray
---
source
ray $r_{m,n}$
---
$z$
---
$x$
---
$y$
---
detector
---
${\bm{\rho}}$
---
pixel
---
$(m,n)$

$R_{m,n}$
---

Figure 1: Sketch of the tomography setup. In reality, the object is 3D and the
detector is a plane. ${\bm{\rho}}\in\mathbb{R}^{N_{1}\times N_{2}\times
N_{3}}$ denotes the density field of the object being imaged, $R_{m,n}$ is the
set of rays entering detector pixel $(m,n)$, and $r_{m,n}$ is the central ray
of $R_{m,n}$.

In a generic radiographic setup (see Figure 1 for a sketch), there is an X-ray
source, one or more objects being imaged, and a detector. Every ray between
the source and detector, $r$, intersects a certain amount of each object. This
is called the areal density of the object along $r$,

${\rho_{A}}_{i}(r)=\int_{-\infty}^{\infty}\rho_{i}(r_{x}(t),r_{y}(t),r_{z}(t))dt,$
(1)

where ${\rho_{A}}_{i}(r)$ has units g/cm2, $r_{x}(t),r_{y}(t),r_{z}(t)$ are a
parameterization of the ray $r$, and $i$ is the object index—each different
material being imaged corresponds to a different $i$. We note a slight abuse
of notation here: a continuous $\rho$ in (1) versus a discrete ${\bm{\rho}}$
mentioned before; we will put aside the details of discretization here as they
are an implementation detail of the reconstruction software discussed in
Section 3.5.

###  2.1 Monoenergetic case

In the absence of scatter and with a monoenergetic source, the number density
of photons reaching the detector along ray $r$ is

${q}(r)={q}_{\text{in}}(r)\exp\left(-\sum_{i}\xi_{i}{\rho_{A}}_{i}(r)\right),$
(2)

where ${q}(r)$ has units 1/cm2, ${q}_{\text{in}}(r)$ is the number density of
the incident beam along ray $r$, and $\xi_{i}\in\mathbb{R}$ is the mass
attenuation coefficient of the $i$th material, with units cm2/g. In our
application, we assume known materials, therefore the mass attenuation
coefficients can be looked up, e.g., in the XCOM database [10].

In reality, what we measure is not the number density of photons at
infinitesimal points, as given by (2), but a discrete, finite image which
represents the number of photons reaching each of a finite grid of detectors.
We refer to this scatter-free measurement data as the direct radiograph,
${\bm{d}}\in\mathbb{R}^{M_{1}\times M_{2}}$. It is given by

${\bm{d}}[m,n]=\int_{R_{m,n}}{q}(r)dr\approx C{q}(r_{m,n})$ (3)

where ${\bm{d}}[m,n]$ is unitless, $R_{m,n}$ is the set of rays entering
detector pixel $(m,n)$, $r_{m,n}$ is the ray through the center of pixel
$(m,n)$, and $C$ is a constant related to the area of a detector pixel and
angle of incidence; the approximation holds when the detector elements are
sufficiently small and regular in size (think pixels). Real measurement
systems also have a per-pixel gain, meaning that we actually measure
${\bm{c}}[m,n]{\bm{d}}[m,n]$. We ignore this gain factor as it is easily
handled via calibration.

###  2.2 Polyenergetic case

In the polyenergetic case, we have

${q}(r)=\int_{0}^{\infty}{q}_{\text{in}}(r,E)\exp\left(-\sum_{i}\xi_{i}(E){\rho_{A}}_{i}(r)\right)dE,$
(4)

where ${q}_{\text{in}}(r,E)$ is the energy profile of the number density of
the incident beam and $\xi_{i}(E)$ gives the mass attenuation coefficient of
the $i$th material at energy $E$. Using the same numerical integration
approximation as in the monoenergetic case and an additional quadrature for
energy gives

${\bm{d}}[m,n]\approx\sum_{e}C{q}(r_{m,n},E_{e}),$ (5)

where $e$ indexes discrete energy bins, $C$ is a constant depending on the
pixel size and energy bin width, and

${q}(r,E)={q}_{\text{in}}(r,E)\exp\left(-\sum_{i}\xi_{i}(E){\rho_{A}}_{i}(r)\right).$
(6)

$\displaystyle-5$$\displaystyle 0$$\displaystyle 5$y
[cm]$\displaystyle-5.0$$\displaystyle-2.5$$\displaystyle 0.0$$\displaystyle
2.5$$\displaystyle 5.0$x [cm]$\displaystyle 0.5$$\displaystyle
1.0$$\displaystyle 1.5$$\displaystyle 2.0$density [g/cc]$\displaystyle\times
10^{1}$

$\displaystyle-5$$\displaystyle 0$$\displaystyle 5$y
[cm]$\displaystyle-5.0$$\displaystyle-2.5$$\displaystyle 0.0$$\displaystyle
2.5$$\displaystyle 5.0$x [cm]$\displaystyle 1$$\displaystyle 2$$\displaystyle
3$transmission$\displaystyle\times 10^{-5}$

$\displaystyle 0$$\displaystyle 1$$\displaystyle 2$$\displaystyle
3$$\displaystyle 4$$\displaystyle 5$x [cm]$\displaystyle 0$$\displaystyle
5$$\displaystyle 10$$\displaystyle 15$$\displaystyle 20$density
[g/cc]truthscatterno scatter

Figure 2: Illustration of how the presence of scatter negatively impacts
quantitative density reconstruction (monoenergetic dataset). Left: Cross
section showing the density of a spherically symmetric, 3D object. Right:
Radiographic measurement of the object on the left, which includes scatter.
Our goal is to reconstruct the left image from the center one. Bottom: 1D
profiles of the true density and reconstruction from radiographs with and
without scatter. Without scatter, reconstruction is nearly perfect (orange,
dotted line overlaps black); with scatter, reconstruction involves density
errors on the order of 1 g/cm3.

In some cases, it is sufficient to stop at the model for the direct radiograph
given by (3) or (5). However, in many cases, there is a nontrivial amount of
scatter, i.e., photons that reach the detector without following a straight
path through the object (Figure 2 illustrates the effect this scatter on
reconstruction). Scatter arises from several different physical processes. A
photon may undergo exactly one coherent scattering interaction with the
object; undergo exactly one Compton scattering interaction with the object;
interact with and scatter from parts of the measurement apparatus, e.g., the
detector itself [11] or the room; or undergo any other interaction(s) not
included above, e.g., multiple scatters or pair production. These processes
are well understood, but they are complicated enough that they need to be
modeled at the individual photon level, a computationally expensive process
and one that is difficult to invert. For our purposes, we consider the
measured radiograph ${\bm{t}}\in\mathbb{R}^{M_{1}\times M_{2}}$ as a sum of
direct photons, ${\bm{d}}$ as described in (3) or (5), and scattered photons,
${\bm{s}}\in\mathbb{R}^{M_{1}\times M_{2}}$, ${\bm{t}}={\bm{d}}+{\bm{s}}$. We
will also refer to ${\bm{t}}$ as the total radiograph, since it contains both
the contribution of the scattered and direct photons. Finding a way to
efficiently and accurately model ${\bm{s}}$ is the main focus of this work.

##  3 Proposed approach—local modeling for descattering

As stated earlier, the goal of this work is to recover a 3D object,
${\bm{\rho}}\in\mathbb{R}^{N_{1}\times N_{2}\times N_{3}}$, from its 2D
tomographic measurements (corrupted by scatter), ${\bm{t}}_{1}$,
${\bm{t}}_{2}$, …, ${\bm{t}}_{V}\in\mathbb{R}^{M_{1}\times M_{2}}$. Our
proposed approach consists of two steps: descattering radiographs and
tomographic reconstruction; we now discuss the first of these. We discuss
related work, introduce our proposed local modeling approach, describe the
specific local models we evaluate, and finish with a description of our
descattering algorithm.

###  3.1 Related work on scatter estimation

Most methods that estimate scatter from a direct radiograph use some form of
the model

${\bm{s}}[m,n]=\sum_{m^{\prime},n^{\prime}}f({\bm{d}}[m^{\prime},n^{\prime}])\bm{H}_{\bm{\rho}}[m^{\prime},n^{\prime},m,n],$
(7)

where ${\bm{s}}[m,n]$ and ${\bm{d}}[m,n]$ are the scatter (only) and direct
(without scatter) radiographs at position $m,n$ and $f$ is a scalar function.
The term $f({\bm{d}}[m,n])$, sometimes called the _scatter potential_ or
_scatter amplitude factor_ , represents how much scatter comes off of a given
pixel. The term $\bm{H}_{\bm{\rho}}[m^{\prime},n^{\prime},m,n]$, called the
_scatter kernel_ , represents how much radiation scatters from position
$(m^{\prime},n^{\prime})$ to position $(m,n)$. In the most general case, it is
dependent on the 3D density field of the object being imaged, ${\bm{\rho}}$.

The model (7) is hardly useful as it is because it is so general. The
challenge is simplifying it to a point where it can be parameterized and the
parameters estimated. One common approach is to fix the scatter amplitude
factor to the parametric form $f(x)=x^{\alpha}|\ln x|^{\beta}$, which can be
motivated by a simple physical model [12]. Another is to constrain the scatter
kernel to be shift invariant,
$\bm{H}_{\bm{\rho}}[m^{\prime},n^{\prime},m,n]=\bm{H}_{\bm{\rho}}[m^{\prime}-m,n^{\prime}-n]$,
or a spatially-varying weighted sum of such kernels, where the weights may
depend on (estimated) object thickness or distance to the object boundary [3].
There also exists recent work that uses a convolutional neural network to map
from a direct radiograph to the corresponding scatter radiograph [13], which
we can view as an overparameterized and nonlinear version of (7).

No matter the form of the scatter model, its parameters must be determined in
order to perform scatter estimation. One approach is to fit the parameters on
simple simulations, e.g., of a pencil beam through a slab [3, 14]. Another is
to approach fitting as an optimization problem on complex reference
simulation(s) [13]. The benefit of the former approach is that it may
generalize well and fits well with physically interpretable models. The
benefit of the latter is that it may provide more accurate fits, because it
explicitly minimizes fit error.

###  3.2 Local modeling formulation

Given a training set of 2D direct radiographs and the corresponding 2D scatter
radiographs, $\\{({\bm{d}}_{t},{\bm{s}}_{t})\\}_{t=1}^{T}$, our goal is to
learn a mapping, $\mathcal{A}:\mathbb{R}^{M_{1}\times
M_{2}}\to\mathbb{R}^{M_{1}\times M_{2}}$, from direct to scatter. (We leave
aside the question of whether such a mapping exists, i.e., whether two objects
can produce the same direct radiograph but different scatter radiographs.)
These training radiographs can be formed by generating synthetic objects,
$\rho_{t}$, and simulating their radiographs via particle transport software,
e.g., MCNP [15] or Geant4 [16].See Section 4.2 for details of the training
sets we used in our experiments.

Our proposed approach is to represent $\mathcal{A}$ via a large number of
simple models, i.e., a collection of $\mathcal{A}_{{\bm{d}}}$’s, each acting
in a different region of $\mathbb{R}^{M_{1}\times M_{2}}$. To do this, we
define the nearest neighbor function, $\mathbb{U}$, to return the set of
indices of the $G$ nearest training pairs to ${\bm{d}}$ in the training set,
i.e., $\mathbb{U}({\bm{d}})=\\{t_{1},t_{2},\dots t_{G}\\},$ such that
$\|{\bm{d}}-{\bm{d}}_{s}\|_{F}^{2}\leq\|{\bm{d}}-{\bm{d}}_{t}\|_{F}^{2}$ for
$s\in\mathbb{U}({\bm{d}})$ and $t\not\in\mathbb{U}({\bm{d}})$, where
$\|\cdot\|_{F}^{2}$ denotes the squared Frobenius norm (sum of the squared
elements). (Other metrics may be used to determine these neighborhoods, e.g.,
the Wasserstein metric, if they are more suitable to the particular
application) Finally, we select a class, $\mathcal{H}$, of models for the
local $\mathcal{A}_{{\bm{d}}}$’s; see Section 3.3 for specific examples.

Our local model $\mathcal{A}$ is then given by

$\mathcal{A}({\bm{d}})=\mathcal{A}_{{\bm{d}}}({\bm{d}}),\quad\text{where}\quad\mathcal{A}_{{\bm{d}}}=\operatorname*{argmin}_{\mathcal{B}\in\mathcal{H}}\sum_{t\in\mathbb{U}({\bm{d}})}\|\mathcal{B}({\bm{d}}_{t})-{\bm{s}}_{t}\|_{F}^{2}.$
(8)

Note the use of a subscript ${\bm{d}}$: $\mathcal{A}_{\bm{d}}$ is a model
(i.e., an image-to-image function) fit on the neighbors of ${\bm{d}}$;
$\mathcal{A}_{\bm{d}}({\bm{d}})$ is that model evaluated at the point
${\bm{d}}$, which is an image. Put simply: to predict ${\bm{s}}$ from
${\bm{d}}$, we fit a model on the nearest neighbors to ${\bm{d}}$ in the
training set and apply that model to ${\bm{d}}$.

What we have described is distinct from the standard approach to fitting a
model on a training set, which we refer to here as _global_ modeling. Using
our formulation, global modeling corresponds to local modeling with the number
of neighbors, $G$, set to the entire training set size, $T$. This simplifies
to

$\mathcal{A}_{{\bm{d}}}=\mathcal{A}=\operatorname*{argmin}_{\mathcal{B}\in\mathcal{H}}\sum_{t=1}^{T}\|\mathcal{B}({\bm{d}}_{t})-{\bm{s}}_{t}\|_{F}^{2}.$
(9)

To see the conceptual advantage of local over global models, consider the case
of linear regression: global modeling fits a single hyperplane though all of
the data, while local modeling fits different hyperplanes in different regions
of the space. Thus, local modeling can accurately describe more complex
relationships than global modeling, with the drawback that it may lead to
worse overfitting.

Related work on local modeling. Local modeling has a long history in a variety
of fields. For example, in statistics, local regression (also called LOESS),
which involves fitting local models to form a smooth interpolator, has been
studied since at least the 1970s [17]. The novelty of our work is to propose
and validate specific, high-dimensional local models for use in the context of
X-ray scatter estimation, descattering, and density reconstruction.

###  3.3 Models

We now describe the four classes of models (i.e., choices of $\mathcal{H}$ in
(8) and (9)) used in our experiments. The first two (single field and
convolutional) are simple models that we propose; the second two (parametric
and multikernel) come from the scatter modeling literature [3].

####  3.3.1 Single field model

The single field model represents the scatter as being a single image,
$\mathcal{A}_{\bm{d}}({\bm{d}})=\hat{{\bm{s}}}_{\bm{d}}\in\mathbb{R}^{M_{1}\times
M_{2}}$. Fitting the single field model requires solving the problem

$\hat{{\bm{s}}}_{\bm{d}}=\operatorname*{argmin}_{{\bm{s}}\in\mathbb{R}^{M_{1}\times
M_{2}}}\sum_{t\in\mathbb{X}}\|{\bm{s}}-{\bm{s}}_{t}\|^{2}_{F},$ (10)

where $\mathbb{X}$ may represent either the entire training set (_global
fitting_) or $\mathbb{U}({\bm{d}})$ (_local fitting_). The solution to (10) is
easily obtained by taking the pixelwise average of the
$\\{{\bm{s}}_{t}\\}_{t\in\mathbb{X}}$. In the case where one neighbor is used,
the single field model is the classical nearest neighbor interpolator. This
model provides a simple baseline for comparison.

####  3.3.2 Convolutional model

The convolution model represents the scatter as a convolution, which is a
simple and computational efficient way of modeling the spatial distribution of
scatter. The form of the model is $\mathcal{A}({\bm{d}})={\bm{k}}_{\bm{d}}\ast
f_{1,1}({\bm{d}})$, where

$f_{\alpha,\beta}({\bm{d}})[m_{1},m_{2}]=-{\bm{d}}[m1,m2]^{\alpha}(\log{\bm{d}}[m_{1},m_{2}])^{\beta}.$
(11)

Here, $\ast$ denotes the 2D, linear convolution between zero-padded images,
where the result is cropped to match the support of ${\bm{d}}$. The size of
the kernel is the smallest such that each point in ${\bm{d}}$ can contribute
to each point in the result, thus
${\bm{k}}_{\bm{d}}\in\mathbb{R}^{(2M_{1}-1)\times(2M_{2}-1)}$. Each pixel in
the kernel is arbitrary, thus, the number of parameters for the convolutional
model is $(2M_{1}-1)\times(2M_{2}-1)$. The nonlinear function $f$ is common in
the descattering literature and can be justified by a physics argument when
$\alpha=\beta=1$ [12]. (To make sense of the negative sign, note that
${\bm{d}}$ is expected to be normalized to lie between zero and one.) The
flexibility provided by $\alpha$ and $\beta$ is used in the next model.

The kernel is determined by solving the optimization problem

${\bm{k}}_{\bm{d}}=\operatorname*{argmin}_{\bm{k}}\sum_{t\in\mathbb{X}}\|{\bm{k}}\ast
f_{1,1}({\bm{d}}_{t})-{\bm{s}}_{t}\|_{F}^{2}$ (12)

where $\mathbb{X}$ is an appropriate set of training indices (local or
global). The problem (12) is a linear least squares problem because the
$f_{1,1}({\bm{d}}_{t})$’s are fixed. For our experiments, we solved it using
the conjugate gradient method [18].

####  3.3.3 Parametric kernel model

The parametric kernel model is a version of the convolution model with a
parametric kernel and nonlinearity. It appears in the recent literature on
descattering, e.g., [3]. The form of the model is
$\mathcal{A}({\bm{d}})={\bm{k}}_{A,B,\sigma_{1},\sigma_{2}}\ast
f_{\alpha,\beta}({\bm{d}}),$ where $f$ is defined as in (11) and all of the
parameters, $A$, $B$, $\sigma_{1}$, $\sigma_{2}$, $\alpha$, and $\beta$ depend
on ${\bm{d}}$. The kernel ${\bm{k}}$ is a sum of two centered Gaussian
functions,

${\bm{k}}_{A,B,\sigma_{1},\sigma_{2}}[m,n]=\frac{A}{\sigma_{1}\sqrt{2\pi}}e^{\left(-\frac{(m-m_{0})^{2}+(n-n_{0})^{2}}{2\sigma_{1}^{2}}\right)}+\\\
\frac{B}{\sigma_{2}\sqrt{2\pi}}e^{\left(-\frac{(m-m_{0})^{2}+(n-n_{0})^{2}}{2\sigma_{2}^{2}}\right)},$
(13)

where $m_{0},n_{0}$ gives the central pixel of the radiograph. Fitting the
parametric kernel model involves the following nonconvex optimization with
respect to the scalar parameters $A,B,\sigma_{1},\sigma_{2},\alpha$, and
$\beta$,

$\operatorname*{argmin}_{A,B,\sigma_{1},\sigma_{2},\alpha,\beta}\sum_{t\in\mathbb{X}}\|\mathcal{A}({\bm{d}})-{\bm{s}}_{t}\|^{2}_{F}.$
(14)

For our experiments, we used the LBFGS algorithm [19] (using the PyTorch
implementation), which we found to be significantly faster than gradient
descent. We note that LBFGS is not guaranteed to find the global minimizer of
$\eqref{eq:parametric_fit}$ because it is nonconvex.

####  3.3.4 Multikernel model

The multikernel model, proposed in [3], is a generalization of the parametric
kernel model wherein different parameters are used in different regions of the
image. To fit $K$ kernels, $K-1$ thresholds,
$\tau_{1},\tau_{2},\dots,\tau_{K-1}$ are established and a partition function
is defined according to

$\bm{r}_{\bm{d}}[m,n]=\begin{cases}1&0\leq-\log({\bm{d}}[m,n])\leq\tau_{1}\\\
2&\tau_{1}<-\log({\bm{d}}[m,n])\leq\tau_{2}\\\ \vdots&\vdots\\\
K&\tau_{K-1}<-\log({\bm{d}}[m,n]).\end{cases}$ (15)

The model is then given by

$\mathcal{A}({\bm{d}})=\sum_{k=1}^{K}-{\bm{k}}_{A_{k},B_{k},\sigma_{1,k},\sigma_{2,k}}\ast(\bm{1}_{\bm{r}_{\bm{d}}=k}\odot
f_{\alpha_{k},\beta_{k}}({\bm{d}})),$ (16)

where each kernel is given by (13), the indicator function is given by

$\bm{1}_{\bm{r}=k}[m,n]=\begin{cases}0,&\bm{r}[m,n]\neq k;\\\
1,&\bm{r}[m,n]=k,\end{cases}$ (17)

and $\odot$ indicates elementwise multiplication. Put simply, the multikernel
model uses a combination of parametric kernel models, each acting in a
different region of the input (specified by the indicator). Fitting the
multikernel model, like the parametric kernel, involves solving a nonconvex
optimization problem over all the scalar parameters, for which we again used
the LBFGS algorithm [19].

###  3.4 Descattering algorithm

We have so far discussed scatter estimation from a known direct radiograph,
but our ultimate goal is descattering, estimating a direct radiograph from a
total radiograph containing scatter. With the scatter models from Section 3.3
in hand, we are ready to descatter by solving
${\bm{t}}=\mathcal{A}({\bm{d}})+{\bm{d}}$ for ${\bm{d}}$. (One could readily
incorporate regularization on ${\bm{d}}$ and terms modeling noise or other
errors in this model; we leave these for future work.) Inverting this equation
may be challenging because it involves additional underlying optimization
problems such as (10), (12), or (14) that can themselves be challenging and
nonconvex. Here, we use an iterative fixed-point method. Fixed-point
algorithms are commonly used for descattering, e.g., see [3]. The wrinkle here
is that, in the local learning case, the scatter model parameters themselves
change as a function of the current estimate of ${\bm{d}}$. The pseudocode for
our descattering algorithm is given in Algorithm 1.

1:function descatter(${\bm{t}}$,
$\\{({\bm{d}}_{t},{\bm{s}}_{t})\\}_{t=1}^{T}$)

2: initialize: ${\bm{d}}\leftarrow{\bm{t}}$,

3: for fixed number of iterations do

4: fit a model, $\mathcal{A}_{\bm{d}}$ $\triangleright$ e.g., via (10), (12),
or (14)

5: ${\bm{d}}\leftarrow{\bm{t}}-\mathcal{A}_{\bm{d}}({\bm{d}})$

6: ${\bm{d}}[{\bm{d}}<0]\leftarrow 0$ $\triangleright$ enforces nonnegativity
of ${\bm{d}}$

7: return ${\bm{d}}$

Algorithm 1 Fixed point descattering.

###  3.5 Monoenergetic and polyenergetic density reconstruction

After descattering, we have an estimate of the set of direct radiographs,
$\\{{\bm{d}}_{v}\\}_{v=1}^{V}$, wherein each ${\bm{d}}_{v}$ is estimated
independently. The next tasks is to reconstruct the density field of the
object, ${\bm{\rho}}$ from this set, which amounts to inverting (3) or (5)
(and, consequently, (2) or (4) and (1)) with respect to the underlying
density. In general, if the inversion is ill-posed (e.g., limited view
imaging), additional regularization can be used during reconstruction [20, 1].
Here, we consider the case when enough views are acquired to adequately
reconstruct the 3D object. In particular, because the objects we are imaging,
which are of interest in important nuclear security applications, are
spherically symmetric, a single view suffices for reconstruction [21]. This
inversion is a well-studied problem and not the main focus of this work; we
therefore leave the details of our approach to Sections 7.1 and 7.2.

##  4 Experiments and results

To evaluate our algorithms, we performed descattering and reconstruction
experiments with simulated data. These experiments elucidate the descattering
performance of the methods, their ability to fit to known scatter, and the
effect the number of local neighbors used. In this section, we describe
implementation details for each comparison method, our datasets, and our
experiments and results.

###  4.1 Preliminaries

Implementation details. In all experiments, scatter models were fit as
described in Section 3.3, and all code was written in Python. The area outside
the object, i.e., beyond a 5 cm radius, was excluded in all cost functionals
because, in practice, there would be no difficulty locating the outer edge of
the object. Full-resolution images were $257\times 257$, but scatter fitting
was performed on $4\times$ downsampled versions (to improve fitting speed and
because the scatter is spatially smooth; this is common in the literature,
e.g., [3]) and scatter estimates were upsampled to the original size using
bilinear interpolation. Data was rescaled to simplify the tuning of fitting
routines. Specifically: each training direct radiograph was divided by the
Frobenius norm of the entire direct training set (square root of the sum of
the squares of all pixels). Each training scatter radiograph was divided by
the Frobenius norm of the entire scatter training set. At testing time, model
inputs were divided by the previously-computed training direct norm and
outputs were multiplied by the previously computed training scatter norm. For
CG (for the convolutional model), the number of iterations was 40. For both
the parametric kernel and multikernel models, the number of LBFGS iterations
was 20 with step size 1.0 and back-tracking line search was enabled; the
initial parameters were
$A=1.0,B=1.0,\sigma_{1}=4,\sigma_{2}=64,\alpha=1.0,\beta=0$. For the
multikernel model, we used $K=3$. While we did not make a careful study of the
runtimes of the proposed methods, fitting global models was slower than
fitting local ones (minutes versus seconds on a consumer-level laptop) and
evaluating the models once fit took under a second for all models. Code for
our method is available upon request.

Inverse Abel transforms were performed with the three-point algorithm of [21]
from the PyAbel Python package. Polyenergetic reconstructions were slower than
monoenergetic ones due to the lookup table, but both took on the order of
seconds on a consumer-level laptop.

Box and whisker plots. Throughout this section, we summarize results using box
and whisker plots, e.g., see Figure 3. Unless otherwise noted, each box
summarizes the ten testing images. Whiskers indicate minimum and maximum
values and boxes indicate first and third quartile. The orange line in the
center of each box indicates the median value.

###  4.2 Datasets

Our experiments use two simulation datasets with properties motivated by our
nuclear security application: large, dense, spherical objects. We chose to use
objects comprising a single material (plus a known collimator) to simplify
reconstruction and focus on the effects of descattering.

We generated synthetic radiographs using Monte Carlo particle transport
simulation software (MCNP6 [15]) that consisted of five concentric shells. The
outermost shell had a fixed radius of 5 cm and the other four shell radii were
randomly varied between 0.25 and 2.5 cm. The material was set to uranium for
the entire object, but the density of each shell was randomly selected from a
set of five densities, and the mass of the object was conserved through a
multiplicative constant. The object center was placed 133 cm from the source
and the detector was placed at 525 cm from the source. A tantalum collimator
(3.5 inches thick, 6.4 inch diameter, 1.4 inch opening, and 2.55 degree half-
angle) was placed at 98 cm from the source. While the collimator design was
not modified for each object, it served to reduce the dynamic range of the
image, in general. We simulated radiographs of these objects with both a
monoenergetic and polyenergetic X-ray source. The monoenergetic single-
material dataset was generated using a 1.5 MeV X-ray source and the
polyenergetic single-material dataset was generated using a bremsstrahlung
with a 19.4 MeV endpoint energy. We generated synthetic radiographs using the
radiography tally in MCNP, which generates and image produced by X-rays. This
tally is a quasideterministic calculation that maps each particle to each
pixel and the final image does not contain stochastic noise. Radiographs of
direct and scattered radiation were generated separately. The direct radiation
calculations utilized 50000 particles while the scattered radiation
calculations utilized 1e6 particles. Additionally, the scattered radiation
calculations included the PDS card in MCNP, which estimates all possible
reaction contributions for each collision and significantly reduces the
appearance of fliers due to rare coherent scatter events. Degradation of the
image due to source blur and detector effects such as scintillator efficiency
and blur were not modeled in this study.

For the monoenergetic dataset, we generated 99 pairs of direct and scatter
radiographs. For the polyenergetic dataset, we generated 100 pairs. These were
the largest datasets we could simulate in a reasonable amount of (order of
days). In both cases, we partition the datasets into a testing set consisting
of ten pairs (used for evaluation of algorithms) and a training set consisting
of the remaining pairs (used for fitting models). These partitions were chosen
to balance the need for computing statistics on the testing set with the
desire to have as large a training set as possible.

###  4.3 Descattering and density reconstruction

$\displaystyle 0.00$$\displaystyle 0.25$$\displaystyle 0.50$$\displaystyle
0.75$$\displaystyle 1.00$median absolute density error [g/cm3]w/ scattersingle
fieldparametricmultikernelconvolutionDSEw/out scattermono dataset

$\displaystyle 0.00$$\displaystyle 0.25$$\displaystyle 0.50$$\displaystyle
0.75$$\displaystyle 1.00$median absolute density error [g/cm3]w/ scattersingle
fieldparametricmultikernelconvolutionDSEw/out scatterpoly dataset

Figure 3: Error in density reconstruction after descattering for several
scatter estimation methods, local and global versions. In each stacked pair of
boxes, the upper, white box is the global version and the lower, grey box is
the local version. For the sake of comparison, the top- and bottom-most boxes
give the error when scatter is not removed and when it is removed perfectly,
respectively. Local methods always outperformed their global counterparts.

(a) global convolution

(b) local convolution

Figure 4: Error magnitude in density reconstruction after descattering for one
representative object in the monoenergetic dataset.

Our main experiment evaluated the ability of each of the models to enable
accurate descattering and density reconstruction. The experimental setup
mimicked our actual application. We fit scatter models on a training set of
$({\bm{d}}_{t},{\bm{s}}_{t})$ pairs and then descattered a set of new total
radiograph using the fit models and Algorithm 1. Finally, we reconstructed the
underlying density based on the estimated direct radiograph. For local models,
we fixed the number of neighbors to three (we explore the effect of the number
of neighbors in Section 4.7).

We quantified the reconstruction errors in terms of the median absolute
density error (MADE) on the central slice of the 3D reconstruction,

$\operatorname*{median}_{m,n}\left(\lvert\hat{{\bm{\rho}}}_{t}[m,n,0]-{\bm{\rho}}_{t}[m,n,0]\rvert\right),$
(18)

in g/cm3, where $\hat{{\bm{\rho}}}$ is the reconstructed density,
${\bm{\rho}}_{t}$ is the ground truth, and we let $m$ and $n$ range over the
area where ${\bm{\rho}}_{t}[m,n,0]$, is nonzero. (Recall that the objects we
are imaging are 3D, hence the three indices for ${\bm{\rho}}$.) We use the
central slice it because matches one way we typically view the
reconstructions, i.e., as a slice; quantifying error on 1D profiles or the
whole volume would also be reasonable. We used the median rather than the mean
because the mean is dominated by errors at the edge of the simulation and at
the interface between materials.

Comparison methods. While our main focus is the comparison between global
fitting and the proposed local fitting, we note that the global versions of
the parametric model and multikernel model represent baseline comparison
methods from the recent literature [3]. As as an additional comparison, we
implemented the deep scatter estimation (DSE) network [13], which is a CNN-
based scatter estimation method (trained globally and incorporated into the
fixed point algorithm in the same manner as the other scatter models).

Results. The results for the descattering experiment are given in Figure 3.
The results show that local models uniformly outperform global ones, which
justifies the use of the proposed local fitting approach. For the
monoenergetic dataset, the best descattering was provided by the local
convolutional model (maximum MADE=0.040 g/cm3), followed by the local
parametric model (0.045 g/cm3) and the DSE (0.050 g/cm3). None of the methods
matched the performance of the reconstruction without scatter (maximum
MADE=.009 g/cm3). For the polyenergetic dataset, the best descattering was
provided by the DSE (maximum MADE=0.253 g/cm3), followed by the local
convolutional model (0.270 g/cm3). Here, the error when reconstructing from
the direct signal (0.243 g/cm3) is higher than for the monoenergetic case:
this is to be expected because the particle transport simulation used to
generate the data and the forward model used to reconstruct it are distinct.
Slight differences in the physical constants (beam spectrum and mass
attenuation coefficients) between the model and simulation cause these
reconstruction errors. We believe that this mismatch is a strength of the
experiments, because it reflects (at least partly) the uncertainties that
exist when reconstructing from experimental data. We show error images for the
convolutional method in Figure 4. This shows that errors are highest at
material interfaces, and that there exist both density under- and overshoots.
We provide reconstruction and error profiles in Figures 7 and 8.

Convergence. We quantified the convergence of the fixed point algorithm by
computing the normalized mean squared error (NMSE) at each iteration,
$\|{\bm{d}}+\mathcal{A}_{\bm{d}}({\bm{d}})-{\bm{t}}\|_{2}^{2}/\|{\bm{t}}\|_{2}^{2}.$
(Note that this is not identically zero because the ${\bm{d}}$s on the left
and right hand sides of Algorithm 1 line 5 represent two successive iterations
and are not in general equal.) Although the value of the NMSE and its rate of
decrease varied from image to image, it generally either plateaued within a
few iterations or dropped to very low values (below $10^{-10}$). In the latter
case, this suggests that the fixed point algorithm finds a ${\bm{d}}$ such
that ${\bm{d}}+\mathcal{A}_{\bm{d}}({\bm{d}})={\bm{t}}$; in the former case,
we suspect there is no solution to the system and we get
${\bm{d}}+\mathcal{A}_{\bm{d}}({\bm{d}})\approx{\bm{t}}$. This behavior
occurred for both local and global models. We provide convergence plots in
Figure 9.

Noise. In a separate, small experiment (described in Section 9) we explored
how descattering accuracy was affected by noise. Results showed that noise
passed without amplification through the descattering process, indicating that
performing descattering does not necessitate additional noise removal beyond
that performed in the reconstruction step (e.g., regularization).

###  4.4 Oracle scatter fitting

To explore the capacity of the models, we evaluated them in an oracle fitting
setting, i.e., we evaluated their ability to map from a single direct
radiograph, ${\bm{d}}$ to its corresponding scatter, ${\bm{s}}$, when both
${\bm{d}}$ and ${\bm{s}}$ are known (provided by some oracle that one would
not have in practice). A model’s performance at this task upper bounds its
descattering performance; however, it does not give a complete picture of a
model’s potential for descattering, since models may overfit the training data
and fail to generalize.

| method
---|---
dataset | w/out scatter | single field | parametric | multikernel | convolution | w/ scatter
mono | 0.017 | 0.017 | 0.168 | 0.061 | 0.192 | 1.268
poly | 0.268 | 0.268 | 0.358 | 0.371 | 0.364 | 1.281
Table 1: Maximum MADE (g/cm3) after oracle scatter fitting.

For each $({\bm{d}}_{t},{\bm{s}}_{t})$ in the monoenergetic and polyenergetic
testing sets, we fit a single field, parametric kernel, multi kernel, and
convolutional model. (This is equivalent to performing local fitting with a
single neighbor if the training set includes the testing set.) We used each
fitted model to make a prediction, $\hat{\bm{s}}_{t}$. We formed an estimated
direct transmission,
$\hat{{\bm{d}}}_{t}={\bm{d}}_{t}+{\bm{s}}_{t}-\hat{{\bm{s}}}_{t}$, and
reconstructed, yielding $\hat{\rho}_{t}$. We compared these estimated
densities to the ground truth in terms of median absolute density error (18).

The scatter fitting results (Table 1) show that all models have some potential
to improve density reconstruction via descattering, giving maximum density
errors of under 0.2 g/cm3 and 0.4 g/cm3 in the monoenergetic and polyenergetic
cases, respectively, as compared to errors of over 1.25 g/cm3 when scatter is
simply ignored. The single field model performs nearly identically to
reconstruction without scatter; this is expected, because each
$\hat{{\bm{s}}}_{t}$ and ${\bm{s}}_{t}$ are identical except for the fact
that, as described in Section 4.1, we downsample by a factor of four during
scatter fitting. We can view the parametric model as a generalization of the
convolutional model with extra flexibility in the tuning of the parameters of
the nonlinearity (11). That the parametric model gives a smaller error than
the convolutional model here but not in the main experiment in Section 4.3
suggests that this extra flexibility leads to some degree of overfitting,
i.e., less generalization to training data. The same observation holds for the
multikernel model.

###  4.5 Scatter estimation

We also assessed the methods’ scatter estimation performance without
descattering involved, i.e., fitting on a training set and predicting the
corresponding ${\bm{s}}$ for a given known ${\bm{d}}$ from the testing set. We
assessed the results in terms of the median absolute density error (18), when
reconstructing from ${\bm{t}}-\mathcal{A}({\bm{d}})$.

The results show that scatter estimation was not significantly easier than
estimation and descattering (never more than 5% lower error). We interpret
this to mean that the limiting factor for descattering performance is
accurately modeling scatter as opposed to solving for ${\bm{d}}$ from
${\bm{t}}$. We provide estimated scatter profiles in Figure 10 and example
local and global kernels in Figure 11.

###  4.6 Cases with increased scatter

To investigate the performance of local descattering as the amount of scatter
increases, we chose a single object from the polyenergetic dataset and
repeated the MCNP simulation with its mass scaled by a factor of 1.5 and 2
(leaving the collimator unchanged). These scaled objects have a much higher
scatter-to-direct ratio than the original set, with a maximum ratio (over
pixels) of 1.57 and 25.60, respectively, compared to 0.34 in the original
dataset.

(a) $1.5\times$ density

(b) $2\times$ density

Figure 5: Effect of increasing object density on reconstruction performance
using a local convolutional model and fixed-point descattering (Algorithm 1).
As the density scales up, the scatter-to-direct ratio increases, increasing
the error.

For each of these scaled simulations, we additionally perturbed the density by
adding or subtracting 1.0 or 2.0 g/cm3 from the density of each shell,
creating four variations (-2.0, -1.0, +1.0, +2.0). We fit a convolutional
model to these four variations and used it to descatter the scaled total
transmission radiograph. Thus, we can think of the four variations as forming
a local training set. We compare the original descattering result to the two
high-scatter cases in Figure 5. The first trend we note is that, as the
density increases, the reconstruction from the total radiograph (no
descattering) deviates more and more from the ground truth, indicating that
scatter has become a larger source of error. While descattering is successful
in the original simulation and in the $1.5\times$ case, it begins to struggle
in the $2\times$ case, reaching density errors of over 10 g/cm3 in some
regions. However, the local approach is still better than ignoring scatter
altogether.

###  4.7 Effect of number of neighbors on reconstruction performance

To evaluate the effect of the number of neighbors, we repeated the
descattering experiment with the local convolutional model with the number of
neighbors varying from one to five. The results (Figure 6) show that the
reconstruction error is robust with respect to the number of neighbors, with
local fitting consistently outperforming global no matter the number of
neighbors. Although we might expect that a small number of neighbors could
lead to overfitting and that increasing the number of neighbors may cause
underfitting, we do not see clear evidence of this here. (Except that the
global model, which is equivalent on our data to a local model with 90
neighbors, performs worse than the local ones.)

Figure 6: Effect of the neighborhood size on density reconstruction error
using the local convolutional model.

##  5 Discussion

Local models vs CNNs. Our experiments showed that, using the proposed local
fitting scheme, the performance of simple models can rival that of a CNN (the
DSE [13]); we now compare and contrast these methods. One advantage of the
proposed method is that it does not require training a CNN, which is slow and
requires challenging hyperparameter tuning. Another is its flexibility in
terms of training: in a setting where we know approximately what we want to
image, we could perform a few scatter simulations, fit a local convolutional
model, and then use it to descatter experimental data. This aspect is
especially appealing for applications wherein creating a training set with
scatter of even tens of objects takes an unacceptable amount of time, due to
the complexity of the object or desired simulation accuracy. Fitting a deep
CNN on the same few scatter simulations would likely lead to overfitting. It
is also typically much easier to simulate the direct radiographs than to
simulate scatter. So in practice, we could simulate a dense set of direct
radiographs for performing the local matching in our algorithm and simulate
scatter only as necessary. The proposed local learning approach affords this
flexibility. The local fitting performed adaptively at reconstruction time
also means that new examples can be added to the training dataset and be used
by the model without retraining. Finally, the proposed model is simple to
describe and easy to interpret, in the sense that the learned kernels can be
viewed as images and intuitively encode the spatial component of scatter. The
main advantage of using a CNN over the proposed approach is that it is faster
at testing time (milliseconds vs seconds to estimate a scatter radiograph).
This may be attractive in settings where data is collected quickly and scatter
estimation needs to be immediate, but may not matter in many practical
settings.

Taking a broader perspective, the local convolutional model is a piecewise
linear function from images to images, i.e., $\mathbb{R}^{N_{1}\times
N_{2}}\to\mathbb{R}^{N_{1}\times N_{2}}$, with the pieces being defined by the
neighborhoods in the training data. ReLU CNNs, including the DSE, are also
piecewise linear functions [22]. The two models differ only in the
parameterization of the linear pieces and their boundaries. These similarities
may partly explain the similarity in performance between the two models.

Limitations. Our simulations did not account for scatter within the detector,
which, depending on the hardware, may contribute significantly to measured
radiographs [11]. Our local fitting method may be able to account for this
type of scatter, but it would need to be included in the training data. Our
testing datasets, although not small by the standards of other descattering
work, are probably not large enough to detect small differences between
methods, especially given the large amount of variation in the data. Our
experimental datasets consisted only of spherically symmetric objects and the
performance of the scatter estimation may be different for different types of
objects. However, we note that spherical objects are themselves of practical
importance in nuclear security applications and, further, there is nothing
about the proposed scatter model that explicitly relies on this symmetry. See
Section 10. for a proof-of-concept experiment on an ellipsoidal object.
Finally, in cases where the direct-to-scatter mapping is one-to-many, no
technique can deterministically estimate scatter from the direct radiograph
perfectly, including the proposed method. However, the proposed technique may
still be practically useful; further, performance might be improved by finding
neighbors based on total transmission, which could be unique even in a case
where the direct signal is not.

##  6 Conclusions and future work

Our experiments show that, when fit locally, even simple models can describe
scatter well, thus avoiding the algorithmic challenges and computational
burden of fitting and using a nonlinear or shift-variant model. Of the four
models we compared, we recommend the local convolutional model, given that it
is easy to fit and provides excellent descattering performance. The proposed
approach offers performance on par with a recent deep learning-based method,
with the advantage that it avoids an expensive network training step.

We see several possible extensions of this work. In the context of X-ray CT,
our method can easily be used for multiview (i.e., not spherically symmetric)
tomography problems, including medical imaging. However, these settings may be
challenging if the variety in the objects to be imaged is large; this would
necessitate gathering a larger training set to ensure good neighbors can be
found. Validating the method in additional settings is an important next step.
We would also like to the study the effect of training set density on
reconstruction performance; ideally, we could find a rule of thumb for how
many simulations are needed to achieve a certain accuracy for a given
application. Another direction may be to study different distance metrics for
neighbor determination, or to learn a good metric from the training data [23].
Finally, we believe that the local modeling approach may have applications
outside of scatter estimation. Wherever we have a training set of image pairs,
we can potentially fit simple, instance-adaptive local models instead of
complex, global models.

##### Acknowledgments.

Erik Skau, Brendt Wohlberg, and Luke Hovey, LANL, provided useful discussions.
Alexander Sietsema, Michigan State University, worked on related experiments
and discussions.

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

##### Supplemental document.

See Supplement 1 for supporting content.

## References

  * [1] S. Ravishankar, J. C. Ye, and J. A. Fessler, “Image reconstruction: From sparsity to data-adaptive methods and machine learning,” Proceedings of the IEEE 108, 86–109 (2020).
  * [2] C. Cohen-Tannoudji, J. Dupont-Roc, and G. Grynberg, _Atom-Photon Interactions: Basic Processes and Applications_ (VCH PUBN, 1998).
  * [3] M. Sun and J. M. Star-Lack, “Improved scatter correction using adaptive scatter kernel superposition,” Physics in Medicine and Biology 55, 6695–6720 (2010).
  * [4] J. P. Stonestrom and A. Macovski, “Scatter considerations in fan beam computerized tomographic systems,” IEEE Transactions on Nuclear Science 23, 1453–1458 (1976).
  * [5] E.-P. Rührnschopf and K. Klingenbeck, “A general framework and review of scatter correction methods in X-ray cone-beam computerized tomography. part 1: Scatter compensation approaches,” Medical Physics 38, 4296–4311 (2011).
  * [6] E.-P. Rührnschopf and K. Klingenbeck, “A general framework and review of scatter correction methods in cone beam CT. part 2: Scatter estimation approaches,” Medical Physics 38, 5186–5199 (2011).
  * [7] J. H. Mason, A. Perelli, W. H. Nailon, and M. E. Davies, “Polyquant CT: direct electron and mass density reconstruction from a single polyenergetic source,” Physics in Medicine & Biology 62, 8739–8762 (2017).
  * [8] I. Elbakri and J. Fessler, “Statistical image reconstruction for polyenergetic X-ray computed tomography,” IEEE Transactions on Medical Imaging 21, 89–99 (2002).
  * [9] A. Kak and M. Slaney, _Principles of Computerized Tomographic Imaging_ , Classics in Applied Mathematics (Society for Industrial and Applied Mathematics, 2001).
  * [10] 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),” .
  * [11] N. Bhatia, D. Tisseur, and J. M. Létang, “Convolution-based scatter correction using kernels combining measurements and monte carlo simulations,” Journal of X-Ray Science and Technology 25, 613–628 (2017).
  * [12] B. Ohnesorge, T. Flohr, and K. Klingenbeck-Regn, “Efficient object scatter correction algorithm for third and fourth generation CT scanners,” European Radiology 9, 563–569 (1999).
  * [13] J. Maier, S. Sawall, M. Knaup, and M. Kachelrieß, “Deep scatter estimation (DSE): Accurate real-time scatter estimation for X-ray CT using a deep convolutional neural network,” Journal of Nondestructive Evaluation 37 (2018).
  * [14] D. Tisseur, N. Bhatia, N. Estre, L. Berge, D. Eck, and E. Payan, “Evaluation of a scattering correction method for high energy tomography,” EPJ Web of Conferences 170, 06006 (2018).
  * [15] 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,” Tech. Rep. LA-UR-18-20808, Los Alamos National Laboratory (2018).
  * [16] Geant4 Collaboration, “Geant4—a simulation toolkit,” Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 506, 250–303 (2003).
  * [17] W. S. Cleveland, “Robust locally weighted regression and smoothing scatterplots,” Journal of the American Statistical Association 74, 829–836 (1979).
  * [18] J. Shewchuk, “An introduction to the conjugate gradient method without the agonizing pain,” Tech. rep., Carnegie Mellon University (1994).
  * [19] D. C. Liu and J. Nocedal, “On the limited memory BFGS method for large scale optimization,” Mathematical Programming 45, 503–528 (1989).
  * [20] M. T. McCann and M. Unser, “Biomedical image reconstruction: From the foundations to deep neural networks,” Foundations and Trends in Signal Processing 13, 283–359 (2019).
  * [21] C. J. Dasch, “One-dimensional tomography: a comparison of abel, onion-peeling, and filtered backprojection methods,” Applied Optics 31, 1146 (1992).
  * [22] G. F. Montufar, R. Pascanu, K. Cho, and Y. Bengio, “On the number of linear regions of deep neural networks,” in _Advances in Neural Information Processing Systems,_ vol. 27 Z. Ghahramani, M. Welling, C. Cortes, N. Lawrence, and K. Q. Weinberger, eds. (Curran Associates, Inc., 2014), pp. 2924–2932.
  * [23] B. Kulis, “Metric learning: A survey,” Foundations and Trends in Machine Learning 5, 287–364 (2013).
  * [24] M. Beister, D. Kolditz, and W. A. Kalender, “Iterative reconstruction methods in X-ray CT,” Physica Medica 28, 94–108 (2012).

##  7 Reconstruction Algorithms

###  7.1 Monoenergetic Case

In the monoenergetic, single-material case, the areal density from a given
angle, ${\rho_{A}}_{v}$, can be easily recovered from the direct radiograph at
that angle, ${\bm{d}}_{v}$. Given a flat field image (an image with no
object), ${\bm{f}}[m,n]\approx C{q}_{\text{in}}(r_{m,n}),$ we have
$\frac{{\bm{d}}_{v}[m,n]}{{\bm{f}}[m,n]}=\exp(-\xi{\rho_{A}}_{v}(r_{m,n})).$
Taking the negative logarithm (pixelwise) and dividing by the mass attenuation
coefficient gives ${\rho_{A}}_{v}(r)$. Going from
$\\{{\rho_{A}}_{v}\\}_{v=1}^{V}$ to ${\bm{\rho}}$ is a standard tomography
problem that can be resolved by any of numerous direct [9], iterative [24], or
learning-based [20, 1] methods.

We also consider the monoenergetic, multimaterial case, where there are
nuisance objects: elements of the scene that we are not interested in
reconstructing. If these can be imaged without the object of interest, their
effect can simply be included in the flat field radiograph and divided out as
described above. If they cannot be imaged alone, but we know the distribution
of their mass, we can compute their linear attenuation map (mass attenuation
coefficient multiplied by areal density) and subtract it after taking the
$\log$. Barring these types of prior information, the reconstruction problem
is much more challenging and beyond the scope of this paper.

###  7.2 Polyenergetic Case

In the single-material, polyenergetic case, we propose to follow the same two-
step inversion as in the monoenergetic case; however this is challenging
because we cannot simply divide to remove the effect of the incident field
${q}_{\text{in}}$ (due to its dependence on energy) and we cannot use a
logarithm to invert the sum of exponentials that appears in (5). To deal with
${q}_{\text{in}}$, we make the simplifying assumption that it is separable
into an (unknown) space-dependant and a (known) energy-dependent term as
${q}_{\text{in}}(r,E)={q}_{\text{in}}(r){q}_{E}(E)$ with
$\sum_{e}{q}_{E}(E_{e})=1$. Under this assumption, division by the flat field
gives

$\frac{{\bm{d}}_{v}[m,n]}{{\bm{f}}[m,n]}=\sum_{e}{q}_{E}(E_{e})\exp\left(-\xi(E_{e}){\rho_{A}}_{v}(r_{m,n})\right).$
(19)

Our solution relies on the fact that (19) is a monotone, scalar function of
${\rho_{A}}_{v}$ and therefore invertible via a lookup table. Our density
reconstruction algorithm is as follows. We consider (19) as a scalar function
of areal mass

$g({\rho_{A}})=\sum_{e}{q}_{E}(E_{e})\exp\left(-\xi(E_{e}){\rho_{A}}\right).$
(20)

We use this expression to compute $g({\rho_{A}}_{0}),g({\rho_{A}}_{1}),\dots
g({\rho_{A}}_{L-1})$, for a set of $L$ areal masses spaced evenly between 0
and some ${\rho_{A}}_{L-1}={\rho_{A}}_{\text{max}}$, the largest areal mass we
expect to encounter (if necessary, this can be determined iteratively). For
each pixel of the direct transmission, we can then perform a lookup in $g$ to
find the corresponding ${\rho_{A}}$ Once the $\\{{\rho_{A}}_{v}\\}$ are
recovered, reconstructing $\rho$ is a standard tomography problem, exactly
like in the monoenergetic case.

In the polyenergetic, multimaterial case, we cannot treat the nuisance objects
as part of the flat field and divide them out, nor can we subtract off a known
linear attenuation. For a single nuisance object with known density field,
e.g., a collimator,

$g({\rho_{A}},m,n)=\\\
\sum_{e}{q}_{E}(E_{e})\exp\left(-\xi(E_{e}){\rho_{A}}-\xi_{\eta}(E_{e}){\rho_{A}}_{\eta}(r_{m,n})\right),$
(21)

where $\xi_{\eta}(E_{e})$ and ${\rho_{A}}_{\eta}$ denote the known mass
attenuation coefficient and areal mass of the nuisance object. Inverting (21)
requires a three dimensional lookup table, which we can interpret as a version
of the lookup table for (20) for each location, $(m,n)$. The same approach can
be used for multiple nuisance objects by adding more terms to the exponential
in (21).

##  8 Additional Results

This section contains additional results from the experiments described in the
main text, Sections 4.3 and 4.5. Please see these sections for further
description.

$\displaystyle{0}$$\displaystyle{2}$$\displaystyle{4}$x
[cm]$\displaystyle{0}$$\displaystyle{5}$$\displaystyle{10}$$\displaystyle{15}$$\displaystyle{20}$$\displaystyle{25}$density
[g/cc]

(a) no descattering

$\displaystyle{0}$$\displaystyle{2}$$\displaystyle{4}$x
[cm]$\displaystyle{0}$$\displaystyle{5}$$\displaystyle{10}$$\displaystyle{15}$$\displaystyle{20}$$\displaystyle{25}$density
[g/cc]

(b) single field

$\displaystyle{0}$$\displaystyle{2}$$\displaystyle{4}$x
[cm]$\displaystyle{0}$$\displaystyle{5}$$\displaystyle{10}$$\displaystyle{15}$$\displaystyle{20}$$\displaystyle{25}$density
[g/cc]

(c) parametric kernel

$\displaystyle{0}$$\displaystyle{2}$$\displaystyle{4}$x
[cm]$\displaystyle{0}$$\displaystyle{5}$$\displaystyle{10}$$\displaystyle{15}$$\displaystyle{20}$$\displaystyle{25}$density
[g/cc]

(d) multiple kernels

$\displaystyle{0}$$\displaystyle{2}$$\displaystyle{4}$x
[cm]$\displaystyle{0}$$\displaystyle{5}$$\displaystyle{10}$$\displaystyle{15}$$\displaystyle{20}$$\displaystyle{25}$density
[g/cc]

(e) convolutional

$\displaystyle{0}$$\displaystyle{2}$$\displaystyle{4}$x
[cm]$\displaystyle{0}$$\displaystyle{5}$$\displaystyle{10}$$\displaystyle{15}$$\displaystyle{20}$$\displaystyle{25}$density
[g/cc]

(f) DSE

Figure 7: Example reconstruction profiles (monoenergetic dataset). Solid:
global model; dotted: local model. See Figure 8 for error plots where
differences are more apparent.

$\displaystyle{0.0}$$\displaystyle{2.5}$$\displaystyle{5.0}$x
[cm]$\displaystyle{-1.5}$$\displaystyle{-1.0}$$\displaystyle{-0.5}$$\displaystyle{0.0}$$\displaystyle{0.5}$$\displaystyle{1.0}$$\displaystyle{1.5}$density
error [g/cm3]

(a) single field

$\displaystyle{0.0}$$\displaystyle{2.5}$$\displaystyle{5.0}$x
[cm]$\displaystyle{-1.5}$$\displaystyle{-1.0}$$\displaystyle{-0.5}$$\displaystyle{0.0}$$\displaystyle{0.5}$$\displaystyle{1.0}$$\displaystyle{1.5}$density
error [g/cm3]

(b) parametric kernel

$\displaystyle{0.0}$$\displaystyle{2.5}$$\displaystyle{5.0}$x
[cm]$\displaystyle{-1.5}$$\displaystyle{-1.0}$$\displaystyle{-0.5}$$\displaystyle{0.0}$$\displaystyle{0.5}$$\displaystyle{1.0}$$\displaystyle{1.5}$density
error [g/cm3]

(c) multiple kernels

$\displaystyle{0.0}$$\displaystyle{2.5}$$\displaystyle{5.0}$x
[cm]$\displaystyle{-1.5}$$\displaystyle{-1.0}$$\displaystyle{-0.5}$$\displaystyle{0.0}$$\displaystyle{0.5}$$\displaystyle{1.0}$$\displaystyle{1.5}$density
error [g/cm3]

(d) convolutional

$\displaystyle{0.0}$$\displaystyle{2.5}$$\displaystyle{5.0}$x
[cm]$\displaystyle{-1.5}$$\displaystyle{-1.0}$$\displaystyle{-0.5}$$\displaystyle{0.0}$$\displaystyle{0.5}$$\displaystyle{1.0}$$\displaystyle{1.5}$density
error [g/cm3]

(e) DSE

Figure 8: Example reconstruction error profiles (monoenergetic dataset).
Solid: global model; dotted: local model. Descattering improves density
reconstruction across the entire profile and local models generally outperform
their global counterparts.

$\displaystyle 5$$\displaystyle 10$iteration$\displaystyle
10^{-22}$$\displaystyle 10^{-18}$$\displaystyle 10^{-14}$$\displaystyle
10^{-10}$$\displaystyle 10^{-6}$NMSE [1]global free kernel

$\displaystyle 5$$\displaystyle 10$iteration$\displaystyle
10^{-22}$$\displaystyle 10^{-18}$$\displaystyle 10^{-14}$$\displaystyle
10^{-10}$$\displaystyle 10^{-6}$NMSE [1]local free kernel

Figure 9: Convergence of the fixed point algorithm using the convolutional
model on the monoenergetic dataset. Each line represents one image in the test
set, all ten are not visible due to overlaps.

$\displaystyle-5$$\displaystyle 0$$\displaystyle 5$$\displaystyle m$
(cm)$\displaystyle 1.0$$\displaystyle 1.5$$\displaystyle 2.0$$\displaystyle
2.5$$\displaystyle\bm{s}[m,0]$$\displaystyle\times 10^{-7}$

(a) single field

$\displaystyle-5$$\displaystyle 0$$\displaystyle 5$$\displaystyle m$ (cm)

(b) parametric kernel

truthgloballocal

$\displaystyle-5$$\displaystyle 0$$\displaystyle 5$$\displaystyle m$ (cm)

(c) multiple kernels

$\displaystyle-5$$\displaystyle 0$$\displaystyle 5$$\displaystyle m$ (cm)

(d) convolutional

Figure 10: Example scatter estimations from direct on the monoenergetic
dataset, where fitting was performed on a separate training set.

(a) global kernel

(b) local kernel

(c) kernel profiles

Figure 11: Qualitative comparison of the global kernel versus the local kernel
learned for a single profile, plotted here as a function of position in
pixels. The local kernels tend to be sharper.

##  9 Descattering in the presence of noise

In this small experiment, we explored whether the proposed descattering
approach magnifies noise. We added various levels of AWGN to the total
transmission and performed descattering. Our results (Figure 12) show that
noise passes though descattering without significant amplification, i.e., the
RMSE of the descattered direct radiograph w.r.t the ground truth direct
radiograph increases linearly with the standard deviation of the noise.

We also verified that, for the levels of noise in the above plot, the nearest
neighbors of the total transmission radiographs are not changing. This
explains the linear behavior with respect to the noise, because when the
neighbors don’t change (with respect to same algorithm in noiseless case), the
scatter model is a linear function.

$\displaystyle{0}$$\displaystyle{5}$$\displaystyle{10}$$\displaystyle{15}$$\displaystyle{20}$$\displaystyle{25}$$\displaystyle{30}$noise
stdev [% of T
stdev]$\displaystyle{0}$$\displaystyle{2}$$\displaystyle{4}$$\displaystyle{6}$$\displaystyle{8}$RMSE
of $\displaystyle D_{hat}$$\displaystyle\times{10^{-7}}{}$ Figure 12: Noise
after descattering as a function of noise in the total transmission. The
relationship is linear, suggesting that noise passes without magnification
through the descattering procedure.

##  10 Nonspherical objects

As a small proof-of-concept, we performed scatter estimation on a simulated
ellipsoidal object using a set of four other ellipsoidal objects as training.
Results from the proposed local convolutional model are given in Figure 13.
While preliminary, these results suggest that the proposed method can work on
nonspherical objects.

$\displaystyle{0}$$\displaystyle{250}$$\displaystyle{500}$$\displaystyle{750}$$\displaystyle{0}$$\displaystyle{200}$$\displaystyle{400}$$\displaystyle{600}$$\displaystyle{800}$$\displaystyle{0}$$\displaystyle{5}$$\displaystyle{10}$$\displaystyle{15}$density
[g/cm3]

(a) Object density ($xy$-slice)

$\displaystyle{0}$$\displaystyle{250}$$\displaystyle{500}$$\displaystyle{750}$$\displaystyle{0}$$\displaystyle{200}$$\displaystyle{400}$$\displaystyle{600}$$\displaystyle{800}$$\displaystyle{2}$$\displaystyle{3}$$\displaystyle{4}$$\displaystyle{5}$percent
error

(b) Error in scatter estimate

$\displaystyle{0}$$\displaystyle{250}$$\displaystyle{500}$$\displaystyle{750}$$\displaystyle{0}$$\displaystyle{200}$$\displaystyle{400}$$\displaystyle{600}$$\displaystyle{800}$$\displaystyle{0}$$\displaystyle{250}$$\displaystyle{500}$$\displaystyle{750}$$\displaystyle{0}$$\displaystyle{200}$$\displaystyle{400}$$\displaystyle{600}$$\displaystyle{800}$$\displaystyle{1.65}$$\displaystyle{1.70}$$\displaystyle{1.75}$$\displaystyle{1.80}$$\displaystyle{1.85}$scatter$\displaystyle\times{10^{-6}}{}$
(c) True (left) and predicted (right) scatter

Figure 13: Proof-of-concept experiment on a nonspherical object.