TRINIDI: Time-of-Flight Resonance Imaging with Neutrons for Isotopic Density Inference

TL;DR

TRINIDI is a novel algorithm for neutron TOF resonance imaging that accurately reconstructs 2D and 3D isotope densities by addressing complex confounding effects through a two-step inverse modeling process.

Contribution

It introduces a two-step inversion method that accounts for non-linear absorption, noise, and background variations in neutron resonance imaging.

Findings

Accurately reconstructs 2D isotope density images.

Successfully extends to 3D volumetric density reconstruction.

Demonstrates effectiveness on both synthetic and real data.

Abstract

Accurate reconstruction of 2D and 3D isotope densities is a desired capability with great potential impact in applications such as evaluation and development of next-generation nuclear fuels. Neutron time-of-flight (TOF) resonance imaging offers a potential approach by exploiting the characteristic neutron absorption spectra of each isotope. However, it is a major challenge to compute quantitatively accurate images due to a variety of confounding effects such as severe Poisson noise, background scatter, beam non-uniformity, absorption non-linearity, and extended source pulse duration. We present the TRINIDI algorithm which is based on a two-step process in which we first estimate the neutron flux and background counts, and then reconstruct the areal densities of each isotope and pixel. Both components are based on the inversion of a forward model that accounts for the highly non-linear absorption, energy-dependent emission profile, and Poisson noise, while also modeling the substantial spatio-temporal variation of the background and flux. To do this, we formulate the non-linear inverse problem as two optimization problems that are solved in sequence. We demonstrate on both synthetic and measured data that TRINIDI can reconstruct quantitatively accurate 2D views of isotopic areal density that can then be reconstructed into quantitatively accurate 3D volumes of isotopic volumetric density.

Tables

Table 1. Table I: Ground truth and estimated values of nuisance parameters for the simulated data sample.

Quantity	Ground Truth
$\left[\begin{array}{c}\mathrm{isotopes}\\ {𝐳}^{\top }{}^{\left[\mathrm{mmol}/{\mathrm{cm}}^2\right]}\end{array}\right]$
${\alpha }_1,{\alpha }_2$	$0.483,0.685$
${𝛉}^{\top }$
Quantity	Estimates
$\left[\begin{array}{c}\mathrm{isotopes}\\ {\widehat{𝐳}}^{\top }{}^{\left[\mathrm{mmol}/{\mathrm{cm}}^2\right]}\end{array}\right]$
${\hat{\alpha }}_1,{\hat{\alpha }}_2$	$0.481,0.687$
${\widehat{𝛉}}^{\top }$

Table 2. Table II: Mean and standard deviation of areal density maps of ground truth ( Z 𝑍 Z ), linear baseline reconstruction  [ 37 ] ( Z lin superscript 𝑍 lin Z^{\mathrm{lin}} ), and TRINIDI reconstruction ( Z ^ ^ 𝑍 \hat{Z} ) in units of [ mmol / cm 2 ] delimited-[] mmol superscript cm 2 [\mathrm{mmol/cm^{2}}] of the simulated data sample.

Isotopes	Ground Truth	Linear Baseline [37]	TRINIDI
${}_{}{}^{238}U$	$5.0\pm 0.0$	$0.825\pm 0.082$	$5.065\pm 1.045$
${}_{}{}^{239}\mathrm{Pu}$	$3.0\pm 0.0$	$3.613\pm 0.336$	$3.009\pm 0.341$
${}_{}{}^{240}\mathrm{Pu}$	$0.2\pm 0.0$	$0.062\pm 0.001$	$0.202\pm 0.021$
${}_{}{}^{181}\mathrm{Ta}$	$4.0\pm 0.0$	$0.631\pm 0.021$	$4.015\pm 0.376$
${}_{}{}^{241}\mathrm{Am}$	$0.5\pm 0.0$	$0.703\pm 0.067$	$0.499\pm 0.039$

Table 3. Table III: Nuisance parameter values for the TaW TaW \mathrm{Ta}\mathrm{W} sample.

Quantity	Estimates
$\left[\begin{array}{c}\mathrm{isotopes}\\ {\widehat{𝐳}}^{\top }{}^{\left[\mathrm{mmol}/{\mathrm{cm}}^2\right]}\end{array}\right]$
${\alpha }_1,{\alpha }_2$	$0.483,0.685$
${𝛉}^{\top }$

Table 4. Table IV: Mean and standard deviation of areal densities and thicknesses vs. ground truth thicknesses of the TaW TaW \mathrm{Ta}\mathrm{W} sample.

${}_{\text{Quantity}}{}^{}\setminus {}^{\text{Isotopes}}$	${}_{}{}^{181}\mathrm{Ta}$	$W$
Areal Density Estimate ${}^{\left[\frac{\mathrm{mmol}}{{\mathrm{cm}}^2}\right]}$	$20.56\pm 1.59$	$21.11\pm 2.09$
Ground Truth Thickness [cm]	$0.242$	$0.175$
Estimated Thickness [cm]	$0.223\pm 0.017$	$0.197\pm 0.020$

Table 5. Table V: Nuisance parameter values for the fuel pellet sample.

Quantity	Estimates
$\left[\begin{array}{c}\mathrm{isotopes}\\ {\widehat{𝐳}}^{\top }{}^{\left[\mathrm{mmol}/{\mathrm{cm}}^2\right]}\end{array}\right]$
${\hat{\alpha }}_1,{\hat{\alpha }}_2$	$2.39,1.02$
${\widehat{𝛉}}^{\top }$

Table 6. Table VI: Average values of volumetric densities, 𝐱 ^ ^ 𝐱 \hat{\mathbf{x}} , of the fuel pellet sample in units of [ mmol / cm 3 ] delimited-[] mmol superscript cm 3 {[\mathrm{mmol/cm^{3}}]} .

Quantity	Estimates
$\left[\begin{array}{c}\mathrm{isotopes}\\ {\widehat{𝐱}}^{\top }\end{array}\right]$

Table 7. Table VII: Mass density estimates the fuel pellet sample compared to mass spectrometry  [ 36 ] (ground truth) in units of g / cm 3 g superscript cm 3 \mathrm{g/cm^{3}} .

${}_{\text{Mass Density}}{}^{}\setminus {}^{\text{Isotopes}}$	${}_{}{}^{237}\mathrm{Np}$	${}_{}{}^{238}U$	${}_{}{}^{239}\mathrm{Pu}$	${}_{}{}^{240}\mathrm{Pu}$	${}_{}{}^{241}\mathrm{Am}$
Ground Truth	0.34	9.01	2.34	0.37	0.32
Estimated	0.331	9.014	2.000	0.435	0.340

Equations

$$E=\frac{1}{2}m{v}^2=\frac{1}{2}m{\left(\frac{L}{t}\right)}^2,$$

$$z=\frac{\ell \rho }{M}.[mol/c{\mathrm{m}}^2]$$

$${\mathbf{q}}_j=\exp \left\{-z{\mathbf{d}}_j\right\}.$$

$$x=\frac{\rho }{M}.[\mathrm{mol}/{\mathrm{cm}}^3]$$

$${\bar{Y}}_{\mathrm{s}}={\alpha }_1\left[\Phi \odot \exp \left\{-ZD\right\}+{\alpha }_{2}B\right],$$

$$R_{i,j}=\mathbb{P}\left\{\text{mbox}TOA={\mathbf{t}}_j^{\mathrm{A}}|\text{mbox}TOF={\mathbf{t}}_i^{\mathrm{F}}\right\}.$$

$${\mathbf{t}}_j^{\mathrm{A}}$$

$${\mathbf{t}}_i^{\mathrm{F}}$$

$${\bar{Y}}_{\mathrm{s}}$$

$${\bar{Y}}_{\mathrm{o}}$$

$$\Phi$$

$$B$$

$$\mathbf{b}(\text{uptheta}{)}^{\top }=\exp \{{\text{uptheta}}^{\top }P\},$$

$$1^{\top }\mathbf{v}=N_{\mathrm{p}},$$

$$(\forall i\in {\Omega }_{\mathrm{z}})Z_{i,\ast }$$

