Autonomous Polycrystalline Material Decomposition for Hyperspectral Neutron Tomography

TL;DR

This paper presents an autonomous algorithm for decomposing and localizing polycrystalline materials in hyperspectral neutron tomography, enabling non-destructive analysis of complex crystalline structures with high accuracy.

Contribution

The paper introduces an autonomous material decomposition algorithm that estimates attenuation spectra and reconstructs 3D material volumes from hyperspectral neutron data, advancing non-destructive crystalline analysis.

Findings

Accurately estimates attenuation spectra from simulated and experimental data.

Successfully decomposes polycrystalline structures in 3D reconstructions.

Demonstrates effectiveness on both simulated and real neutron tomography data.

Abstract

Hyperspectral neutron tomography is an effective method for analyzing crystalline material samples with complex compositions in a non-destructive manner. Since the counts in the hyperspectral neutron radiographs directly depend on the neutron cross-sections, materials may exhibit contrasting neutron responses across wavelengths. Therefore, it is possible to extract the unique signatures associated with each material and use them to separate the crystalline phases simultaneously. We introduce an autonomous material decomposition (AMD) algorithm to automatically characterize and localize polycrystalline structures using Bragg edges with contrasting neutron responses from hyperspectral data. The algorithm estimates the linear attenuation coefficient spectra from the measured radiographs and then uses these spectra to perform polycrystalline material decomposition and reconstructs 3D material volumes to localize materials in the spatial domain. Our results demonstrate that the method can accurately estimate both the linear attenuation coefficient spectra and associated reconstructions on both simulated and experimental neutron data.

Equations

$$y_{r,c,k}^o=\sum _{i=-\infty }^{\infty }\sum _{j=-\infty }^{\infty }{\bar{y}}_{i,j,k}^o{h}_{r-i,c-j},$$

$$p_{v,r,c,k}=-\log\bigg{(}\frac{y_{v,r,c,k}}{y_{r,c,k}^{o}}\bigg{)}\ .$$

$$p_{v,r,c,k}=p_{v,r,c,k}^{\text{meas}}-b_{v,k}.$$

$$(V^s,D^s)=\underset{(V\ge 0,D\ge 0)}{\mathrm{argmin}}\{\Vert p-V{D}^t{\Vert }_F^2\}.$$

$${\hat{x}}_j^s=\underset{x_j^s}{\mathrm{argmin}}\{f(x_j^s)+h(x_j^s)\},$$

$$f(x_j^s)=\frac{1}{2{\sigma }_v^2}\Vert V_j^s-A{x}_j^s{\Vert }^2,$$

$$D^m=D^s{T}^t.$$

$$V^m=\underset{V\ge 0}{\mathrm{argmin}}\{\Vert p-V(D^m{)}^t{\Vert }_F^2\}.$$

$${\hat{x}}_j^m=\underset{x_j^m}{\mathrm{argmin}}\{f(x_j^m)+h(x_j^m)\},$$

Full text

We introduce an autonomous material decomposition (AMD) algorithm to automatically characterize and localize polycrystalline structures using Bragg edges with contrasting neutron responses from hyperspectral data. The algorithm estimates the linear attenuation coefficient spectra from the measured radiographs and then uses these spectra to perform polycrystalline material decomposition and reconstructs 3D material volumes to localize materials in the spatial domain. Our results demonstrate that the method can accurately estimate both the linear attenuation coefficient spectra and associated reconstructions on both simulated and experimental neutron data.

**Index Terms— ** neutron computed tomography, material decomposition, non-negative matrix factorization, clustering

1 Introduction

Hyperspectral neutron computed tomography (HSnT) is a powerful technique that has enabled 3D imaging of important material properties such as crystallographic phases [1] and isotopic compositions [2]. HSnT is typically carried out at pulsed neutron source facilities by rotating a sample in the path of the beam, measuring the transmitted projection using a high-efficiency time-of-flight detector corresponding to each orientation, and processing the resulting data using a reconstruction algorithm. The combination of the pulsed source and time-of-flight detector enables resolving each projection into thousands of wavelength bins depending on the travel time of the measured neutrons [3].

