Crystalline phase discriminating neutron tomography using advanced reconstruction methods

The measurement model in ToF neutron CT is given by the Beer–Lambert law, which relates the attenuation properties of the material to the measured intensity:

$$\begin{equation} I = I_0~\exp \left( -\int_L \mu(x, \lambda) \mathrm{d}x \right), \end{equation} \tag{ 1 }$$

where I₀ and I corresponds to the beam intensity incident on the object and on the detector element, respectively, L is a linear path through the object and µ(x, λ) is the wavelength-dependent attenuation coefficient at the physical position x in the object for the given wavelength λ. The probability of neutron-matter interaction is a function of neutron energy (wavelength) and is given by the microscopic total cross-section $\sigma_{\textrm{tot}} (\lambda),~(\mathrm{cm}^2)$ of the nucleus. When neutrons travel through material, the probability of interaction depends not only on the microscopic total cross-section $\sigma_{\textrm{tot}} (\lambda)$, but also on the number $N,~(\mathrm{atoms\,cm}^{-3})$ of nuclei within a unit volume of material. The macroscopic total cross-section $\Sigma_{\textrm{tot}} (\lambda),~(\mathrm{cm}^{-1})$ defines the probability of neutron-matter interaction per unit distance travelled of the neutron [29], i.e. the attenuation coefficient:

$$\begin{equation} \mu(\lambda) = \Sigma_{\textrm{tot}} (\lambda) = N \sigma_{\textrm{tot}} (\lambda). \end{equation} \tag{ 2 }$$

The microscopic total neutron cross-section $\sigma_{\textrm{tot}} (\lambda)$ is a linear combination of several contributions, defining the probability of elastic, inelastic, coherent and incoherent scattering, and absorption [30]. All interaction processes contribute to the decrease in transmitted intensity, however only coherent elastic scattering is responsible for characteristic abrupt increases in the transmitted intensity at 2d_hkl allowing the detection of specific lattice planes. In figure 1 we show the wavelength-dependent macroscopic cross-section $\Sigma_{\textrm{tot}} (\lambda)$ for materials employed in this study.

Zoom In Zoom Out Reset image size

For this study, a sample was constructed comprising six 6.3 mm diameter thin-walled cylindrical containers formed from aluminium foil. Five were filled with metal powders, namely copper (Cu), aluminium (Al), zinc (Zn), iron (Fe) and nickel (Ni); one cylinder was left empty. The containers were sealed and affixed around a hollow aluminium cylinder in a hexagonal close packed arrangement (figure 2). Powders were chosen for the sample to reduce the effect of texture which can strongly affect the shape of the Bragg edges [32]. This provides idealised Bragg edge spectra likely to be in good agreement with theoretical calculations of neutron transmission (figure 1).

Zoom In Zoom Out Reset image size

The data [33] was acquired at the imaging and materials science & engineering (IMAT) beamline operating at the ISIS spallation neutron source (Rutherford Appleton Laboratory, UK) [34, 35]. IMAT operates in a ToF measurement mode. Neutrons generated through spallation are slowed down by the L-H2 moderator, so they become 'thermalised'. Neutrons reach the sample and then the detector at different times according to their energy. Neutrons with shorter wavelength have higher velocity and are recorded first followed by increasingly less energetic neutrons. Consequently, given the distance travelled between the source and the detector, and the elapsed time from the pulse leaving the source the energy, and hence wavelength, can be calculated based on the de Broglie equation.

IMAT is installed on ISIS TS2 (2nd target station) which operates at 10 Hz repetition rate and the flight path between the source (more specifically, the neutron moderator) and the detector is ≈56.4 m giving an effective wavelength range around 6 Å. The IMAT beamline is equipped with a borated-microchannel plate (MCP) detector combined with Timepix chip [18]. The detector consists of a 2 × 2 array of 256 × 256 readout chips, resulting in 512 × 512 active pixels. The pixel size is 0.055 mm giving a field of view of approximately $28 \times 28~\mathrm{mm}^2$. The ToF spectrum recorded at each pixel can have up to 3100 time (energy) bins. There are small gaps and a slight misalignment between the readout chips [18], but these are not expected to have any significant implications for the current study.