$$(\forall i\in {\Omega }_0)Z_{i,\ast }$$

$$\displaystyle Y_{\mathrm{s}}\sim\mathrm{Poisson}(\bar{Y}_{\mathrm{s}})\$$

$$\displaystyle Y_{\mathrm{o}}\sim\mathrm{Poisson}(\bar{Y}_{\mathrm{o}})\$$

$$\hat{\mathbf{v}}=N_{\mathrm{p}}\frac{Y_{\mathrm{o}}1}{1^{\top }{Y}_{\mathrm{o}}1}.$$

$$\mathbf{q}(\mathbf{z}{)}^{\top }=\exp \left\{-{\mathbf{z}}^{\top }D\right\}R,$$

$$\mathbf{q}(0{)}^{\top }=1^{\top },$$

$$({\text{upomega}}_0{)}_i$$

$$({\text{upomega}}_{\mathrm{z}}{)}_i$$

$${\bar{\mathbf{y}}}_{\mathrm{o}}^{\top }$$

$${\bar{\mathbf{y}}}_{\mathrm{sz}}^{\top }$$

$${\bar{\mathbf{y}}}_{s0}^{\top }$$

$${\bar{\mathbf{y}}}_{\mathrm{o}}^{\top }$$

$$=\frac{1^{\top }\mathbf{v}}{1^{\top }\mathbf{v}}{\left(\text{upvarphi}+\mathbf{b}(\text{uptheta})\right)}^{\top }={\left(\text{upvarphi}+\mathbf{b}(\text{uptheta})\right)}^{\top }.$$

$$\text{upvarphi}={\bar{\mathbf{y}}}_{\mathrm{o}}-\mathbf{b}(\text{uptheta}).$$

$${\bar{\mathbf{y}}}_{\mathrm{sz}}^{\top }$$

$$=\frac{{\alpha }_1}{{\text{upomega}}_{\mathrm{z}}^{\top }\mathbf{v}}{\text{upomega}}_{\mathrm{z}}^{\top }\left[\Phi \odot \left(\exp \left\{-1{\mathbf{z}}^{\top }D\right\}R\right)+{\alpha }_{2}B\right]$$

$$=\frac{{\alpha }_1}{{\text{upomega}}_{\mathrm{z}}^{\top }\mathbf{v}}{\text{upomega}}_{\mathrm{z}}^{\top }\left[(\mathbf{v}{\text{upvarphi}}^{\top })\odot \left(1\mathbf{q}(\mathbf{z}{)}^{\top }\right)+{\alpha }_2\mathbf{vb}(\text{uptheta}{)}^{\top }\right]$$

$$={\alpha }_1{\left[\text{upvarphi}\odot \mathbf{q}(\mathbf{z})+{\alpha }_2\mathbf{b}(\text{uptheta})\right]}^{\top }$$

$$={\alpha }_1{\left[({\bar{\mathbf{y}}}_{\mathrm{o}}-\mathbf{b}(\text{uptheta}))\odot \mathbf{q}(\mathbf{z})+{\alpha }_2\mathbf{b}(\text{uptheta})\right]}^{\top }.$$

$${\bar{\mathbf{y}}}_{s0}^{\top }$$

$$={\alpha }_1{\left[{\bar{\mathbf{y}}}_{\mathrm{o}}+({\alpha }_2-1)\mathbf{b}(\text{uptheta})\right]}^{\top }.$$