A typical workflow for HSnT involves processing the data corresponding to each wavelength bin using a reconstruction algorithm like filtered back projection (FBP), extracting the resulting spectra for each voxel in 3D [1], and identifying homogeneous regions based on their spectral signature. For polycrystalline materials, these spectra correspond to the linear attenuation coefficient (or $\mu$-spectra) of a certain phase and are characterized by sharp dips in the attenuation values (called Bragg edges) at certain wavelengths. However, due to the low signal-to-noise ratio in each bin that is typical for HSnT data sets, the resulting reconstructions can be extremely noisy. Furthermore, it is inefficient to process each wavelength independently when the underlying sample is only composed of a few phases, each with its unique spectral signature.

In order to address some of the challenges with existing approaches, dimensionality reduction methods have been applied in the context of HSnT. In [2, 4], the authors have used a few known spectral signatures as a basis to decompose the measured HSnT data (in the epithermal range) and performed the reconstruction of this lower dimensional dataset, thereby significantly reducing the computation when compared to processing each of the wavelength bins independently. A recent work adopted a method of computing the $\mu$-spectra for polycrystalline materials utilizing known material locations in the hyperspectral reconstructions and using them for material decomposition [5]. However, none of these approaches are fully automated, as they use some form of supervision (known spectra or material locations).

In this paper, we present a novel algorithm for autonomous material decomposition (AMD) on hyperspectral neutron data without prior information. A key component of our algorithm is a two-stage subspace extraction procedure that transforms high-dimensional hyperspectral measurements into low-dimensional subspace projections. Stage one maps into an intermediate dimensional subspace that is designed to promote accurate, meaningful segmentation while rejecting most of the inherent noise and simultaneously drastically reducing the computation needed for reconstruction. Stage two uses this segmentation to produce the actual material decomposition. A second key component of AMD is an unsupervised clustering procedure that automatically segments the material regions used for $\mu$-spectra estimation. We test on both simulated and measured data and provide quantitative and qualitative performance evaluations.

2 Methodology

Our AMD algorithm calculates $\mu$-spectra and material reconstructions from hyperspectral neutron radiographs. It uses low-noise open-beam images (also known as blank scans or flat fields) to normalize object projections. The algorithm decomposes these normalized hyperspectral projections using subspace basis vectors and reconstructs within this subspace, then performs clustering operations to spatially segment material regions. The individual $\mu$-spectra are calculated based on mean voxel values for each material region, after which the projection domain material decomposition and material reconstruction are performed. The entire procedure is illustrated in Figure 1.

2.1 Open-Beam Processing

To achieve the highest signal-to-noise ratio (SNR) possible for the open-beam images, we first average over multiple open-beam sets to reduce the noise. Then we filter the open-beam images at each wavelength using a 2D normalized Hamming kernel $h$ for further smoothing. This operation can be expressed as 2D convolution given by

[TABLE]

where $y^{o}\in\mathbb{R}^{N_{r}\times N_{c}\times N_{k}}$ is the smooth open-beam array, $\bar{y}^{o}\in\mathbb{R}^{N_{r}\times N_{c}\times N_{k}}$ is the averaged open-beam array, $N_{r}$ is the number of rows in the detector, $N_{c}$ is the number of columns in the detector, and $N_{k}$ is the number of energies/wavelengths.

2.2 Projection Density Computation

To convert neutron counts to projection densities, we need to compute the negative log of the ratios of the object projections and the corresponding open-beam images. Ideally, the open-beam images $y^{o}$ would have the same neutron dosage rate as the object projections $y\in\mathbb{R}^{N_{v}\times N_{r}\times N_{c}\times N_{k}}$, where $N_{v}$ is the number of views. Then the ideal projection density array $p\in\mathbb{R}^{N_{v}\times N_{r}\times N_{c}\times N_{k}}$ can be expressed as