The MCP detector functions in event timing mode, in which the arrival time of each neutron at each pixel is measured with respect to some external trigger (pulse). Unfortunately, the detector can register only one event per pixel per time frame. Consequently, slower (lower energy) neutrons have a lower probability of being detected, as some pixels are already occupied by faster neutrons. This effect is typically referred to as detector dead time. It is possible to set-up several time frames (also referred to as shutter intervals) within one pulse with arbitrary length and bin width for each time frame, with data read out in the end of each shutter interval. Data readout takes 320 µs and introduces gaps in the measured ToF signal. The introduction of several shutter intervals within one pulse reduces the spectral signal distortions due to detector dead time, but cannot fully eliminate it. In [36] the authors proposed an algorithm called 'overlap correction' to compensate for the counts lost. The efficiency of the algorithm was confirmed in [37]. For the present study, the detector shutter intervals were selected in such a way that data readout take place between theoretical Bragg edges for all measured materials (table 1). With this configuration 2843 energy channels were measured for every projection having wavelength resolutions between 0.7184 × 10⁻³ and 2.8737 × 10⁻³ Å.

The sample (figure 2) was placed and clamped onto the rotary stage. The sample was moved as close as possible to the MCP detector using the kinematic system. An alignment laser was used to visually align the vertical axis of the sample with respect to the vertical edge of the detector and to ensure that the object was within the field of view for all rotation angles. A set of spectral projections were acquired at 120 equally-spaced angular positions over 180^∘ rotation (1.5^∘ angular increments). Each projection was acquired with 15 min exposure time.

We acquired two sets of flat field (also referred to as open beam) images (before and after acquisition) to compensate for detector imperfections and decrease of beam intensity over time which was observed in other experiments. Also, flat field averaging is generally recommended to reduce ring artifacts. Therefore four spectral flat field images were acquired before, and four after, sample acquisition (eight in total), each one with 15 min exposure.The MCP detector has very low dark current noise [38]. Therefore dark field correction is not needed [34].

Following [36, 39], MCP detector related corrections (including both overlap correction and scaling to the same number of incident neutrons) were performed both for the set of 120 projections and the set of eight flat field images. We performed three different types of flat-field correction:

• Using an average of four flat field images acquired before sample data acquisition.
• Using two averaged flat field images, where the first flat field image corresponds to an average of four flat field images acquired before the data acquisition and the second one corresponds to an average of four flat field images acquired afterwards. Then, the actual flat field image for every projection is calculated based on pixel-wise channel-wise linear interpolation of intensity values between two averaged flat field images.
• Using an average of all eight flat field images. The intensity of the averaged flat field image is then multiplied by the quotient of the mean intensity in unoccupied detector columns in a projection image and mean intensity in the averaged flat field image itself. The scaling factor is calculated individually for each projection and each channel. We refer to this approach as flux normalisation.

To demonstrate the effect of the three different flat field correction approaches we show a flat-field corrected white beam (sum of all wavelength) transmission sinogram for a single slice of the acquired dataset (figure 3). Noticeable horizontal stripes are present for the first two flat field correction approaches which are eliminated with the flux normalisation. Profile lines (marked as a white solid line in the surrounding air in figure 3(a)) spotlight the difference between the first two approaches: linear interpolation compensates for the upwards trend in image intensity (figure 3(b)). Finally, the flux normalisation not only de-trends the signal but also drastically reduces the overall fluctuations. As a result, the flux normalisation approach was used in the present study.

Zoom In Zoom Out Reset image size

As a relatively short exposure time was used for the acquisition, spectral images were down-sampled in the spectral dimension to increase counting statistics and produce reasonable spectra. Down-sampling was performed for each shutter interval independently by taking the average of each 16, 8, 8 and 4 channels for the respective shutter intervals (table 1). As a result, each preprocessed projection contains 339 channels having a uniform bin width of 11.5 × 10⁻³ Å over a wavelength range between 1 and 5 Å.