$$(\hat{\mathbf{z}},{\hat{\alpha }}_1,{\hat{\alpha }}_2,\widehat{\text{uptheta}})=\underset{\mathbf{z},{\alpha }_1,{\alpha }_2,\text{uptheta}}{\mathrm{argmin}}\{$$

Full text

We present the TRINIDI algorithm which is based on a two-step process in which we first estimate the neutron flux and background counts, and then reconstruct the areal densities of each isotope and pixel. Both components are based on the inversion of a forward model that accounts for the highly non-linear absorption, energy-dependent emission profile, and Poisson noise, while also modeling the substantial spatio-temporal variation of the background and flux. To do this, we formulate the non-linear inverse problem as two optimization problems that are solved in sequence. We demonstrate on both synthetic and measured data that TRINIDI can reconstruct quantitatively accurate 2D views of isotopic areal density that can then be reconstructed into quantitatively accurate 3D volumes of isotopic volumetric density.

I Introduction

Neutrons, unlike X-rays, protons, or electrons, primarily interact with the nucleus rather than the electron shell of atoms, and thus neutron radiography can offer complementary contrast mechanisms to more conventional radiography probes [1]. While X-rays generally get attenuated more by heavier isotopes, this simple relationship does not hold for neutron attenuation. Neutron imaging thus enables the non-destructive evaluation of nuclear fuel, bulk objects, and other materials formed of heavy nuclei [2], while revealing (X-ray-transparent) features such as plant roots and soil water [3, 4]. Applications range from direct methods such as transmission radiography and tomography as well as indirect methods that leverage neutron diffraction and small angle scattering to determine crystal structures, shapes of biological molecules as well as atomic resolution techniques such as neutron holography [5, 2, 6]. In addition, simultaneous X-ray and neutron tomography can provide multi-modal attenuation mapping such as the mechanical degradation (via X-rays) as well as chemical diffusion process (via neutrons) in lithium batteries [7, 8].

Neutron time-of-flight (TOF) imaging is a promising imaging technique that exploits the property that neutrons travel with a velocity that depends on their energy. At advanced neutron source facilities, such as ISIS [9], LANSCE [10], SNS [11], and J-PARC [12] spallation neutrons are pulsed and moderated down to yield a broad spectrum of energies. These types of neutron sources allow for neutron TOF measurements, where neutron energies are determined by the time it takes to traverse a known distance between the source and the detector. Thus, pulsed neutron TOF measurements are inherently hyperspectral, and can be used to measure neutron absorption or diffraction as a function of energy [13, 14, 15, 16] facilitating, for example, 2D temperature mapping [17, 18] through Doppler broadening of neutron absorption resonances as well as neutron diffraction imaging techniques for identification of crystal structure and material strains [19, 20].

This paper deals with the problem of reconstructing isotopic densities from neutron TOF transmission measurements given the knowledge of isotope-specific neutron cross section spectra. The neutron energies suitable for this modality range from approximately $1\,\mathrm{eV}$ to $10\,\mathrm{keV}$, categorized as resonance neutrons [16]. Absorption resonances are isotope specific and increase the attenuation of neutrons by orders of magnitude within a narrow energy range [21]. Quantitatively accurate isotope reconstruction is very challenging for a number of reasons. TOF neutron imaging measurements are typically exceptionally noisy, with average neutron counts as low as only $2$ counts per measurement bin after several hours of measuring. This is mainly because each measurement bin only corresponds to a single small pixel ($\lesssim 100\,\upmu$m) and a short TOF window ($\lesssim 100\,\mathrm{ns}$) and because even the most advanced neutron sources emit several orders of magnitude fewer particles than typical X-rays sources [22, 23]. In addition, the measurements contain substantial background counts which are typically about as strong as the signal itself [24]. Finally, the measurement physics are highly non-linear due to the non-linear absorption, the large dynamic range of the spectral neutron responses of the nuclei, and the effects of the extended pulse duration of the neutron beam, an energy-dependent phenomenon that is caused by the moderation process[16].

The usual method for low-noise measurements [25, 26, 27, 28] is to first compute a transmission spectrum (i.e. open beam normalized transmission measurement) from the TOF data and then to process it with software packages such as SAMMY [29], CONRAD [30] or REFIT [31]. However, these software packages require relatively low-noise estimates of the transmission spectrum as input. In addition, since these software packages are primarily designed for the analysis of individual neutron spectra, they would be very slow when applied to a large array of pixels.

In an imaging application low-noise spectra can be obtained by selecting and averaging regions of interest which can then be processed as a 1D signal [32, 24]. However, the averaging does not yield density estimates on a per-pixel basis, so the problem is shifted to the accurate estimation of transmission spectra at each pixel. It is due to this difficulty that the vast majority of papers on neutron resonance imaging is limited to 2D mapping of isotopes rather than quantification of isotope densities. Our work here aims to overcome these shortcomings.

The transmission spectrum is the expected ratio of the sample and open beam scans, both of which must be corrected for the severe background and beam non-uniformity. There are a number of reasons for the difficulty of accurately estimating these transmission spectra. Although explicit models of the background have been proposed [26], there are no widely accepted methods for jointly estimating background and beam non-uniformity. Moreover, these problems are exacerbated by the need to separately estimate the transmission at each pixel. This means that the spatio-temporal non-uniformity in the background and beam-profile must be accounted for along with the very high level of neutron counting noise. Importantly, since the transmission spectrum estimate is a ratio of two very noisy signals, the resulting spectrum is even noisier.

Notably among the early efforts of estimating 2D and 3D isotope densities are the work of Sato, Kamiyama, et al. [17, 33]. Here a large $8\!\times\!8$ pixel-coarse imaging detector is used to yield low-noise measurements. In addition, the effects of background events were neglected. In [34] a similarly sized detector was used where the estimation was based on the assumption that the isotope densities are simply proportional to the area of the transmission resonance dips. Losko, Vogel et al. [35, 36] accomplished density estimates with a pixel-dense imaging detector. Similar to the previous cases, the transmission spectra are relatively low noise likely due to a combination of much larger TOF bins or spatial denoising pre-processing. Such measures limit spatial and temporal resolution while allowing for only restricted energy range to only about $10\,\mathrm{eV}$. Our work here focuses on low count, high noise measurements.

Previously we proposed a method [37] based on a linearized model and the assumption of opaque “black resonances” [38] to estimate the background counts. While this research demonstrates fast and semi-quantitative results, it has a number of limitations. First, black resonances are not always available or entirely opaque after the blur caused by the extended pulse shape (see Appendix A). Second, the linearized model may not be sufficiently accurate as it neglects to model the extended pulse duration and model the counting noise adequately. Third, this approach implicitly still first estimates the very noisy transmission spectra and as a consequence densities of highly attenuating isotopes can be strongly under-estimated. (See the ${}^{238}\mathrm{U}$ and ${}^{240}\mathrm{Pu}$ estimates in Table 2 of [37].)

In this paper, we present the Time-of-flight Resonance Imaging with Neutrons for Isotopic Density Inference (TRINIDI) algorithm. An open-source Python implementation [39] is available in the spirit of reproducible research. Our model-based iterative reconstruction (MBIR [40, 41, 42, 43]) algorithm is based on a two-step process in which we first estimate the neutron flux and background, and then estimate the areal density of each isotope. Importantly, the initial computation of transmission spectra is avoided by modeling the neutron measurement in their native counting domain. Both components are based on the inversion of a forward model that accounts for the highly non-linear absorption, extended pulse duration, and Poisson noise, while also modeling the very substantial space and energy varying background. To do this, we formulate the non-linear inverse problem as two optimization problems that are solved in sequence.

We perform numerical experiments on both simulated and measured data that demonstrate the ability of TRINIDI to reconstruct quantitatively accurate isotopic areal density from 2D data. Finally, we use a series of 2D views of a nuclear fuel sample that were experimentally measured at the LANSCE facility [10] to compute a quantitatively accurate 3D reconstruction of the fuel sample’s material composition.

II The Neutron Imaging System

Figure 1 illustrates the major components of a neutron beamline. Each pulse contains neutrons with a wide range of energies. The flux neutrons are (quasi-)exponentially attenuated as they pass through a sample, and finally, they are detected by a high speed imaging detector that collects approximately $2000$ equi-spaced frames at a rate of $30\,\mathrm{ns}$ per frame. The background events are also counted by the detector and are assumed to derive from a variety of scattering events in the flight path.

Neutrons, unlike photons, have non-zero mass, so their energy is a function of their velocity, which is in turn a function of their TOF. Consequently, each TOF corresponds to a distinct neutron energy, and the set of images represent a spectral decomposition of the detected neutrons. The neutron velocities are sufficiently low for relativistic effects to be ignored, so the kinetic energy of the neutron is given by

[TABLE]

where $m\approx 1.045\!\times\!10^{-8}\,\mathrm{eV\,s^{2}/m^{2}}$ is the mass111We use $E=mc^{2}$ to express the unit of mass in terms of $\mathrm{eV}$., $v$ the velocity, and $t$ the TOF of a neutron and $L$ is the flight path length.222This means that $t\propto\nicefrac{{1}}{{\sqrt{E}}}$, and we use TOF and energy interchangeably

Figure 2 illustrates the propagation of a neutron pulse through an $L=10\,\mathrm{m}$ flight path. For this distance, a $1\,\mathrm{eV}$ neutron corresponds to a TOF of $720\,\mathrm{\upmu s}$, and a $100\,\mathrm{eV}$ neutron corresponds to a TOF of $72\,\mathrm{\upmu s}$. Notice that the higher energy neutrons travel faster, and a $100\times$ increase in energy results in a $10\times$ decrease in TOF.

Notice that the initial pulse is not instantaneous. This means that neutrons with a particular energy may arrive over a range of TOFs. Equivalently, the arrival time does not precisely determine the energy of the neutron. A wider pulse has the effect of blurring absorption features in the time domain and thus limiting the temporal resolution of the imaging system333similar to a pinhole camera where a wider aperture causes spatial blurring.. A longer flight path reduces the blurring444The absolute temporal blur is independent of the flight path length. However the longer the flight path the more the time-of-arrival is dominated by the velocity (energy) of the neutron as opposed to the uncertainty in the time the neutron enters the flight path. This effect reduces the relative blur observed at the detector plane., but also reduces the overall flux via the inverse-square law. There is therefore an inherent trade-off between temporal resolution and neutron flux. Lastly, notice that the duration of the pulse is energy-dependent555The energy-dependence of the pulse length comes from the moderation process, a mechanism to create a wide spectrum from a high-energy pulse. Lower energy neutrons are emitted with a longer pulse since there is a larger number of possible collisions in the moderator to arrive at that lower energy.. Accurate modeling of this effect will prove to be very important in solving our problem.

Our experiments will focus on neutron energies in the range of approximately $1\,\mathrm{eV}$ to $100\,\mathrm{eV}$, i.e. resonance region neutrons. At resonances, the interaction probability of an incident neutron can increase by orders of magnitude and cause strong, isotope specific attenuation peaks. Below $1\,\mathrm{eV}$ isotopes usually do not have neutron resonances, and thus do not provide sufficient isotope-specific characteristics. At energies above about $100\,\mathrm{eV}$ the neutron resonances are so dense that the inherent temporal blur due to the neutron pulse renders individual resonances unresolvable. Furthermore, at such high energies, the neutron counts are dominated by background, which makes the very noisy signal less useful.

Figure 3(a) shows a neutron spectrum as a function of the TOF. The blue line shows the number of counts averaged over many pixels of the measured spectrum through a sample. Notice that the spectrum shows a great deal of fine structure, i.e. the neutron resonances, that can be used to distinguish different isotopes. Since neutrons interact with the nuclei, these characteristic resonances can be used to distinguish between different isotopes of the same element, such as the fissile ${}^{235}\mathrm{U}$ and the non-fissile ${}^{238}\mathrm{U}$ uranium isotopes. The red line shows the background counts which are quite large, and can result from a wide range of scattering effects in the flight path and sample. Finally, the yellow line shows the open beam measurement that results from a combination of the flux and the background events. Importantly, the background is large and the flux has a great deal of variation, so both must be accurately modeled in order to produce an accurate reconstruction of a sample.

Figure 3(b) shows a plot of a spectral neutron measurement of a single pixel. Since the measurement consists of neutron counts, it is discrete, and it is clear that for individual pixels the signal-to-noise (SNR) ratio is very low. Neutron TOF imaging is challenging due to the very low SNR.

Finally, Figure 3(c) illustrates the cross section of various isotopes in the sample, in units of $\mathrm{cm^{2}/mol}$, as a function of neutron TOF. The concept of a neutron cross section is used to express the likelihood of interaction between an incident neutron and a target nucleus. As the cross section of the material increases, it more strongly attenuates the neutrons and reduces the transmission. For the purposes of this research, these cross section spectra are considered known666Cross section spectra are tabulated in databases such as ENDF/B-VIII.0[44], which are frequently updated through more accurate measurements. for particular isotopes that make up the sample.

As an example, imagine a sample made of a single isotope with a known cross section of $\mathbf{d}_{j}\,^{[\mathrm{cm^{2}/mol}]}$ at TOF bin $j$. Furthermore, assume that the sample is of uniform thickness $\ell\,^{[\mathrm{cm}]}$, with mass density $\rho\,^{[\mathrm{g/cm^{3}}]}$, and molar mass $M\,^{[\mathrm{g/mol}]}$, then it has a uniform areal density of

[TABLE]

Using Beer’s Law, we would expect the transmission of the neutrons through the material at TOF bin, $j$, to be

[TABLE]

Our first goal will be to estimate the areal density, $z$, for each pixel and each material in the sample. The areal density is predominantly desired with a quasi-2D, flat sample such as a foil or a physical slice of a larger sample.

Our second goal will be to estimate the volumetric density for each isotope in the sample. This is important when the sample is an arbitrarily shaped 3D object, such as the nuclear fuel sample used in our experiments. In this case, our goal is to reconstruct the 3D volumetric density, $x$, where

[TABLE]

We can estimate $x$ for each voxel and material in the sample by extracting the areal densities for each rotational view followed by tomographic reconstruction.

III Single-Radiograph Forward Model

In this section, we derive our forward model for the hyperspectral neutron measurements of a radiograph. Our goal will be to estimate a matrix $Z\in\mathbb{R}^{N_{\mathrm{p}}\times N_{\mathrm{m}}}$ where $N_{\mathrm{p}}$ is the number of scalar projections through a sample (pixels) and $N_{\mathrm{m}}$ is the number of unique isotopes that make up the sample. Then $Z_{i,m}$ represents the areal density in units of $\mathrm{mol}/\mathrm{cm}^{2}$ of the $i^{\mathrm{th}}$ projection through a sample for the $m^{\mathrm{th}}$ isotope.

The matrix $D\in\mathbb{R}^{N_{\mathrm{m}}\times N_{\mathrm{F}}}$ denotes the cross section dictionary in units of $\mathrm{cm}^{2}/\mathrm{mol}$ for $N_{\mathrm{F}}$ neutron energies. I.e. $D_{m,j}$ is the cross section for material $m$ at TOF $j$. Notice that each resonant isotope in the sample must be accounted for because distinct isotopes typically have distinct spectral signatures.

Assuming the conventional Beer-Lambert’s law attenuation model, the expected measured neutron counts would be

[TABLE]

where $\alpha_{1}$ and $\alpha_{2}$ are positive scalars, the symbol $\odot$ is the Hadamard product operator, $\varPhi_{i,j}$ is the neutron flux, and $B_{i,j}$ is the background for projection $i$ and TOF $j$, respectively.

The flux, $\varPhi$, and background, $B$, are important scan parameters that vary as a function of position and energy. In X-ray imaging, the background is relatively small and thus the flux can be accurately measured with an “open beam” calibration procedure. However, in neutron imaging, the background tends to be quite large, somewhat sample dependent, and difficult to measure separately. Consequently, it is also very difficult to measure the flux directly from a single open beam measurement. Therefore, a major challenge will be to formulate an imaging methodology that allows us to recover both the flux and background from available measurements.

From Section II, we know that the neutron TOF varies with energy. However, in practice, the measured time of arrival (TOA) does not exactly determine the theoretical time of flight (TOF) of the neutron due to uncertainty in the time the neutron is emitted. The neutrons are emitted in a causal pulse distribution with a long tail that has significant time duration as compared to the sampling rate of the detector camera. A neutron with a specific energy ($\leftrightarrow$ TOF) therefore has finite probability of being detected over a range of times ($\leftrightarrow$ TOA).

In order to model this effect, we introduce the resolution function or resolution operator, $R\in\mathbb{R}^{N_{\mathrm{F}}\times N_{\mathrm{A}}}$, representing the conditional probability

[TABLE]

The measurement bins are uniformly spaced in time, so we denote the vectors of TOF and TOA bins by $\mathbf{t}^{\mathrm{F}}\in\mathbb{R}^{N_{\mathrm{F}}}$ and $\mathbf{t}^{\mathrm{A}}\in\mathbb{R}^{N_{\mathrm{A}}}$, respectively, such that

[TABLE]

where $\Delta_{\mathrm{t}}$ is the TOF bin size, $t_{0}$ is the time of the first measured bin, and $i_{0}$ is a positive integer offset. Importantly, the set of TOF bins is taken to be sufficiently large that it includes any neutron that could possibly contribute to a measurement. This means that the set of TOA bins will be a contiguous subset of the set of TOF bins, and that $N_{\mathrm{F}}>N_{\mathrm{A}}$. Equivalently, $R$ is a sparse matrix that is taller than it is wide. Finally, we assume that the columns of $R$ sum to 1, so that $(\forall j)\sum_{i}R_{i,j}=1$. In Appendix A we show that this is because the uncertainty in TOA is modeled as a linear combination of temporal convolutions with the neutron pulse shape.

Using this notation, the final equation for the expected measured neutron flux is given by

[TABLE]

Note that both $D$ and $R$ are assumed to be known for the materials of interest and the particular neutron beamline. The unknowns that must be estimated are therefore $Z$, $\varPhi$, $B$, $\alpha_{1}$, and $\alpha_{2}$. It is clear that the dimension of the unknowns far exceeds twice the dimension of the measurements, $\bar{Y}_{\mathrm{s}}$. Consequently, the direct inversion of (5) is very ill-posed.

In order to make the inversion problem well posed, we make three additional assumptions. First, assume that we have made an additional open beam measurement, $\bar{Y}_{\mathrm{o}}$, with the sample removed. Second, assume that the flux and background can be modeled with a low rank approximation. And finally, assume there are two regions in our sample measurement, one in which the sample is not present, and another in which the sample has constant, non-zero areal density.

First, we will assume that a second open beam measurement is made with the sample removed. In this case, $Z=\mathbf{0}$, and we have that

[TABLE]

which uses the result (A) from Appendix A.

Second, in order to further improve the conditioning of the inversion problem, we will introduce the constraints that the neutron flux and background, $\varPhi$ and $B$, can be modeled with a rank-one approximation given by

[TABLE]

where $\mathbf{v}\in\mathbb{R}^{N_{\mathrm{p}}}$, $\boldsymbol{\upvarphi}\in\mathbb{R}^{N_{\mathrm{A}}}$, and

[TABLE]

where $\boldsymbol{\uptheta}\in\mathbb{R}^{N_{\mathrm{b}}}$ is a low dimensional parameter vector, and $P\in\mathbb{R}^{N_{\mathrm{b}}\times N_{\mathrm{A}}}$ is a matrix of basis functions used to model the spectral characteristics of the background counts. The specific form of $P$ is given in detail in Appendix B. In order to make the decomposition unique, we will also require that

[TABLE]

where $\mathbb{1}$ denotes a column vector of $1$’s, i.e. $\mathbf{v}$ is normalized so that its components have an average value of $1$.

Third, we will assume that in the sample measurement field of view there are two regions of pixels, denoted by $\Omega_{\mathrm{0}}$, the open beam region, and $\Omega_{\mathrm{z}}$, the uniformly dense region. In $\Omega_{\mathrm{0}}$ no part of the sample is covering the beam. In $\Omega_{\mathrm{z}}$ the sample has uniform cross section. More specifically, we assume

[TABLE]

where $\mathbf{z}\in\mathbb{R}^{N_{\mathrm{m}}}$ is an unknown vector of areal densities and $\mathbb{0}$ is a vector of zeros. These regions must be selected by the user based on knowledge of the sample properties, which is feasible in most applications. This can be achieved by leaving a part of the field of view unobscured to derive the $\Omega_{\mathrm{0}}$ region. One possible method to pick the $\Omega_{\mathrm{z}}$ region is to select a region in good faith, perform a material decomposition and then check whether the resulting density, $Z$, is actually approximately uniform in $\Omega_{\mathrm{z}}$. If the sample does not readily provide a region of approximately equal areal density, a calibration sample could be added to the field of view to derive the $\Omega_{\mathrm{z}}$ region.

Finally, the actual measurements $Y_{\mathrm{s}}$ and $Y_{\mathrm{o}}$ are assumed to be conditionally Poisson distributed, which is expressed as

[TABLE]

IV Estimation of Nuisance Parameters

In this section, we describe our procedure for estimating the nuisance parameters, $\mathbf{v}$, $\boldsymbol{\upvarphi}$, $\boldsymbol{\uptheta}$, $\alpha_{1}$, and $\alpha_{2}$. This is a challenging problem because it is equivalent to estimation of the background, which tends to be large and spatially variable in neutron imaging.

We start by estimating $\mathbf{v}$ by observing that columns of the open beam scan, $\bar{Y}_{\mathrm{o}}$, must all be proportional to $\mathbf{v}$, and that $\mathbb{1}^{\top}\hat{\mathbf{v}}=N_{\mathrm{p}}$. From this we obtain the estimate

[TABLE]

Intuitively, this estimate simply averages the columns of $Y_{\mathrm{o}}$ and scales the result so that (10) holds. While this estimator is not optimal for Poisson noise, we have found that it works well given the large number of pixels.

As stated in (11) and (12), the measurement contains a uniformly dense region, $\Omega_{\mathrm{z}}$, and a open beam region, $\Omega_{\mathrm{0}}$. Then for each pixel in the region $\Omega_{\mathrm{z}}$, the theoretical transmission spectrum through the sample is given by

[TABLE]

where $\mathbf{q}(\mathbf{z})\in\mathbb{R}^{N_{\mathrm{A}}}$ is the vector of transmission resulting from the material areal densities, $\mathbf{z}$. For each pixel in the region $\Omega_{\mathrm{0}}$ we observe 100% transmission in the open beam region,

[TABLE]

where we use the fact that $\mathbb{1}^{\top}R=\mathbb{1}^{\top}$, shown in Appendix A.

Let $\boldsymbol{\upomega}_{\mathrm{0}}\in\mathbb{R}^{N_{\mathrm{p}}}$ and $\boldsymbol{\upomega}_{\mathrm{z}}\in\mathbb{R}^{N_{\mathrm{p}}}$ be a weighting vectors that take the averages over pixels in those regions, $\Omega_{\mathrm{0}}$ and $\Omega_{\mathrm{z}}$, respectively. More precisely, let

[TABLE]

Then define the following three averages:

[TABLE]

By using (17), (6), (7), and (8), we then have that

[TABLE]

Rearranging terms then yields the following result

[TABLE]

Next, by using (18), (5), (7), (8), and (15) along with the expression for $\boldsymbol{\upvarphi}$ from (21), we can derived the following expression for the average $\bar{\mathbf{y}}_{\mathrm{sz}}^{\top}$ given by

[TABLE]

Similar to (IV), an expression for $\bar{\mathbf{y}}_{\mathrm{s0}}^{\top}$ can be found using (19) and (16),

[TABLE]

By substituting in the measured neutron counts ($Y_{\mathrm{s}}$, $Y_{\mathrm{o}}$) for the expected measurements ($\bar{Y}_{\mathrm{o}}$, $\bar{Y}_{\mathrm{s}}$) in (IV) and (IV), we can formulate a loss function. Then we can estimate the nuisance parameters by minimizing this loss function:

[TABLE]

where the forward operators are defined by,

[TABLE]

and the scalar $\beta$ is a weighting parameter that balances the associated costs between solving (IV) and (IV). Since $\boldsymbol{\upomega}_{\mathrm{0}}^{\top}$ and $\boldsymbol{\upomega}_{\mathrm{z}}^{\top}$ perform averaging over many pixels, the Poisson distribution approximates a Gaussian via the Central Limit Theorem, which we use as a justification for the use of Least Squares estimators in (IV).

This optimization yields the estimates $\hat{\mathbf{z}}$, $\hat{\alpha}_{1}$, $\hat{\alpha}_{2}$, and $\hat{\boldsymbol{\uptheta}}$. The final parameter, $\boldsymbol{\upvarphi}$, is then estimated via (21) as

[TABLE]

and the estimates for the flux and background

[TABLE]

V Areal Density Reconstruction

In this section, we describe how we reconstruct the isotope areal density, $Z$, based on a statistical model of the neutron counting statistics. The estimated nuisance parameters from Section IV are considered known at this point.

We assume that the actual neutron sample measurements, $Y_{\mathrm{s}}$, are independent Poisson distributed with conditional mean $\bar{Y}_{\mathrm{s}}$. More precisely, if we define the operator

[TABLE]

then we assume that

[TABLE]

and the conditional mean of the observations given the unknown $Z$ is therefore

[TABLE]

For a Poisson random matrix $Y$ with mean $\varLambda$, the negative log-likelihood is

[TABLE]

Then our maximum likelihood estimate for the material areal density is given by

[TABLE]

VI Volumetric Density Reconstruction

In this section we describe a 3D tomographic extension to the 2D case that is possible when measurements of multiple rotational views are available. This will be a two step approach involving first reconstructing the areal densities, $Z$, for each view, and then performing tomographic reconstruction of the volumetric densities, $X$.

Let $\{Y_{\mathrm{s}}^{(k)}\}_{k=1}^{N_{\mathrm{r}}}$ be the set of $N_{\mathrm{r}}$ measurements at distinct views. Corresponding to these measurements, we will reconstruct $N_{\mathrm{r}}$ areal density views $\{Z^{(k)}\}_{k=1}^{N_{\mathrm{r}}}$ using the material decomposition of (30). We assume all nuisance parameters except $\alpha_{1}$ are independent of view index and thus only need to be computed once from a single view. The $\{\alpha_{1}^{(k)}\}_{k=1}^{N_{\mathrm{r}}}$ set of parameters are assumed to be dependent on the view index, $k$, and are estimated using the closed form expression777This expression is a corollary of (IV) and its derivation is omitted.

[TABLE]

where $[\mathbf{y}_{\mathrm{s0}}^{(k)}]^{\top}=\nicefrac{{([\boldsymbol{\upomega}_{\mathrm{0}}^{(k)}]^{\top}Y_{\mathrm{s}}^{(k)})}}{{([\boldsymbol{\upomega}_{\mathrm{0}}^{(k)}]^{\top}\hat{\mathbf{v}})}}$ correspond to the user selected open beam regions, $\Omega_{\mathrm{0}}^{(k)}$, while $\Omega_{\mathrm{z}}^{(k)}$ regions need not necessarily be available for every rotational view.

The estimated set of areal densities is stacked, and we assume the following decomposition

[TABLE]

where the areal density estimates, $\hat{Z}\in\mathbb{R}^{(N_{\mathrm{r}}N_{\mathrm{p}})\times N_{\mathrm{r}}}$, are in units of $\mathrm{mol}/\mathrm{cm}^{2}$, the tomographic projector, $A\in\mathbb{R}^{(N_{\mathrm{r}}N_{\mathrm{p}})\times N_{\mathrm{v}}}$, is in units of $\mathrm{cm}$, and the volumetric densities, $X\in\mathbb{R}^{N_{\mathrm{v}}\times N_{\mathrm{m}}}$ are in units of $\mathrm{mol}/\mathrm{cm}^{3}$ and comprised of $N_{\mathrm{v}}$ voxels. Then $X_{j,m}$ represents the volumetric density of the $j^{\mathrm{th}}$ voxel for the $m^{\mathrm{th}}$ isotope in the sample, and $A_{i,j}$ integrates along the tomographic projection line of the $j^{\mathrm{th}}$ voxel and the $i^{\mathrm{th}}$ measurement projection.

Finally, computing the volumetric densities, $X$, boils down to solving (32) which can be decomposed by columns:

[TABLE]

Consequently for each isotope the volumetric density $X_{*,m}$ is reconstructed from $\hat{Z}_{*,m}$ as the solution of

[TABLE]

using a tomographic reconstruction algorithm, where $r_{m}(\cdot)$ is a 3D regularization term of the $m^{\mathrm{th}}$ isotope. Although (33) specifies the use of MBIR reconstruction, other tomographic reconstruction methods can be used. However, we note that neutron tomography typically must use sparse views due to the slow view acquisition times, and that traditional CT reconstruction methods, such as filtered back projection, result in severe artifacts for the sparse view case. So we believe that MBIR reconstruction is a good choice for this application.

VII Optimization Methods

The key steps in the reconstruction of (IV) and (30) require the solution of non-convex optimization problems. In this section, we discuss the techniques we use to implement robust strategies for computing their solutions.

VII-A Preconditioning

A practical issue in optimization is that both (IV) and (30) incorporate loss functions with exponential terms such as $\mathbf{b}(\boldsymbol{\uptheta})^{\top}=\exp\{\boldsymbol{\uptheta}^{\top}P\}$ and $\mathbf{q}(\mathbf{z})^{\top}=\exp\left\{-\mathbf{z}^{\top}D\right\}R$. These terms can have large dynamic range and can therefore be sensitive to small changes in their arguments.

In order to limit the dynamic range of these exponential terms, we defined $P$ in (51) so that its rows have unit norm. We can do the same for the dictionary matrix $D$ by defining

[TABLE]

where $C$ is a diagonal matrix such that $C_{i,i}=\|D_{i,*}\|$. However, this also requires that we define the new functions

[TABLE]

that are used to replace the functions of (15) and (29).

After the modified optimization problems are solved, the solution of the original problem can be found using

[TABLE]

VII-B Initialization

Since the optimization problems are not convex, it is important to select good initializations that lead to solutions that are close to the global minimum. We first focus on the initialization of the optimization in (IV) of the nuisance parameter estimation. The initial values that we use are

[TABLE]

where $(\cdot)^{\dagger}$ denotes the Moore–Penrose pseudoinverse and

[TABLE]

is the average measured transmission, given the initial nuisance parameters, $\alpha_{1}$, $\alpha_{2}$, and $\boldsymbol{\uptheta}$. This initialization is motivated by (IV) and the following crude approximations:

[TABLE]

For the optimization of (30), we initialize

[TABLE]

where the measured transmission

[TABLE]

is computed given the final nuisance parameters.

VIII Experimental Results

In this section we present results of the TRINIDI algorithm applied to both simulated and experimentally measured data. Section VIII-A outlines the experimental settings and notes about the evaluation of the results. In Section VIII-B, we demonstrate the method using 2D simulated radiographs. In Section VIII-C, we provide similar 2D results for experimentally measured data using a known phantom. Finally in Section VIII-D, we present results for fully tomographic 3D reconstruction of nuclear fuel pellets that are verified using independent mass spectroscopy measurements.

VIII-A Methods

Experiments in this work were performed at the Flight-Path-5 beamline at LANSCE [10]. In this facility, high energy spallation neutrons are slowed down by a high intensity ambient temperature water moderator resulting in a $20\,\mathrm{Hz}$ pulsed neutron beam with a wide energy spectrum. A detailed description of this facility including a characterization of the incident neutron spectrum can be found in [22].

The neutrons travel along a flight path of length of $L=10.4\,\mathrm{m}$ through the sample onto a TOF imaging detector. The neutron sensitive multi-channel plate detector works in conjunction with four Timepix readout chips, described in more detail in [45]. The detector has $N_{\mathrm{p}}=512\!\times\!512$ pixels with a $55\!\times\!55\,(\mathrm{\upmu m})^{2}$ pixel pitch i.e. covering a field of view of approximately $28\!\times\!28\,(\mathrm{mm})^{2}$, and frames are recorded at a rate of, $\nicefrac{{1}}{{\Delta_{\mathrm{t}}}}\approx\nicefrac{{1}}{{30\mathrm{ns}}}\approx 33.3\,\mathrm{MHz}$.

The resolution function is modeled as an energy dependent weighted sum of two chi-squared probability density functions based on the work in [46]. As described in Appendix A, the energy dependent pulse shape is approximated by a weighted sum of convolutions of $K$ kernels. We choose $K=5$ since we found that this gave sufficiently accurate representation of the resolution function.

In order to assess the accuracy of the nuisance parameter estimation of (IV), we plot $\mathbf{y}_{\mathrm{s0}}$ (the measurement in the open beam region) and the fit of $\mathbf{y}_{\mathrm{s0}}$ (the effective open beam) given by

[TABLE]

where (44) is the expected measurement when there is $100\,\%$ transmission (i.e. zero material). We also plot $\mathbf{y}_{\mathrm{sz}}$ (the measurement in the uniformly dense region) and the fit of $\mathbf{y}_{\mathrm{sz}}$ given by

[TABLE]

where (45) is the expected measurement when there is a uniform material density, $\mathbf{z}$, attenuating the beam. Lastly, we plot the effective background given by

[TABLE]

where (46) is the expected measurement when there is there is $0\,\%$ transmission (i.e. infinite material).

We display normalized density images for more accessible interpretation of our density estimates. The scales of the densities differ widely from one another so that instead of directly displaying the all $Z_{*,m}$ with a different color map, we choose to display all $Z_{*,m}/\mathbf{z}_{m}$ with the same color map. The scalars $\mathbf{z}_{m}$ are the ground truth average density of the $m^{\mathrm{th}}$ isotope (or some other adequate normalization). Thus, the normalized densities, $Z\mathop{\mathrm{diag}}(\mathbf{z})^{-1}$, are unitless, where $1$ corresponds to the values of $\mathbf{z}$ i.e. ground truth. A similar treatment is applied to the volumetric densities, $X$.

Finally, for the cross section dictionary, $D$, we use tabulated data, of the neutron total cross section at a temperature of $293.6\ \mathrm{K}$ from ENDF/B-VIII.0[44].

We compare our reconstruction results to the reconstruction using the linear model from [37] which we refer to as linear baseline reconstruction or $Z^{\mathrm{lin}}$. In an attempt to further compare our results to existing methods, we tried processing the transmission spectra from our experiments using SAMMY [29]. However, due to the significant noise, we were unable to get the software to function properly.

All tomographic reconstructions were performed using the super-voxel model-based iterative reconstruction (SVMBIR) software package [47], with parameters selected to maximize subjective image quality, and all optimizations were performed using the accelerated proximal gradient method (APGM) with robust line search888This algorithm was selected due to the convenient availability of a reliable implementation rather than any expectation of optimum suitability for this problem, and it is expected that other optimization algorithms may provide faster convergence. In particular, after completion of the computational experiments reported here, we observed that the Broyden–Fletcher–Goldfarb–Shannon (BFGS) algorithm appears to converge significantly faster than APGM. as part of the scientific computational imaging code (SCICO) software package [48]. The $\beta$ parameter from (IV) is chosen to be equal to $1$ in Section VIII-B and Section VIII-D and chosen to be equal to [math] in Section VIII-C because in this experiment, there is no $\Omega_{\mathrm{0}}$ region available. The number of background basis functions is chosen to be $N_{\mathrm{b}}=5$.

VIII-B 2D Reconstruction from Simulated Data

In this experiment, we reconstruct a 2D phantom from Monte-Carlo simulated data. Figure 4 illustrates the structure of the simulated 2D phantom which is composed of a series of overlapping disks each formed by a distinct material isotope. The other parameters of the simulation were chosen to approximately mimic a typical experimental scenario.

The ground truth areal densities were chosen to represent $N_{\mathrm{m}}=5$ distinct isotopes, where each areal density map, $Z_{*,m}$, consists of a disk of constant density $\mathbf{z}_{m}$ with a surrounding area of zero material. Each map has $N_{\mathrm{p}}=128\!\times\!128$ pixels. The numerical values of the densities, $\mathbf{z}$, can be found in Table I and were chosen to yield a visible distribution of resonances in the measurement spectra. The isotopes used are ${}^{238}\mathrm{U}$, ${}^{239}\mathrm{Pu}$, ${}^{240}\mathrm{Pu}$, ${}^{181}\mathrm{Ta}$, ${}^{241}\mathrm{Am}$ and their respective cross sections are shown in Figure 6(a).

The simulated measurement counts, $Y_{\mathrm{s}}$ and $Y_{\mathrm{o}}$, were generated by first computing their sample means using (5) and (6), and then generating pseudo-random Poisson samples corresponding to those means. In order to make the simulation representative of a typical physical experiment, the values for $\mathbf{v}$, $\boldsymbol{\uptheta}$, $\alpha_{1}$, $\alpha_{2}$, $\boldsymbol{\upvarphi}$, and the TOF sampling interval are chosen to be the same as the estimates from the real data covered in Section VIII-C.

Figure 5 shows the summed areal densities over all isotopes, the TOF integrated measurement ratio $\nicefrac{{Y_{\mathrm{s}}\mathbb{1}}}{{Y_{\mathrm{o}}\mathbb{1}}}$, and the regions $\Omega_{\mathrm{z}}$ and $\Omega_{\mathrm{0}}$ used for this experiment. The set $\Omega_{\mathrm{z}}$ was taken to be the region in the center of the phantom in which all the disks overlap, and the set $\Omega_{\mathrm{0}}$ was taken to be the region surrounding the disks that is not covered by any material.

Table I lists the results of the nuisance parameter estimation using (IV). Notice that the parameters of interest, $\hat{\mathbf{z}}$, are accurately estimated within approximately $3\%$. The estimation of the nuisance parameters, $\alpha_{1}$, $\alpha_{2}$ and $\boldsymbol{\uptheta}$, is more variable; however, we note that accurate estimation of nuisance parameters may not always be necessary as long as the effective open beam and background are accurately modeled.

Figure 6(b) shows plots corresponding to the fitted curves of the nuisance parameter estimation. The orange and blue curves are the average sample measurements, $\mathbf{y}_{\mathrm{s0}}$ and $\mathbf{y}_{\mathrm{sz}}$, and the dashed lines of the same respective colors are their corresponding fits. Note that the measurements and fits visually align very well. The red dashed line is the resulting effective background.

Figure 7 shows a direct comparison to the ground truth spectra. The effective open beam (yellow lines) and background (red lines) are compared directly their respective ground truth spectra (dashed lines). Notice that all three fits are quite accurate, even though the underlying nuisance parameter estimates in Table I where not as accurate.

Figure 8 shows the objective of the optimization as a function of iteration number. Although the optimization takes many iterations and is not monotone, the objective eventually converges to approximately the same value as the ground truth.

Figure 9(a) shows the normalized ground truth areal densities, Figure 9(b) shows the linear baseline [37] reconstructed areal density, and Figure 9(c) shows the corresponding TRINIDI reconstructed areal density. The linear baseline reconstruction deviates significantly from the ground truth. Notice that using TRINIDI each isotope map is accurately reconstructed from the simulated data. Also, even though the disks significantly overlap with one another, especially in the center of the field of view, there are no visible artifacts, which is not the case for the linear baseline.

It is worth pointing out that the cross section dictionary, $D$, contains all five isotopes, however, most pixels in the ground truth do not have all five isotopes present. For example, in the very top-center region only ${}^{238}\mathrm{U}$ is present. The TRINIDI reconstructions show that ${}^{238}\mathrm{U}$ is reconstructed accurately as non-zeros while the other isotopes are reconstructed accurately as (close to) zeros. This indicates that overpopulation of the cross section dictionary will likely not result in degradation of the reconstruction, however, we suggest that unnecessary overpopulation should be avoided.

Table II shows the mean and standard deviation of the density estimates, $\hat{Z}_{*,m}$ for each of the disks of material individually compared to the ground truth ($Z$) and the linear baseline ($Z^{\mathrm{lin}}$). Note that even though there is significant noise, the mean estimates using TRINIDI are very close to the ground truth densities. Also note that the different isotopes have different estimation accuracy. In comparison, the linear baseline estimates are not quantitatively accurate.

Figure 10 shows the convergence behaviour of the optimization of (30). The final objective is slightly better than the ground truth objective after only about 50 iterations and then converges to its final value. This indicates likely convergence to near the global minimum, however, with possible slight overfitting. Note that the convergence takes significantly fewer iterations as compared to the nuisance parameter estimation.

VIII-C 2D Reconstruction from Experimental Data

This experiment is similar to the 2D Monte-Carlo simulated experiment of Section VIII-B, but it uses measured data. We use a known, well-defined sample in order to demonstrate the quantitative accuracy of our method with experimentally measured data.

The sample is a stack of uniform tantalum ($\mathrm{Ta}$) and tungsten ($\mathrm{W}$) plates. The plates have thicknesses of $\ell_{\mathrm{Ta}}=2.42\,\mathrm{mm}$ and $\ell_{\mathrm{W}}=1.75\,\mathrm{mm}$. We refer to this sample as the $\mathrm{Ta}\mathrm{W}$ sample. The $N_{\mathrm{A}}=2260$ TOF bins span an interval from $70.11\,\mathrm{\upmu s}$ up to $739.1\,\mathrm{\upmu s}$ which corresponds to an energy range of $115.0\,\mathrm{eV}$ down to $1.04\,\mathrm{eV}$, respectively. We use a cropped region of the detector of $N_{\mathrm{p}}=128\!\times\!128$ pixels.

Figure 11 shows the TOF integrated measurement ratio, $\nicefrac{{Y_{\mathrm{s}}\mathbb{1}}}{{Y_{\mathrm{o}}\mathbb{1}}}$ of the sample. The sample covers the entire field of view and thus there is no set of pixels that is suitable for the open beam region $\Omega_{\mathrm{0}}$. For this reason we choose the parameter $\beta=0$ in (IV), effectively dropping the constraint of jointly solving (IV) in the nuisance parameter estimation. Since the sample is assumed to be constant across the field of view, the set $\Omega_{\mathrm{z}}$ is chosen to be the entire field of view.

Figure 12(a) shows the cross section dictionary, $D$, used for this experiment. Since tantalum in natural abundance consists of practically $100\,\%$ of the ${}^{181}\mathrm{Ta}$ isotope, we use this isotope for the cross section dictionary. In contrast, natural tungsten consists of several isotopes, so for its cross section we use a proportionally weighted sum of the cross section of its isotopes, ${}^{182}\mathrm{W}$, ${}^{183}\mathrm{W}$, ${}^{184}\mathrm{W}$, and ${}^{186}\mathrm{W}$.

Figure 12(b) shows the plots corresponding to the nuisance parameter estimation of (IV). Note that (compared to Figure 6), the fit of $\mathbf{y}_{\mathrm{s0}}$ is not shown since there is no $\Omega_{\mathrm{0}}$ region. Nevertheless, the optimization seems to result in a good fit and plausible spectra of the effective flux and background.

Table III shows the numerical values of the nuisance parameters. Note that the $\alpha_{2}$ parameter is relatively different from $1$, indicating that the effects of the sample on the background are important to model.

Figure 13 shows the areal density estimates, $\hat{Z}\mathop{\mathrm{diag}}(\mathbf{z})^{-1}$. We observe that, apart from some noise, the estimates are visually relatively constant as is expected from a spatially constant sample. The noise level appears to be similar to the simulated data of Section VIII-B, however, it seems to be more spatially correlated, indicating possible detector cross talk.

Since the $\mathrm{Ta}\mathrm{W}$ sample matches the idealized example from Section II we can use (2) and the areal density estimates, to make a thickness estimate for each of the plates. Table IV shows the mean and standard deviation computed from the areal density estimates, $\hat{Z}$, the resulting estimated thicknesses and the ground truth thicknesses of the plates. The estimated thicknesses are relatively close to the actual thicknesses with an error of less than $13\%$.

We believe that the slight inaccuracies in Table III may be primarily inaccuracies in the assumed values of $R$ and $D$. We rely on previous measurements by [46] for the parameters of $R$, which may not be sufficiently accurate. In addition, for the cross sections, $D$, we use tabulated data from ENDF/B-VIII.0[44], which is empirically compiled and may not be sufficiently accurate for our experiment. In order to support our conjecture, Figure 14 highlights a region (red box) in the spectrum where there are likely errors in $D$. In this region, the fitted resonance energy (blue dashed) in ${}^{181}\mathrm{Ta}$ near $200\,\mathrm{\upmu s}$ has a slightly lower TOF than the measured resonance (blue solid). This type of mismatch can significantly affect the density estimates. There are likely other errors in the cross sections, which may however not be as obvious to spot as this one.

VIII-D 3D Reconstruction from Experimental Data

In this experiment we demonstrate the tomographic reconstruction extension that is described in Section VI on measured experimental data of a transmutation nuclear fuel pellet sample. We will also compute average volumetric mass densities and validate them through independent mass spectroscopy measurements.

Figure 16 shows the TOF integrated measurement ratio, $\nicefrac{{Y_{\mathrm{s}}^{(k)}\mathbb{1}}}{{Y_{\mathrm{o}}\mathbb{1}}}$, of a single rotational view. The cylindrical nuclear fuel pellet is surrounded by a double wall steel cladding. We choose the regions surrounding the sample as the $\Omega_{\mathrm{0}}$ region and a thin axial region the center of the sample as the $\Omega_{\mathrm{z}}$ region in the radiograph selected for the nuisance parameter estimation.

Figure 17 shows the cross section dictionary, $D$. We assume that the sample contains the resonant isotopes ${}^{237}\mathrm{Np}$, ${}^{238}\mathrm{U}$, ${}^{239}\mathrm{Pu}$, ${}^{240}\mathrm{Pu}$, ${}^{241}\mathrm{Am}$. We also assume there are non-resonant materials such as for example $\mathrm{Zr}$ in the fuel and steel in the double wall cladding. For this reason, the cross section dictionary, $D$, not only includes the resonant isotopic cross sections but also the mostly constant ${}^{1}\mathrm{H}$ (hydrogen) cross section. The reconstructed areal density associated with ${}^{1}\mathrm{H}$ is interpreted as hydrogen-equivalent areal density, regardless whether the sample contains any hydrogen.

The sample is imaged at $N_{\mathrm{r}}=101$ different rotational views. The $N_{\mathrm{A}}=2440$ TOF bins span an interval from $47.6\,\mathrm{\upmu s}$ up to $791.7\,\mathrm{\upmu s}$ which corresponds to an energy range of $249.8\,\mathrm{eV}$ down to $0.902\,\mathrm{eV}$, respectively. We use a cropped region of the detector of $N_{\mathrm{p}}=256\!\times\!512$ pixels. We applied a $4\!\times$ binning because the measurements suffer from severe defects such as many malfunctioning detector rows and columns. We understand that in general binning may not be the best way to remove such defects, however, for simplicity sake we deemed this procedure appropriate.

Figure 17(b) shows the fits resulting from the nuisance parameter estimation of one single view. The good fits of the respective measurements indicate successful estimation of the parameters. Table V lists the estimated parameters.

Figure 18 shows the normalized areal density reconstruction, $\hat{Z}^{(k)}\mathop{\mathrm{diag}}(\hat{\mathbf{z}})^{-1}$, of a single view from the material decomposition. Each of the resonant isotopes, ${}^{237}\mathrm{Np}$, ${}^{238}\mathrm{U}$, ${}^{239}\mathrm{Pu}$, ${}^{240}\mathrm{Pu}$, ${}^{241}\mathrm{Am}$, clearly resolves the cylindrical fuel pellet while the ${}^{1}\mathrm{H}$ equivalent density is mostly concentrated in the cladding.

Figure 19 shows a single axial slice of the normalized volumetric density reconstruction, $\hat{X}\mathop{\mathrm{diag}}(\hat{\mathbf{x}})^{-1}$, from (33). We computed the vector $\hat{\mathbf{x}}\in\mathbb{R}^{N_{\mathrm{m}}}$ as the average estimated volumetric density in the mask shown in Figure 16. The the values of $\hat{\mathbf{x}}$ are shown in Table VI. In the displayed slice the clear separation between fuel and cladding is visible. Also, there is a void visible in the 7-o’clock position in all resonant isotope maps, indicating what is probably a void in the sample.

Finally, Table VII shows the average reconstructed mass density estimates, $\hat{\boldsymbol{\uprho}}\,^{[\mathrm{g/cm^{3}}]}$, computed from the volumetric densities, $\hat{\mathbf{x}}\,^{[\mathrm{mol}/\mathrm{cm}^{3}]}$ and the relationship in (3). The table also lists the values of ground truth measurements obtained using mass spectroscopy. Note that the reconstructed estimates are very close to the ground truth values which provides strong experimental validation of our method.

IX Conclusions

We present a novel method for both 2D and 3D reconstruction of material areal and volumetric density from TOF neutron images. Our method, TRINIDI, accounts for a number of important measurement effects including spatially varying background, non-uniform flux, finite pulse duration, and Poisson counting noise. Especially, when including the finite pulse duration in the forward model, the problem cannot be linearized in the same fashion that is common with Beer-Law attenuation models. TRINIDI is based on a 2-step process in which we first estimate nuisance parameters associated with background and flux, and we then reconstruct the unknown isotopic areal density. For the 3D case, we then perform tomographic reconstruction of the views to form a 3D volumetric density estimate. Both our simulated and experimental results indicate that TRINIDI can reconstruct quantitatively accurate estimates in both 2D and 3D. Notably, comparisons of our 3D material decomposed reconstructions for nuclear fuel pellets are accurate when compared to mass spectroscopy measurements.

Appendix A Resolution Function

The resolution operator, $R$, models the uncertainty of emission time of the neutrons given their known neutron energy. Equivalently, it models the conditional probability of TOA’s given the neutron’s TOF’s, stated more precisely in (4). Since the pulse shape changes with every neutron energy, we choose to approximate the resolution operator using a weighted sum of $K$ convolution operations, so that

[TABLE]

The $R^{k}$ matrix corresponds to the convolution operation with kernel $\mathbf{r}^{k}$, which in turn corresponds to the pulse shape at the TOF bin with index $\lfloor k\frac{N_{\mathrm{F}}-1}{K-1}\rfloor$. The $W^{k}$ matrices are diagonal weighting matrices with diagonals $\mathbf{w}^{k}$ which blend together the blurring operators, $R^{k}$.

Figure 20 illustrates the application of the resolution operator with $K=5$ components in a typical scenario. In Figure 20(a) the dashed lines are the diagonals, $\mathbf{w}^{k}$ and the solid lines are the corresponding kernels, $\mathbf{r}^{k}$ centered at the TOF where they are most active. The blue line in Figure 20(b) shows a typical non-blurred transmission spectrum, $\exp(-\mathbf{z}^{\top}D)$. The orange line shows the corresponding blurred transmission spectrum, $\exp(-\mathbf{z}^{\top}D)R$. Note that the blurring is weaker in the low-TOF region and stronger in the high-TOF region and that black resonances in the non-blurred signal are mostly not opaque after the blurring.

Next we show that $\mathbb{1}^{\top}R=\mathbb{1}^{\top}$. First, we select the weights so every TOA bin they sum to $1$, i.e.

[TABLE]

Second, we assume that $\sum_{j}\mathbf{r}^{k}_{j}=1$ so that every emitted neutron is measured as some TOA bin. As explained in Section III this is achieved by choosing the number of rows of $R^{k}$ to be sufficiently large so that

[TABLE]

so we have that

[TABLE]

Finally, using (47), (48), and (49) we see that

[TABLE]

Appendix B Spectral Basis Matrix of the Background

In this work, we choose the background basis to be

[TABLE]

where $a_{n}=\sqrt{\sum_{k=0}^{N_{\mathrm{A}}-1}[\left(\log\left[k\Delta+k_{0}\right]\right)^{n}]^{2}}$ scales the rows of $P$ so they have unit norm and $0\leq n<N_{\mathrm{b}}$, $0\leq k<N_{\mathrm{A}}$, and we have found that scalar $\Delta=\frac{\mathrm{e}^{1}-\mathrm{e}^{-1}}{N_{\mathrm{A}}-1}$ and the offset $k_{0}=\mathrm{e}^{-1}$ work well across a range of cases.