Crystalline phase discriminating neutron tomography using advanced reconstruction methods

Time-of-flight neutron imaging offers complementary attenuation contrast to X-ray computed tomography (CT), coupled with the ability to extract additional information from the variation in attenuation as a function of neutron energy (time of flight) at every point (voxel) in the image. In particular Bragg edge positions provide crystallographic information and therefore enable the identification of crystalline phases directly. Here we demonstrate Bragg edge tomography with high spatial and spectral resolution. We propose a new iterative tomographic reconstruction method with a tailored regularisation term to achieve high quality reconstruction from low-count data, where conventional filtered back-projection (FBP) fails. The regularisation acts in a separated mode for spatial and spectral dimensions and favours characteristic piece-wise constant and piece-wise smooth behaviour in the respective dimensions. The proposed method is compared against FBP and a state-of-the-art regulariser for multi-channel tomography on a multi-material phantom. The proposed new regulariser which accommodates specific image properties outperforms both conventional and state-of-the-art methods and therefore facilitates Bragg edge fitting at the voxel level. The proposed method requires significantly shorter exposure to retrieve features of interest. This in turn facilitates more efficient usage of expensive neutron beamline time and enables the full utilisation of state-of-the-art high resolution detectors.


Introduction
Neutron radiography [1] and later neutron computed tomography (CT) [2,3] have been available at neutron sources for some time. As uncharged particles, neutrons probe the nucleus rather than the electron cloud and so in contrast to X-ray CT [4] where the attenuation contrast increases with atomic number, neutron attenuation can give strong contrast between neighbouring elements in the periodic table (e.g. Cu and Ni) or even different isotopes [5]. In addition to conventional attenuation contrast, the transmitted beam contains important information regarding coherent scattering from the crystalline lattice. Indeed, for many materials the scattering cross-section is comparable to the attenuation. Since the thermal neutrons produced by neutron spallation sources have wavelengths comparable to interplanar lattice spacings, polycrystalline materials will scatter the incident beam elastically according to Bragg's law. Since Bragg scattering can occur only for wavelengths shorter than twice the spacing between the lattice planes (d hkl ), the transmitted neutron spectrum exhibits characteristic abrupt increases in the transmitted intensity at these wavelengths. As a result the transmitted spectrum has a characteristic signature displaying distinct jumps in transmitted intensity corresponding to the Bragg edges for all the crystallographic phases in the material.
For pulsed spallation neutron sources it is possible to infer the energy and hence wavelength of each detected neutron from its time of flight (ToF) from the source to the detector. The development of pixelated ToF detectors enabled the Bragg edges to be imaged [6,7]. The relation between spectral fingerprint and crystalline properties allows identification of polycrystalline materials and characterisation of their properties [7,8,9] such as phase [10,11], texture and strain [6,12,13,14]. Combined with sample rotation, three-dimensional Bragg edge neutron CT is a natural extension of two-dimensional Bragg edge neutron imaging, as has been already demonstrated in several studies [15,16,17]. However, as Kockleman et al. [8] pointed out in their early preliminary neutron energy selective imaging experiments, the signal, which is usually integrated over the white beam, becomes particularly low when distributed across many hundreds of ToF (energy) channels making tomographic imaging slow and noisy. In recent years the size of pixelated ToF discriminating detectors has improved significantly (from a 10 × 10 array of 2 × 2 mm 2 [6] to a 512 × 512 array of 0.055 × 0.055 mm 2 [18]), thereby improving the potential resolution, but further decreasing the flux received by each pixel and further confounding the reconstruction challenge.
The conventional way of reconstructing tomographic datasets is filtered back-projection (FBP) which is a fast and well-established method, but demanding in terms of input data quality. As a result of the low count rates and hence slow acquisition, energy dispersive spectra tend to be heavily down-sampled prior to reconstruction and, even with reduced resolution, averaging over a relatively large region-of-interest in reconstructed images is still required to improve signal-to-noise ratio and characterise Bragg edges. Artificially induced low resolution thus currently hinders application of the technique.
Here we present advanced reconstruction methods that make energy-dispersive neutron tomography a practical proposition. We explore the application of dedicated multi-channel reconstruction techniques with inter-channel correlation to improve reconstruction quality without compromising valuable resolution. Contrary to previous work [16,17] where each channel is treated as a separate reconstruction problem, we use iterative methods with prior information to jointly reconstruct spectral images.
Iterative methods formulate reconstruction as an optimisation problem consisting of a data fidelity term and a regulariser. The latter term encodes prior knowledge about the image and helps to find a unique solution for the ill-posed tomographic problem by encouraging desired properties in reconstructed images. There is no unique regulariser which performs well for all tomographic problems and it has to be carefully selected depending on the underlying tomographic data properties. In this study we tailored a regulariser specifically for Bragg edge neutron CT based on a combination of Total Variation (TV) regularisation [19,20] in the spatial dimension and Total Generalised Variation (TGV) regularisation [21,22] in the spectral dimension. TV preserves edges and suppresses noise by encouraging piece-wise constant regions in reconstructed spatial images [19,20], while TGV regularisation promotes characteristic piece-wise smooth behaviour in the spectral dimension [21]. For brevity, henceforth we refer to the proposed method as TV-TGV.
We compare TV-TGV against conventional FBP reconstruction and a state-of-the-art regulariser for multi-channel CT images -Total Nuclear Variation (TNV). TNV is a recent regulariser which enforces common edges across all channels in multi-channel images [23,24,25]. Wider implementation of these advanced reconstruction methods is possible through the open source CCPi Core Imaging Library (CIL) reconstruction framework [26,27], for example one may employ TV only or TGV only if desired. The performance of all methods is demonstrated on a multi-material sample comprising aluminium cylinders filled with various metallic powders of high purity.
Experimental details with a particular focus on preprocessing are presented in section 2. In sec-tion 3, we provide a general overview of reconstruction methods used in the present study. In section 4, a comparison between all the reconstruction methods in this study is presented. A Bragg edge fitting procedure is further used to detect and locate Bragg edges in TNV and TV-TGV reconstructions and extract crystallographic information. We also demonstrate decomposition of reconstructed spectral images into individual material maps. Discussion and conclusions are given in section 5.