In figure 4(a), we show the preprocessed attenuation sinograms for individual wavelength channels. Even after spectral downsampling, strong noise is present in the measured tomographic data. Recorded spectra (figure 4(b)) taken for a path through surrounding air (P1), iron (P2), both iron and nickel (P3) and copper (P4) show that noise dominates over valuable Bragg edge information and makes Bragg edge characterisation for individual pixels barely feasible.

Zoom In Zoom Out Reset image size

It is also important to recognise that the incident spectrum on IMAT is a function of the neutron moderation process. The incident spectrum has a crude 'bell shape' with a peak around 2.6 Å [34, 35]. The recorded spectrum is also affected by the detector dead time. As a result the number of counts, and consequently noise level, in each channel depends on wavelength and position of shutter intervals. The elevated noise level is noticeable in the wavelength channels below 1.5 Å and above 4.5 Å, where IMAT has the lowest flux, and between 2.5 and 3 Å where a combination of high flux and detector dead time causes distortions in the recorded spectrum. The wavelength-noise dependence has direct implications on the reconstruction quality of each individual channel.

Short of further downsampling of the recorded multi-channel images, advanced reconstruction methods are needed to handle noisy multi-channel neutron CT data. A potential solution to this low-count CT problem is to effectively exploit available prior information, such as the anticipated image properties and structural correlation between channels.

Tomographic reconstruction aims to recover a map of sample attributes from a set of integral measurements acquired at various angles. In every wavelength channel, Bragg edge neutron CT can be well approximated by the standard absorption tomography model (Radon transform). Consequently, algorithms developed for x-ray CT can be conveniently used to reconstruct the neutron attenuation map for every wavelength channel. Here, we focus on the two-dimensional reconstruction problem in the spatial domain; extension to the third spatial dimension is straightforward.

Given a multi-channel tomographic dataset with K channels, a conventional approach is to reconstruct each channel independently using the FBP algorithm. The FBP algorithm is derived from the Fourier Slice theorem which relates line integral measurements to two-dimensional Fourier transform of an object's slice. In FBP-type reconstruction methods, projections are filtered independently and then back-projected onto the plane of the tomographic slice. FBP reconstruction is very fast but requires high quality input data and dense angular sampling to achieve good results.

Alternatively the inverse problems framework can be used to reconstruct low-count CT data. Consider a single channel k, k = 0, 1..., K − 1, in a spectral dataset consisting of K wavelength channels, and let $\overline{b}(\lambda_k)$ be a recorded discrete sinogram with P × D elements, with P being the total number of projections and D being the number of pixels in a detector row. The sinogram $\overline{b}(\lambda_k)$ is vectorised as a column vector with M = PD elements (obtained by stacking the columns in the original two-dimensional sinogram). Let $\overline{u}(\lambda_k)$ be a vectorised two-dimensional array of material attributes we want to reconstruct with N = D² elements (voxels). The discrete version of the Radon transform is then given by the projection operator $\overline{A}$ containing M × N elements. If i, i = 0, 1, ..., M − 1 and j, j = 0, 1, ..., N − 1, then $\overline{A}_{(i,j)}$ is the length of intersection of the ith ray with the jth voxel. The reconstruction problem then takes the form:

$$\begin{equation} \overline{b}(\lambda_k) = -\textrm{ln} \left( \frac{I(\lambda_k)}{I_0(\lambda_k)} \right) \approx \overline{A} \overline{u}(\lambda_k), \end{equation} \tag{ 3 }$$

where $I_0(\lambda_k)$ and I(λ_k ) corresponds to the wavelength-dependent flux measured without (open beam), and with, a sample in the field of view, respectively.