[TABLE]

However, in actual experiments, the dosage rate might vary for each view and energy, making the measured projection density $p^{\text{meas}}\in\mathbb{R}^{N_{v}\times N_{r}\times N_{c}\times N_{k}}$ deviate from the ideal. In order to correct for this, we subtract a view and wavelength-dependent offset $b\in\mathbb{R}^{N_{v}\times N_{k}}$ that makes the background approximately zero. This gives

[TABLE]

We estimate each $b_{v,k}$ by taking the spatial average of $p^{\text{meas}}_{v,k}$ in a region outside the sample, where $p_{v,k}$ is assumed to be 0.

2.3 Subspace Extraction

To reduce the effect of noise and promote accurate material segmentation, we first use non-negative matrix factorization (NMF) [6] to reduce from full hyperspectral measurements to a lower dimensional subspace representation. More precisely, we seek to approximate $p=VD^{t}$, where $p\in\mathbb{R}^{N_{p}\times N_{k}}$ is the matrix obtained by using one hyperspectral projection density per row, $D\in\mathbb{R}^{N_{k}\times N_{s}}$ is a matrix whose columns are non-negative basis vectors for a subspace, and $V\in\mathbb{R}^{N_{p}\times N_{s}}$ gives the representation of each projection in subspace coordinates. Here $N_{p}=N_{v}\times N_{r}\times N_{c}$ is the number of projections, and $N_{s}$ is the dimension of the subspace. Taking $N_{m}$ to be the number of materials, we choose $N_{s}$ so that $N_{m}<N_{s}<<N_{k}$. The resulting decomposition can be obtained by solving the optimization problem given by

[TABLE]

Since $N_{s}<<N_{k}$, this reconstruction is significantly faster than full hyperspectral reconstruction.

2.4 Subspace Reconstruction

Space domain analysis requires us to perform volumetric reconstructions. The subspace representation of $p$ allows us to perform high-quality reconstructions within a reasonable time using SVMBIR, a sophisticated Python package for super-voxel-based iterative CT reconstruction [7, 8]. The reconstructed subspace volume $x^{s}\in\mathbb{R}^{N_{r}\times N_{c}\times N_{c}\times N_{s}}$ is obtained by restricting to each subspace component, $x_{j}^{s},j=1,\ldots,N_{s}$, one at a time and solving the associated optimization problem

[TABLE]

where $f(x_{j}^{s})$ is the forward model, and $h(x_{j}^{s})$ is the prior model. SVMBIR is used to perform this optimization using Q-Generalized Gaussian Markov random field (QGGMRF) as the prior model [9] and a forward model of the form

[TABLE]

where $A$ is the linear projection operator from a 3D volume to a rasterized sinogram, $V_{j}$ is the $j$th column of $V$, and $\sigma_{v}$ is the assumed standard deviation of the noise.

2.5 Material Decomposition and Reconstruction

To convert from the subspace to a physically meaningful material basis, we first use GMCluster [10], a Gaussian mixture model-based clustering package, to segment the voxels in $x^{s}$ to identify regions corresponding to each material. Within each segmented region, we compute the vector mean to obtain a row of an $N_{m}\times N_{s}$ matrix $T$. With $T$, we have the relationship between subspace and material basis dictionaries given by

[TABLE]

We scale $D^{m}$ by a factor of $\frac{1}{\delta}$ to obtain the set of $\mu$-spectra, where $\delta$ is the voxel thickness/pixel pitch.

To perform material decomposition, we need to rerun NMF on the projection densities. However, this time, we calculate only $V^{m}$ using fixed $D^{m}$. The optimization problem is given by

[TABLE]

Finally, we estimate the desired material decomposition reconstruction $\hat{x}^{m}$ using

[TABLE]

where now $f$ is defined as in (6) with $V^{m}$ in place of $V^{s}$.

3 Results