Measurement model
The measurement model in ToF neutron CT is given by the Beer-Lambert law, which relates attenuation properties of material to measured intensity: where I 0 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 σ tot (λ), [cm 2 ] of the nucleus. When neutrons travel through material, the probability of interaction depends not only on the microscopic total cross-section σ tot (λ), but also on the number N, [atoms/cm 3 ] of nuclei within unit volume of material. The macroscopic total cross-section Σ tot (λ), [cm −1 ] defines the probability of neutron-matter interaction per unit distance travelled of the neutron [28], i.e. the attenuation coefficient: The microscopic total neutron cross-section σ tot (λ) is a linear combination of several contributions, defining the probability of elastic, inelastic, coherent and incoherent scattering, and absorption [29]. All interaction processes contribute to the decrease of 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 fig. 1 we show the wavelength-dependent macroscopic cross-section Σ tot (λ) for materials employed in this study.

Sample design
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 ( fig. 2). Powders were chosen for the sample to reduce the effect of texture which can strongly affect the shape of the Bragg edges [31]. This provides idealised Bragg edge spectra likely to be in good agreement with theoretical calculations of neutron transmission ( fig. 1).

Instrument settings
The data was acquired at the Imaging and Materials Science & Engineering (IMAT) beamline operating at the ISIS spallation neutron source (Rutherford Appleton Laboratory, U.K.) [32,33]. IMAT operates in a ToF measurement mode. Neutrons generated through spallation are slowed down by the L-H2 moderator on IMAT, 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 : Theoretical neutron spectra for materials employed in this study (Calculated using the NXS software package [30]).
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, neutron moderator) and the detector is ≈ 56.4 m giving an effective wavelength range around 6 Å. The IMAT beamline is equipped with a boratedmicrochannel 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 × 28 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 [34] authors proposed an algorithm called "overlap correction" to compensate for the counts lost. The efficiency of the algorithm was confirmed in [35]. 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 −3 Å and 2.8737 · 10 −3 Å.