Data acquisition in Bragg edge neutron CT is very time consuming since neutron fluxes are typically low (compared to synchrotron x-ray sources) and because the detected neutrons are shared among multiple ToF (energy) channels. Hence, the acquisition mode does not inherently provide sufficient data for a numerically stable solution of equation (3), similar to the low-dose tomography problem in medical x-ray CT imaging. The resulting problem is said to be ill-posed in a mathematical sense, i.e. the naive solution of equation (3) will barely produce any useful result. A common treatment for ill-posed problems is to design a surrogate problem which is consistent with the original problem, but is well-posed and computationally tractable. In other words, we obtain a stable approximate solution by incorporating some prior knowledge about the problem.

Multi-channel CT data can be considered as a stack of two- or three-dimensional datasets, where every voxel in the reconstructed volume contains a vector with a spectral material response, spanning an additional dimension. In this respect, every data point in the material response vector belongs in the same physical structure, i.e. channels share structural information. Then, inter-channel correlations can be utilised to improve reconstruction quality. The spectral CT reconstruction takes a form:

$$\begin{equation} b = A u, \end{equation} \tag{ 4 }$$

where b and u are obtained by stacking K column vectors $\overline{b}(\lambda_k)$ and $\overline{u}(\lambda_k)$, respectively, and $A = I_{(K \times K)} \otimes \overline{A}$, $\otimes$ is the Kronecker product, and $I_{(K \times K)}$ is the identity matrix of order K.

The reconstruction problem is then constructed as an optimisation problem:

$$\begin{equation} \mathop{\textrm{arg}}{\textrm{min}}_u \{F(u) = f(b,Au) + \alpha g(u) \}, \end{equation} \tag{ 5 }$$

where f(b, Au) is a data fidelity metric which measures the discrepancy between the projection of solution u and the acquired data b. The regularisation term g(u) imposes certain prior assumptions typically expressed in terms of desired image properties. The scalar parameter α controls the trade-off between the fit to the acquired data u and the regularisation. The choice of regulariser depends on the underlying problem since there is no unique regulariser which performs best in all problems. The art of choosing the right regulariser and finding a balance between the two terms is a challenging task, as regularisation inevitably introduces bias into the solution. The bias is a price one pays for solving a problem which is not solvable by other means [40].

TV is the most commonly used regulariser. TV encourages a sparse image gradient and consequently favours piece-wise constant images with sharp boundaries [19, 20]. The TV model describes well the spatial dimension of images of the type to be reconstructed in this study, but it is known to introduce 'staircase' artifacts for piece-wise affine signals, i.e. a ramp or features similar to the ones observed in wavelength dependent attenuation coefficient (figure 1) might be reconstructed as piece-wise constant, reminiscent of a staircase. TV is also known to suffer from intensity reductions and, as a result, low contrast regions might be lost. Multiple modifications of TV have been proposed to cure the above shortcomings.

The TNV approach proposed in [25] explicitly promotes reconstructions with common edges across all channels. The idea is based on the fact that having shared gradient directions is equivalent to have a rank-one Jacobian of a multi-channel image. Therefore, TNV penalises the singular values of the Jacobian. In the spatial dimension, TNV has similar properties to TV regularisation, i.e. it also favours a sparse image gradient. Consequently, TNV correlates channels and improves reconstruction quality by promoting common structures in multichannel images. Application of TNV for reconstruction of multi-channel images has been successfully demonstrated in [25, 26, 41]. The reconstruction problem is then formulated as,

$$\begin{equation} F(u) = \| Au-b \|^2_2 + \alpha \mathrm{TNV}(u). \end{equation} \tag{ 6 }$$

Similar to TV, TNV suffers from a loss of contrast. Secondly, TNV does not allow the decoupling of regularisation parameters for the spatial and spectral dimensions which makes it impossible to balance the level of regularisation between dimensions.

Here we propose a novel tailored regulariser which treats the spatial and spectral dimensions separately. As the TV model captures piece-wise image properties in the spatial dimension, we rely on another regulariser to support reconstruction in the spectral dimension. A natural remedy for the staircase artifact is to use higher order derivatives in the regularisation term. Here, we use TGV [21] which allows for both sharp changes in spectral signal and gradual intensity changes. TGV also effectively exploits inter-channel correlation. In this case the reconstruction problem is formulated as,