We implemented our AMD algorithm on both simulated and measured data. Each dataset contains neutron radiographs with 1200 wavelength bins ranging between 1.5 to 4.5 Å. We used $64\times 64$ radiographs for simulated data and $512\times 512$ for real data. The numbers of projection angles are 32 and 27 for simulated and measured data, respectively.

The scanned physical phantom in Fig. 6 is formed from powdered nickel (Ni) and powdered copper (Cu) in an aluminum (Al) structure. Data simulation was performed using a digital phantom that roughly approximated the physical phantom. The neutron measurements were then simulated using a Poisson noise model and attenuation calculations based on the object geometry and the theoretically known $\mu$-spectra computed using the Bragg-edge modeling (BEM) library [11]. Heuristically, we found that 3 times the number of materials worked well as the subspace dimension for the cases we considered. So for all results, we used $N_{m}=3$ to represent the three materials (Ni, Cu, and Al) and $N_{s}=N_{m}\times 3=9$ for subspace extraction. For the purpose of clustering, we considered the empty region surrounding the phantom as an additional material but discarded it later.

We compared our AMD results to a baseline in which each of the $N_{k}=1200$ wavelengths is reconstructed using the FBP method, after which the $\mu$-spectra are computed by averaging over a manually selected region of that material; we call this method reconstruction domain material decomposition (RDMD). We note that RDMD is much more computationally expensive since it requires the computation of $N_{k}+N_{m}=1203$ reconstructions rather than the $N_{s}+N_{m}=12$ reconstructions required for our proposed AMD method. Due to the high computational cost, we reconstructed only a single slice for each of the hyperspectral projections ($N_{k}$), while all slices for the material projections ($N_{m}$) in the implementation of the RDMD method.

Figure 2 shows the results of $\mu$-spectra estimation for Ni, Cu, and Al in both the simulated and measured cases. Fig. 2(a) shows the ground truth used in simulation; Fig. 2(b) shows the spectra estimated from simulated data using AMD; Fig. 2(c) shows the spectra estimated from measured data using AMD; and Fig. 2(d) shows the spectra estimated from simulated data using the alternative RDMD method.

The AMD estimates of $\mu$-spectra accurately match the theoretical ground truth in the simulated case with an NRMSE of under 1% for Ni and Cu; and match up to a scaling factor in the measured case, which we believe may be due to compaction of the powdered material; this scaling factor also appears in the RDMD results. The AMD method extracts much cleaner and more accurate curves than the RDMD method. In particular, the Bragg edge locations in the AMD estimates of $\mu$-spectra closely match the ground truth.

Fig. 3(a) shows some ground truth slices, and Fig. 3(b) shows the associated AMD reconstructions from simulated data; Fig. 4 shows 3D renderings of these reconstructions. Notice that the AMD material decomposition of the phantom for simulated data is close to accurate. Fig. 5(a) and (b) show the reconstructions from measured data using AMD and the alternative RDMD. The AMD results are substantially more accurate than the RDMD results, likely due to the more accurate estimates of the $\mu$-spectra. Again, Fig. 6 shows a 3D rendering of the AMD reconstructions. For both simulated and measured data, the artifacts in the Al reconstructions can be attributed to the weak presence of Al in neutron images making it difficult to separate Al from the background.

4 Conclusion

We introduce an autonomous material decomposition (AMD) method based on a 2-stage subspace extraction approach in the projection domain. The method identifies regions of homogeneous polycrystalline material compositions and extracts accurate linear attenuation coefficients spectra for each material resulting in accurate volumetric reconstructions. By dramatically reducing the dimensionality of the hyperspectral data, AMD both reduces computation and improves the quality of reconstruction relative to baseline. While AMD is primarily devoted to Bragg edge imaging, we note that with some modifications, it may be useful for other hyperspectral modalities, such as neutron resonance imaging.

5 Acknowledgment

G. Buzzard was partially supported by NSF CCF-1763896, and C. Bouman was partially supported by the Showalter Trust. This research used resources at the Spallation Neutron Source, a DOE Office of Science User Facility operated by the Oak Ridge National Laboratory.