Data acquisition
The sample ( fig. 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 4 spectral flat field images were acquired before and 4 after sample acquisition (8 in total), each one with 15 min exposure. The MCP detector has no dark current noise. The raw data is available at [36].

Data preprocessing
Following [34,37], 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 8 flat field images. We performed three different types of flat-field correction: • Using an average of 4 flat field images acquired before sample data acquisition.
• Using 2 averaged flat field images, where the first flat field image corresponds to an average of 4 flat field images acquired before the data acquisition and the second one corresponds to an average of 4 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 2 averaged flat field images.
• Using an average of all 8 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 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 ( fig. 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 fig. 3a) spotlight the difference between the first two approaches: linear interpolation compensates for the upwards trend in image intensity ( fig. 3b). 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.
As 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 −3 Å over a wavelength range between 1 Å and 5 Å.
In fig. 4a, 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 ( fig. 4b) 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.
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 Å [32,33]. 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
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 b(λ 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 b(λ k ) is vectorised as a column vector with M = P D elements (obtained by stacking the columns in the original twodimensional sinogram). Let u(λ k ) be a vectorised two-dimensional array of material attributes we want to reconstruct with N = D 2 elements (voxels). The discrete version of the Radon transform is then given by the projection operator A containing M × N elements. If i, i = 0, 1, . . . M − 1 and j, j = 0, 1, . . . , N − 1, then A (i,j) is the length of intersection of the i.th ray with the j.th voxel. The reconstruction problem then takes a form: where I 0 (λ 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 time of flight (energy) channels. Hence, the acquisition mode does not inherently provide sufficient data for a numerically stable solution of eq. (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 eq. (3) will barely produce any useful result. A common treatment for illposed 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: where b and u are obtained by stacking K column vectors b(λ k ) and u(λ k ), respectively, and A = I (K×K) ⊗ A, ⊗ is the Kronecker product, and I (K×K) is the identity matrix of order K. The reconstruction problem is then constructed as an optimisation problem: 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 fit with 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 [38]. 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 ( fig. 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 [24] 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 [24,25,39]. The reconstruction problem is then formulated as, 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, Here we use TV x,y to designate a TV operator over two spatial dimensions x and y, whereas TGV c operates over channel dimension. TV x,y is a channel-by-channel operator applied to all channels individually and then summed over all channels. 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 [40,41,42,43,44] including time resolved magnetic resonance imaging (MRI), energy-dispersive X-ray spectroscopy (EDXS) and high-angle annular dark field (HAADF) 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 eq. (5), the less it is suitable for parallelisation and software acceleration.

Numerical implementation
A common way to solve the large-scale optimisation problem in eq. (5) is by means of the Fast Iterative Shrinkage-Thresholding Algorithm (FISTA) [45] or the Primal-Dual Hybrid Gradient (PDHG) algorithm [46]. Here we used the implementation of both FISTA and PDHG in the CCPi Core Imaging Library (CIL) [26,27]. 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 [47,48,49] to perform forward-and back-projection operations and provides a set of various regularisers through the CCPi Regularisation Toolkit [50]. TNV was solved using FISTA employing the CIL plugin for the CCPi Regularisation Toolkit [50] 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.

Post-processing and analysis
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 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 [37] authors developed a preprocessing and Bragg edge analysis tool called BEAn for wavelengthresolved 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 [51] 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 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 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 [52,53]. The second group of methods employs simultaneous reconstruction and material decomposition [54]. A third approach is to perform material decomposition in the image domain, i.e. after the reconstruction step [55,56]. 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 [57,58], 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û be a K × D 2 matrix version of u and let M be the material basis K × S matrix built up from the predicted wavelength-dependent neutron spectra ( fig. 1), with S = 6 being the number of materials employed in this study (5 metals + air). Then, volume fractions of every material v, consisting of S row vectors with D 2 elements, can be obtained by solving: To solve eq. (8), we formulate the following optimisation problem: where · 2 F is the Frobenius norm and 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 [59,60] with the MOSEK solver to solve this optimisation problem.

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 [26]) of the white beam sinogram (sum of all wavelength channels), as shown in fig. 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 crosssection 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.  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 fig. 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).

The spatial dimension
As expected, at these signal levels the FBP reconstruction is barely interpretable ( fig. 6, top row). Both TNV and TV-TGV demonstrate drastic improvement in reconstruction quality and noise suppression ( figs. 6 and 7). A region-of-interest between two cylinders highlights significant smearing of features in the TNV reconstruction ( fig. 6, middle row); the same features appear sharper in the TV-TGV reconstruction ( fig. 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 ( fig. 7, bottom row). Fine features inside the copper cylinder are also partially preserved in the TV-TGV reconstruction ( fig. 7, top row).

The spectral dimension
The ground truth for the spectral dimension is taken to be the predicted neutron cross-section for the selected materials ( fig. 1). Individual spectra reconstructed for one 0.055 3 mm 3 voxel in each of the 5 materials are plotted in fig. 8 alongside the theoretical predictions. The voxel locations inside the cylinders were chosen arbitrarily (marked in fig. 6).
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  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.

Mapping of crystallographic information
Mapping of crystallographic information was performed on individual spectra reconstructed for one 0.055 3 mm 3 voxel in each of the 5 materials (the same voxels as in the previous section, marked in fig. 6). Following the procedure outlined in section 3.3, the Bragg edges were detected and characterised based on fitting of a dedicated model function ( fig. 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 λ hkl = 2d hkl for the corresponding materials to assess the mapping of the crystalline structural information (table 2). We also show an absolute error ε = λ hkl − λ hkl in the estimated d-spacing λ 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.
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 acqui- sition 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.

Material maps
Finally, in fig. 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 ( fig. 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.

Discussion and conclusions
The unique features of the pulsed time of flight neutron source allow quantification of material properties in bulk samples. Low flux and slow acquisition result in long scan times being required to record sufficient counts across hundreds of wavelength channels, making material characterisation at the full attainable resolution infeasible in practice. Joint reconstruction which exploits inter-channel correlations in multi-channel images is an effective way to shift feature recovery from exposure to reconstruction. In this work we have studied the performance of advanced reconstruction methods with the dedicated multi-channel regularisation techniques for Bragg edge neutron CT. We have demonstrated that the tailored TV-TGV regularisation technique, which favours specific image properties in the spatial and spectral dimensions, allows retrieval of crystallographic properties at a resolution previously unattainable through conventional reconstruction methods using the same exposure time.   [30]. 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 "-".
of low-contrast features. Furthermore, quantitative comparisons were also used to evaluate the performance of both methods. We extracted crystallographic properties from the reconstructed spectra based on detection and characterisation of Bragg edges. TV-TGV facilitated extraction of interplanar spacings for all materials employed in this study; no Bragg edges could be characterised for the low-attenuation materials in the TNV reconstruction.
Our study is closely related to, and builds on, the proof-of-concept study presented in [17] using the IMAT beamline, a similar physical phantom (containing Ni, Fe, Cu, Ti, Pb, Al) and channel-wise  211)). In the current study however the dataset was acquired with three times shorter exposure time, two times higher spectral resolution and 25 times smaller scattering volume (0.275 × 0.275 × 0.055 mm 3 regionof-interest vs. 0.055 3 mm 3 voxel) for Bragg edge characterisation. This represents counting statistics some 150 times worse than in the proof-of-concept study yet comparable image quality and edge detection is observed. This is what is achievable by replacing the conventional FBP reconstruction methods with TNV or our proposed TV-TGV method. For the Bragg edges characterised in [17], the absolute error in Bragg edge location in the present study is within 0.03 Å; for other less prominent Bragg edges, which were not characterised in [17], error is within 0.1 Å. There are two well-known limitations of iterative reconstruction techniques. First, they present a significant computational burden. Depending on implementation and available hardware, reconstruction of a four-dimensional dataset can take a few days. However for ToF neutron CT, where typical exposure time per projection is between 30 minutes and an hour, and where the beamtime cost for every measurement hour is extremely high, slow reconstruction is an acceptable price to pay for increased sample throughput. Secondly, choice of a regularisation parameter, which balances the fitting and regulatisation terms in the reconstruction procedure, is a topic of extensive research in the inverse problems community. Here, we manually tuned regularisation parameters based on visual inspection. Automated procedures to define the regularisations parameters are needed to run the proposed methods on an everyday basis.
The method proposed in this paper is by no means limited to ToF neutron CT but is applicable to other spectral CT modalities. We expect that further quality improvements can be achieved if the full four-dimensional spatio-spectral volume is reconstructed as regularisation along the axial direction will be used to penalise image variations and suppress noise. Furthermore, we plan to join reconstruction and material decomposition, such that the spatial volume fraction maps of basis materials are reconstructed directly, instead of reconstruction of hundreds of channels followed by the material decomposition in the image domain. This approach allows a significant reduction of required computations and storage and also have the potential to further improve image quality through the additional regularisation.
There is an ongoing work on providing the reconstruction methods described in this paper to IMAT scientists and users through MANTID Imaging [61] -a graphical toolkit for processing neutron imaging and tomography data. Relevant CIL modules [26,27] will be integrated as MANTID Imaging plugins, supporting not only high-resolution ToF neutron CT investigations and a more efficient usage of valuable neutron beamtime but also bridging a gap between theory and applications of advanced reconstruction methods.

Funding
This work was funded by EPSRC grants "A Reconstruction Toolkit for Multichannel CT" (EP/P02226X/1), "CCPi: Collaborative Computational Project in Tomographic Imaging" (EP/M022498/1 and EP/T026677/1). We gratefully acknowledge beamtime RB1820541 at the IMAT Beamline of the ISIS Neutron and Muon Source, Harwell, UK. EA was partially funded by the Federal Ministry of Education and Research (BMBF) and the Baden-Württemberg Ministry of Science as part of the Excellence Strategy of the German Federal and State Governments. JSJ was partially supported by The Villum Foundation (grant no. 25893). WRBL acknowledges support from a Royal Society Wolfson Research Merit Award. PJW acknowledges support from the European Research Council grant No. 695638 CORREL-CT.