$$\begin{equation} F(u) = \| Au-b \|^2_2 + \beta \mathrm{TV}_{x,y} (u) + \gamma \mathrm{TGV}_c (u). \end{equation} \tag{ 7 }$$

Here we use $\mathrm{TV}_{x,y}$ to designate a TV operator over two spatial dimensions x and y, whereas $\mathrm{TGV}_c$ operates over the channel dimension. $\mathrm{TV}_{x,y}$ is a channel-by-channel operator applied to all channels individually and then summed over all channels. $\mathrm{TGV}_c$ is a one-dimensional operator applied in each individual voxel across all channels simultaneously.

TGV regularisation has been already successfully applied to multi-channel tomographic data [42–46] including time resolved magnetic resonance imaging, energy-dispersive x-ray spectroscopy and high-angle annular dark field data. To the best of our knowledge, this is the first attempt to combine TV and TGV for spectral tomographic problems, in this case Bragg edge neutron CT data. The decoupling of regularisation into the spatial and spectral dimension allows us to balance the two regularisation terms appropriately and promote the desired properties in both dimensions independently. However, this greater flexibility comes at higher computational cost. The more terms that are included into the optimisation problem equation (5), the less it is suitable for parallelisation and software acceleration. Unfortunately, currently we cannot give a quantitative indicator here because our current implementation of the regularization methods described in the paper was not optimized from a performance point of view. We have quite significantly accelerated the TNV method but TV-TGV is still in the pipeline. Therefore a direct comparison in terms of execution times would give a biased estimate representing our software optimization efforts rather than provide a measure of the performance of the method itself.

A common way to solve the large-scale optimisation problem in equation (5) is by means of the fast iterative shrinkage-thresholding algorithm (FISTA) [47] or the primal-dual hybrid gradient (PDHG) algorithm [48]. Here we used the implementation of both FISTA and PDHG in the CCPi CIL [27, 28]. CIL provides a highly modular Python library for prototyping of reconstruction methods for multi-channel data such as spectral and dynamic (time-resolved) CT. CIL wraps a number of third-party libraries with hardware-accelerated building blocks for advanced tomographic algorithms. CIL relies on the ASTRA toolbox [49–51] to perform forward- and back-projection operations and provides a set of various regularisers through the CCPi Regularisation Toolkit [52]. TNV was solved using FISTA employing the CIL plugin for the CCPi Regularisation Toolkit [52] for the TNV proximal operator. TV-TGV was implemented directly in CIL and solved using PDHG.

Regularisation parameters were carefully chosen to achieve both noise suppression and feature preservation in both spatial and spectral dimensions. For TNV reconstruction, the regularisation parameter α was set to 0.01, while TV-TGV parameters were set to β = 0.0075 and γ = 0.3. In both cases, we ran 2000 iterations of the reconstruction algorithm.

Quantitative analysis of the reconstructed spectra commonly includes detection and characterisation of the Bragg edges. The Bragg edge positions (in terms of d-spacing) allows compositional mapping, as each crystalline structure will have a unique set of lattice spacings and hence fingerprint in the neutron transmitted spectrum. The shape of the detected Bragg edges, i.e. deviation from abrupt step-like theoretical edge, supports characterisation of crystallographic properties.

In [39] authors developed a preprocessing and Bragg edge analysis tool called BEAn for wavelength-resolved neutron transmission data. We adopted Bragg edge detection and fitting functionality from BEAn for the present study. The resulting Bragg edge fitting routine consists of the following steps. (1) BEAn has been developed for transmission data, therefore prior to the analysis we convert reconstructed attenuation data to transmission. (2) The spectrum is then smoothed using the Savitzky–Golay filter, differentiated, smoothed again and passed to a peak-finding routine to detect possible edge locations. Note, smoothing is performed only to support the peak-finding routine and all further steps are performed on unsmoothed spectra. (3) For each detected Bragg edge, a model function derived in [53] is fitted using the non-linear least squares fitting method. Since the model function is highly non-linear, the fitting procedure is very sensitive to initialisation parameters. Therefore we use a brute force approach to help the fitting procedure to achieve acceptable fit. That is, initialisation parameters are drawn from uniform random distributions chosen to cover a realistic range for each parameter in the model function. The best fit is chosen based on the minimum root mean squared error between the fitted function and experimental data.

Additionally, the spatial distribution of individual materials in a sample can be obtained by decomposing the reconstructed spectral images into individual material maps. Material decomposition in spectral CT is an ill-posed problem by its own and various approaches have been proposed in the literature, mainly for spectral x-ray CT. The first group of methods performs material decomposition in the sinogram domain, following by reconstruction of material specific sinograms [54, 55]. The second group of methods employs simultaneous reconstruction and material decomposition [56]. A third approach is to perform material decomposition in the image domain, i.e. after the reconstruction step [57, 58]. Limitations and merits of the individual methods are beyond the scope of this research. As a proof of principle, we adapt the later approach here. In particular, we use the so called volume conservation principle [59, 60], where each voxel in the obtained material maps corresponds to a dimensionless volume fraction occupied by the corresponding material with the full voxel corresponding to a unit volume. This is a valid assumption here as the materials employed in the current study do not mix. Under the volume conservation assumption, the voxel-wise sum of all material maps has to be equal to 1 in each voxel.

Let $\hat{u}$ be a K × D² matrix version of u and let M be the material basis K × S matrix built up from the predicted wavelength-dependent neutron spectra (figure 1), with S = 6 being the number of materials employed in this study (five metals + air). Then, volume fractions of every material v, consisting of S row vectors with D² elements, can be obtained by solving:

$$\begin{equation} \hat{u} \approx Mv. \end{equation} \tag{ 8 }$$

To solve equation (8), we formulate the following optimisation problem:

$$\begin{align} \min_{v} \quad & \| \hat{u} - Mv \|^2_F \nonumber\\ \textrm{s.t.} \quad & v \geq 0 \nonumber\\ & \mathbf{1}_S^T v = \mathbf{1}_{D^2}^T \end{align} \tag{ 9 }$$

where $\| \cdot \|^2_F$ is the Frobenius norm and $\mathbf{1}^T$ is a row vector of all ones of appropriate dimension. The first constraint is an element-wise non-negativity constraint enforcing the volume fractions v to be greater than or equal to 0; the second constraint expresses the volume preservation principle. We use CVX [61, 62] with the MOSEK solver to solve this optimisation problem.

4.1.1. White beam reconstruction.

In order to give a visual reference for the spatial dimension, we first perform conventional FBP reconstruction (implemented as an eponymous processor in CIL [27]) of the white beam sinogram (sum of all wavelength channels), as shown in figure 5. All cylinders are seen clearly defined and with different contrast for each element. The copper powder used in this study is known to have an average particle size comparable to the voxel size. Therefore reconstruction of the copper cylinder appears inhomogeneous, i.e. some reconstructed voxels contain a mixture of copper and air, while some are fully occupied by larger copper particles. There are also noticeable roundness deviations in the reconstructed cross-section of the cylinders. These deviations are caused by the simple fabrication process used for making the phantom. There are also small 'bumps' where a small portion of powder penetrated overlapping foil layers.

Zoom In Zoom Out Reset image size

4.1.2. The spatial dimension.

Figure 6 shows two-dimensional slices for selected (individual) wavelength channels reconstructed using the conventional FBP approach and using the regularised iterative methods discussed in this paper. Linescans across the reconstructed cylinders in figure 7 offer a detailed comparison between the reconstruction methods. We only perform two-dimensional reconstructions in the spatial dimension in this study but results can be generalised to the third dimension. We have chosen a slice roughly in the middle of the upper half of the detector plate in order to avoid detector rows near to the top of the detector or the gap between the detector read-out chips (section 2.3).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

As expected, at these signal levels the FBP reconstruction is barely interpretable (figure 6, top row). Both TNV and TV-TGV demonstrate drastic improvement in reconstruction quality and noise suppression (figures 6 and 7). A region-of-interest between two cylinders highlights significant smearing of features in the TNV reconstruction (figure 6, middle row); the same features appear sharper in the TV-TGV reconstruction (figure 6, bottom row). TNV suffers from the same limitations as conventional TV-based regularisation: it might oversmooth images. Therefore, TNV suppresses noise and ring artifacts visible in FBP reconstruction but produces blurred and enlarged rings especially prominent in shorter and longer wavelength channels, where counts are much lower. Aluminium has very low neutron attenuation and is invisible in the TNV reconstruction due to contrast loss; another known drawback of both TV and TNV regularisation methods.

The results demonstrate that the TV-TGV reconstruction is more effective for low contrast features. Furthermore, the TV-TGV reconstruction does not suffer from ring artefacts and the Al cylinder is visible in the reconstructed slices and profile lines (figure 7, bottom row). Fine features inside the copper cylinder are also partially preserved in the TV-TGV reconstruction (figure 7, top row).

4.1.3. The spectral dimension.

The ground truth for the spectral dimension is taken to be the predicted neutron cross-section for the selected materials (figure 1). Individual spectra reconstructed for one $0.055^3~\mathrm{mm}^3$ voxel in each of the five materials are plotted in figure 8 alongside the theoretical predictions. The voxel locations inside the cylinders were chosen arbitrarily (marked in figure 6).

Zoom In Zoom Out Reset image size

The limited counting statistics acquired in our experiment make the FBP reconstructed spectra uninterpretable. By contrast, both the TNV and TV-TGV reconstructions show drastic improvements in reconstruction quality. The reconstructed spectra for Fe, Ni and Cu closely follow the Bragg edge spectra for both TNV and TV-TGV reconstructions. For the low-attenuation materials (Zn and Al), the TV-TGV clearly outperforms TNV. The amplified noise visible in the TNV reconstructions between 2.5 and 3 Å is caused by an increase in noise level in the input data, as was highlighted in section 2.5.

For the high-attenuation materials (Fe, Ni, Cu) both TNV and TV-TGV show comparable performance. While the TV-TGV produces a much smoother spectra, the Bragg edges appear to be less prominent due to smearing (for instance, small edges around 2 Å in Fe). In the case of the TNV regularisation, the noise dominates over smaller Bragg edges. For the low-attenuation materials (Al and Zn), TV-TGV shows better performance as some Bragg edges are visible in the reconstructed spectrum, but are lost in the noise in the case of TNV.

4.2.1. Mapping of crystallographic information.

Mapping of crystallographic information was performed on individual spectra reconstructed for one $0.055^3~\mathrm{mm}^3$ voxel in each of the five materials (the same voxels as in the previous section, marked in figure 6). Following the procedure outlined in section 3.3, the Bragg edges were detected and characterised based on fitting of a dedicated model function (figure 9). Fitting was performed individually around each detected Bragg edge. The fitted function is shown as the orange solid line. Since the locations of the Bragg edges are directly related to interplanar spacing d_hkl , we use the reference $\lambda_{hkl} = 2d_{hkl}$ for the corresponding materials to assess the mapping of the crystalline structural information (table 2). We also show an absolute error $\varepsilon = \| \lambda_{hkl} - \lambda^{^{\prime}}_{hkl}\|$ in the estimated d-spacing $\lambda^{^{\prime}}_{hkl}$. For some Bragg edges, the edge detection routine identified edges, but the least squares fitting procedure did not converge to a good solution. We use '(-)' to designate such edges.

Zoom In Zoom Out Reset image size

Table 2. Results of mapping of crystallographic information. Estimated $\lambda^{^{\prime}}_{hkl}$ is compared vs. reference $\lambda^{^{\prime}}_{hkl} = 2 d_{hkl}$, where d_hkl were extracted from NXS code [31]. For some Bragg edges, the edge detection routine identified edges but the least squares fitting procedure did not converge to a good solution. We use '(−)' to designate such edges. Edges that could not be detected by the automated procedure are marked with '–'.

		TNV		TV-TGV
	λhkl	$\lambda^{^{\prime}}_{hkl}$	ε	$\lambda^{^{\prime}}_{hkl}$	ε
Fe(110)	4.054	4.055	0.001	4.077	0.023
Fe(200)	2.866	2.819	0.047	2.853	0.013
Fe(211)	2.340	2.366	0.026	2.332	0.008
Fe(220)	2.026	–		(−)
Fe(310)	1.812	(−)		–
Fe(222)	1.654	–		–
Fe(321)	1.532	–		1.492	0.04
Ni(111)	4.064	4.073	0.009	4.076	0.012
Ni(200)	3.520	3.533	0.013	3.545	0.025
Ni(220)	2.490	2.522	0.032	2.512	0.022
Ni(311)	2.122	2.122	0.000	2.136	0.014
Ni(222)	2.032	–		–
Ni(400)	1.760	–		–
Ni(331)	1.616	(−)		1.609	0.007
Cu(111)	4.168	4.216	0.048	4.151	0.017
Cu(200)	3.610	(−)		3.662	0.052
Cu(220)	2.552	(−)		2.627	0.075
Cu(311)	2.176	2.163	0.013	2.215	0.039
Cu(222)	2.084	–		–
Cu(400)	1.804	–		–
Cu(331)	1.656	–		1.629	0.027
Al(111)	4.676	–		4.729	0.053
Al(200)	4.050	–		4.040	0.01
Al(220)	2.864	–		(−)
Al(311)	2.442	–		2.435	0.007
Al(222)	2.338	–		–
Al(400)	2.024	–		–
Al(331)	1.858	–		1.982	0.124
Zn(002)	4.950	–		–
Zn(100)	4.608	(−)		4.648	0.040
Zn(101)	4.178	(−)		4.209	0.031
Zn(102)	3.372	(−)		3.391	0.019
Zn(103)	2.682	(−)		2.671	0.011

For the high-attenuation materials in this study, i.e. Fe and Ni, both TNV and TV-TGV reconstructions show comparable results in terms of error and uncertainty in estimated interplanar spacing. The lower the material's attenuation coefficient, the more the TV-TGV reconstruction surpasses the TNV results. Thus, in the TNV reconstruction we can estimate only some interplanar distances; no Bragg edges can be detected and characterised for Al. Remarkably, the TV-TGV reconstruction is able to locate some Bragg edges even for the low-attenuation materials such as Zn and Al. Accuracy in estimated interplanar distances deteriorate with prominence of corresponding Bragg edges. Also, Bragg edge fitting performs better for isolated edges because of edges smearing.

Of course one could increase the sensitivity to the low scattering materials by increasing the acquisition period or the number of projections, but only at the expense of longer total experiment time, which may be prohibitive. In this study we also used an automated procedure to detect and fit the model function (section 3.3). A more careful manual procedure might improve the Bragg edge fitting results, especially for the noisier TNV reconstruction.

4.2.2. Material maps.

Finally, in figure 10 we show spatial material maps calculated on the basis of the material decomposition method described in section 3.3. The difference between the material maps calculated from TNV and TV-TGV reconstructed images is quite subtle, especially for high-attenuation materials. Advantages of the TV-TGV approach are most pronounced in the Al material map. Compared with TNV, TV-TGV does a better job in preserving the fine Al foil containers and the Al cylinder appears more uniform. The material maps for both reconstructions exhibit some artifacts, for instance, the faint Fe cylinder is visible in the Ni map and similarly Zn is visible in the Cu map. The artifacts are caused by the fact that the most prominent Bragg edges almost coincides in the Fe and Ni spectra, as well as in the Cu and Zn spectra (figure 1). In this study we used a fairly generic method for material decomposition. As we have already mentioned in section 3.3, the material decomposition problem is an ill-posed problem and various regularisation methods can be used to improve material decomposition depending on underlying image properties in the spatial and spectral dimensions.

Zoom In Zoom Out Reset image size