Crystalline phase discriminating neutron tomography using advanced reconstruction methods
Evelina Ametova, Genoveva Burca, Suren Chilingaryan, Gemma Fardell, Jakob S. Jørgensen, Evangelos Papoutsellis, Edoardo Pasca, Ryan Warr, Martin Turner, William R. B. Lionheart, Philip J. Withers
CCrystalline phase discriminating neutron tomography usingadvanced reconstruction methods
Evelina Ametova , Genoveva Burca , Suren Chilingaryan , Gemma Fardell ,Jakob S. Jørgensen , Evangelos Papoutsellis , Edoardo Pasca , Ryan Warr ,Martin Turner , William R. B. Lionheart , and Philip J. Withers Henry Royce Institute, Department of Materials, The University of Manchester, M13 9PL, United Kingdom Laboratory for Application of Synchrotron Radiation, Karlsruhe Institute of Technology, Germany ISIS Pulsed Neutron and Muon Source, STFC, UKRI, Rutherford Appleton Laboratory, United Kingdom Department of Mathematics, The University of Manchester, M13 9PL, United Kingdom Institute for Data Processing and Electronics, Karlsruhe Institute of Technology, Germany Scientific Computing Department, STFC, UKRI, Rutherford Appleton Laboratory, United Kingdom Department of Applied Mathematics and Computer Science, Technical University of Denmark, Denmark Research IT Services, The University of Manchester, M13 9PL, United Kingdom * Correspondence should be addressed to: [email protected]
February 16, 2021
Abstract
Time-of-flight neutron imaging offers complementary attenuation contrast to X-ray computedtomography (CT), coupled with the ability to extract additional information from the variationin 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 theidentification of crystalline phases directly. Here we demonstrate Bragg edge tomography with highspatial and spectral resolution. We propose a new iterative tomographic reconstruction methodwith 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 separatedmode 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 againstFBP 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 bothconventional and state-of-the-art methods and therefore facilitates Bragg edge fitting at the voxellevel. 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 thefull utilisation of state-of-the-art high resolution detectors.
Neutron radiography [1] and later neutron computed tomography (CT) [2, 3] have been availableat neutron sources for some time. As uncharged particles, neutrons probe the nucleus rather thanthe electron cloud and so in contrast to X-ray CT [4] where the attenuation contrast increases withatomic number, neutron attenuation can give strong contrast between neighbouring elements in theperiodic table ( e.g.
Cu and Ni) or even different isotopes [5]. In addition to conventional attenuationcontrast, the transmitted beam contains important information regarding coherent scattering fromthe crystalline lattice. Indeed, for many materials the scattering cross-section is comparable to theattenuation. Since the thermal neutrons produced by neutron spallation sources have wavelengths1 a r X i v : . [ phy s i c s . m e d - ph ] F e b omparable to interplanar lattice spacings, polycrystalline materials will scatter the incident beamelastically according to Bragg’s law. Since Bragg scattering can occur only for wavelengths shorterthan twice the spacing between the lattice planes ( d hkl ), the transmitted neutron spectrum exhibitscharacteristic abrupt increases in the transmitted intensity at these wavelengths. As a result thetransmitted spectrum has a characteristic signature displaying distinct jumps in transmitted intensitycorresponding 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 ofeach detected neutron from its time of flight (ToF) from the source to the detector. The developmentof pixelated ToF detectors enabled the Bragg edges to be imaged [6, 7]. The relation between spectralfingerprint and crystalline properties allows identification of polycrystalline materials and characterisa-tion of their properties [7, 8, 9] such as phase [10, 11], texture and strain [6, 12, 13, 14]. Combined withsample rotation, three-dimensional Bragg edge neutron CT is a natural extension of two-dimensionalBragg edge neutron imaging, as has been already demonstrated in several studies [15, 16, 17]. How-ever, as Kockleman et al. [8] pointed out in their early preliminary neutron energy selective imagingexperiments, the signal, which is usually integrated over the white beam, becomes particularly lowwhen distributed across many hundreds of ToF (energy) channels making tomographic imaging slowand noisy. In recent years the size of pixelated ToF discriminating detectors has improved significantly(from a × array of × [6] to a × array of . × .
055 mm [18]), therebyimproving the potential resolution, but further decreasing the flux received by each pixel and furtherconfounding 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 aresult of the low count rates and hence slow acquisition, energy dispersive spectra tend to be heavilydown-sampled prior to reconstruction and, even with reduced resolution, averaging over a relativelylarge region-of-interest in reconstructed images is still required to improve signal-to-noise ratio andcharacterise Bragg edges. Artificially induced low resolution thus currently hinders application of thetechnique.Here we present advanced reconstruction methods that make energy-dispersive neutron tomogra-phy a practical proposition. We explore the application of dedicated multi-channel reconstructiontechniques with inter-channel correlation to improve reconstruction quality without compromisingvaluable resolution. Contrary to previous work [16, 17] where each channel is treated as a separatereconstruction problem, we use iterative methods with prior information to jointly reconstruct spectralimages.Iterative methods formulate reconstruction as an optimisation problem consisting of a data fidelityterm and a regulariser. The latter term encodes prior knowledge about the image and helps to finda unique solution for the ill-posed tomographic problem by encouraging desired properties in recon-structed images. There is no unique regulariser which performs well for all tomographic problemsand it has to be carefully selected depending on the underlying tomographic data properties. In thisstudy we tailored a regulariser specifically for Bragg edge neutron CT based on a combination ofTotal 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 byencouraging piece-wise constant regions in reconstructed spatial images [19, 20], while TGV regulari-sation 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 regulariserfor multi-channel CT images – Total Nuclear Variation (TNV). TNV is a recent regulariser whichenforces common edges across all channels in multi-channel images [23, 24, 25]. Wider implementationof these advanced reconstruction methods is possible through the open source CCPi Core ImagingLibrary (CIL) reconstruction framework [26, 27], for example one may employ TV only or TGV onlyif desired. The performance of all methods is demonstrated on a multi-material sample comprisingaluminium cylinders filled with various metallic powders of high purity.Experimental details with a particular focus on preprocessing are presented in section 2. In sec-2ion 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 fittingprocedure is further used to detect and locate Bragg edges in TNV and TV-TGV reconstructions andextract crystallographic information. We also demonstrate decomposition of reconstructed spectralimages into individual material maps. Discussion and conclusions are given in section 5. The measurement model in ToF neutron CT is given by the Beer-Lambert law, which relates attenu-ation properties of material to measured intensity: I = I exp (cid:18) − (cid:90) L µ ( x, λ )d x (cid:19) , (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 attenua-tion coefficient at the physical position x in the object for the given wavelength λ . The probability ofneutron-matter interaction is a function of neutron energy (wavelength) and is given by the microscopictotal cross-section σ tot ( λ ) , [cm ] of the nucleus. When neutrons travel through material, the proba-bility of interaction depends not only on the microscopic total cross-section σ tot ( λ ) , but also on thenumber N, [atoms / cm ] of nuclei within unit volume of material. The macroscopic total cross-section Σ tot ( λ ) , [cm − ] defines the probability of neutron-matter interaction per unit distance travelled of theneutron [28], i.e. the attenuation coefficient: µ ( λ ) = Σ tot ( λ ) = N σ tot ( λ ) (2)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 coherentelastic scattering is responsible for characteristic abrupt increases in the transmitted intensity at d hkl allowing the detection of specific lattice planes. In fig. 1 we show the wavelength-dependent macroscopiccross-section Σ tot ( λ ) for materials employed in this study. For this study, a sample was constructed comprising six 6.3 mm diameter thin-walled cylindricalcontainers 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 containerswere 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 affectthe shape of the Bragg edges [31]. This provides idealised Bragg edge spectra likely to be in goodagreement with theoretical calculations of neutron transmission (fig. 1).
The data was acquired at the Imaging and Materials Science & Engineering (IMAT) beamline operatingat the ISIS spallation neutron source (Rutherford Appleton Laboratory, U.K.) [32, 33]. IMAT operatesin a ToF measurement mode. Neutrons generated through spallation are slowed down by the L-H2moderator on IMAT, so they become “thermalised”. Neutrons reach the sample and then the detectorat different times according to their energy. Neutrons with shorter wavelength have higher velocity andare recorded first followed by increasingly less energetic neutrons. Consequently, given the distance3 .5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 λ, [Å] Σ t o t ( λ ) , [ c m − ] FeNiCuZnAl
Figure 1:
Theoretical neutron spectra for materials employed in this study (Calculated using the NXS softwarepackage [30]). travelled between the source and the detector, and the elapsed time from the pulse leaving the sourcethe 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 andthe flight path between the source (more specifically, neutron moderator) and the detector is ≈ . mgiving 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 × array of × readout chips, resulting in × active pixels. The pixel size is 0.055 mm givinga field of view of approximately ×
28 mm . The ToF spectrum recorded at each pixel can haveup to 3100 time (energy) bins. There are small gaps and a slight misalignment between the readoutchips [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 ateach pixel is measured with respect to some external trigger (pulse). Unfortunately, the detector canregister only one event per pixel per time frame. Consequently, slower (lower energy) neutrons have alower probability of being detected, as some pixels are already occupied by faster neutrons. This effectis typically referred to as detector dead time. It is possible to set-up several time frames (also referredto as shutter intervals) within one pulse with arbitrary length and bin width for each time frame, withdata read out in the end of each shutter interval. Data readout takes 320 µ s and introduces gaps inthe measured ToF signal. The introduction of several shutter intervals within one pulse reduces thespectral signal distortions due to detector dead time, but cannot fully eliminate it. In [34] authorsproposed an algorithm called “overlap correction” to compensate for the counts lost. The efficiency ofthe algorithm was confirmed in [35]. For the present study, the detector shutter intervals were selectedin 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 havingwavelength resolutions between . · − Å and . · − Å.4 igure 2:
Sketch showing the test sample employed in the study comprising Al, Fe, Cu, Ni and Zn powders.
Beginning End Bin width Number of channelsms Å ms Å µ s Å15.0000 1.0524 26.6800 1.8719 10.2400 . · − . · − . · − . · − Table 1:
MCP detector shutter intervals chosen for the experimental study.
The sample (fig. 2) was placed and clamped onto the rotary stage. The sample was moved as closeas possible to the MCP detector using the kinematic system. An alignment laser was used to visuallyalign the vertical axis of the sample with respect to the vertical edge of the detector and to ensurethat the object was within the field of view for all rotation angles. A set of spectral projections wereacquired at 120 equally-spaced angular positions over ◦ rotation ( . ◦ angular increments). Eachprojection was acquired with 15 min exposure time.We acquired two sets of flat field (also referred to as open beam) images (before and after ac-quisition) to compensate for detector imperfections and decrease of beam intensity over time whichwas observed in other experiments. Also, flat field averaging is generally recommended to reduce ringartifacts. Therefore 4 spectral flat field images were acquired before and 4 after sample acquisition (8in total), each one with 15 min exposure. The MCP detector has no dark current noise. The raw datais available at [36]. Following [34, 37], MCP detector related corrections (including both overlap correction and scaling tothe same number of incident neutrons) were performed both for the set of 120 projections and the setof 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 of4 flat field images acquired before the data acquisition and the second one corresponds to anaverage of 4 flat field images acquired afterwards. Then, the actual flat field image for every5
Single flat field A n g l e , [ ∘ ] Inter olated flat field
Pixel
Flux normalisationP1 P2 P3 P4 T r a n s m i ss i o n , I / I (a) Rotation po ition, [ ∘ ] T r a n m i i o n , I / I Single flat fieldInterpolated flat fieldFlux normali ation (b)
Figure 3:
Effect of flat field correction. a) Flat-field corrected white beam (sum of all wavelength) transmissionsinogram for a single slice of the acquired dataset. The white vertical line indicates the location of the profileline. b) Intensity profile along vertical profile lines. P1-P4 label pixels chosen for examination of acquiredspectra in the next figure. projection is calculated based on pixel-wise channel-wise linear interpolation of intensity valuesbetween 2 averaged flat field images.• Using an average of all 8 flat field images. The intensity of the averaged flat field image is thenmultiplied by the quotient of the mean intensity in unoccupied detector columns in a projectionimage and mean intensity in the averaged flat field image itself. The scaling factor is calculatedindividually 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-fieldcorrected white beam (sum of all wavelength) transmission sinogram for a single slice of the acquireddataset (fig. 3). Noticeable horizontal stripes are present for the first two flat field correction approacheswhich are eliminated with the flux normalisation. Profile lines (marked as a white solid line in thesurrounding air in fig. 3a) spotlight the difference between the first two approaches: linear interpolationcompensates for the upwards trend in image intensity (fig. 3b). Finally, the flux normalisation notonly de-trends the signal but also drastically reduces the overall fluctuations. As a result, the fluxnormalisation approach was used in the present study.As relatively short exposure time was used for the acquisition, spectral images were down-sampled inthe spectral dimension to increase counting statistics and produce reasonable spectra. Down-samplingwas performed for each shutter interval independently by taking the average of each 16, 8, 8 and4 channels for the respective shutter intervals (table 1). As a result, each preprocessed projectioncontains 339 channels having a uniform bin width of . · − Å over a wavelength range between1 Å 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. Recordedspectra (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 Braggedge characterisation for individual pixels barely feasible.It is also important to recognise that the incident spectrum on IMAT is a function of the neutronmoderation 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 Å, where6
Pixel Σ t o t ( λ ) , [ c m − ] R o t a t i o n p o s i t i o n , [ ∘ ] (a) −0.250.000.25 P1 - Air P2 - Fe P3 - Fe+Ni
Fe Ni1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 λ, [Å] P4 - Cu Σ t o t ( λ ) , [ c m − ] (b) Figure 4:
Sinograms (a) and recorded spectra (b). The spectra have been normalised by the transmission pathlength through the corresponding material. P1-P4 in (a) denote the location of pixels chosen for examination ofrecorded spectra. Vertical dotted lines in the bottom subplot in (b) mark the channels chosen for visualisationin (a).
IMAT has the lowest flux, and between 2.5 Å and 3 Å where a combination of high flux and detectordead time causes distortions in the recorded spectrum. The wavelength-noise dependence has directimplications on the reconstruction quality of each individual channel.Short of further downsampling of the recorded multi-channel images, advanced reconstruction meth-ods are needed to handle noisy multi-channel neutron CT data. A potential solution to this low-countCT problem is to effectively exploit available prior information, such as the anticipated image proper-ties and structural correlation between channels.
Tomographic reconstruction aims to recover a map of sample attributes from a set of integral mea-surements acquired at various angles. In every wavelength channel, Bragg edge neutron CT can bewell approximated by the standard absorption tomography model (Radon transform). Consequently,algorithms developed for X-ray CT can be conveniently used to reconstruct the neutron attenuationmap for every wavelength channel. Here, we focus on the two-dimensional reconstruction problem inthe spatial domain; extension to the third spatial dimension is straightforward.Given a multi-channel tomographic dataset with K channels, a conventional approach is to recon-struct each channel independently using the FBP algorithm. The FBP algorithm is derived from theFourier Slice theorem which relates line integral measurements to two-dimensional Fourier transform ofan object’s slice. In FBP-type reconstruction methods, projections are filtered independently and thenback-projected onto the plane of the tomographic slice. FBP reconstruction is very fast but requireshigh 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. Con-sider a single channel k, k = 0 , . . . , K − , 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 number7f projections and D being the number of pixels in a detector row. The sinogram b ( λ k ) is vectorisedas a column vector with M = P D elements (obtained by stacking the columns in the original two-dimensional sinogram). Let u ( λ k ) be a vectorised two-dimensional array of material attributes wewant to reconstruct with N = D elements (voxels). The discrete version of the Radon transformis then given by the projection operator A containing M × N elements. If i, i = 0 , , . . . M − and j, j = 0 , , . . . , N − , then A ( i,j ) is the length of intersection of the i .th ray with the j .th voxel. Thereconstruction problem then takes a form: b ( λ k ) = − ln (cid:18) I ( λ k ) I ( λ k ) (cid:19) ≈ Au ( λ k ) . (3)where I ( λ 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 typicallylow (compared to synchrotron X-ray sources) and because the detected neutrons are shared amongmultiple time of flight (energy) channels. Hence, the acquisition mode does not inherently providesufficient data for a numerically stable solution of eq. (3), similar to the low-dose tomography problemin 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 ill-posed problems is to design a surrogate problem which is consistent with the original problem but iswell-posed and computationally tractable. In other words, we obtain a stable approximate solution byincorporating some prior knowledge about the problem.Multi-channel CT data can be considered as a stack of two- or three-dimensional datasets, whereevery voxel in the reconstructed volume contains a vector with a spectral material response, spanningan additional dimension. In this respect, every data point in the material response vector belongs in thesame physical structure, i.e. channels share structural information. Then, inter-channel correlationscan be utilised to improve reconstruction quality. The spectral CT reconstruction takes a form: b = Au, (4)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: arg min u { F ( u ) = f ( b, Au ) + αg ( u ) } , (5)where f ( b, Au ) is a data fidelity metric which measures the discrepancy between the projection ofsolution u and the acquired data b . The regularisation term g ( u ) imposes certain prior assumptionstypically expressed in terms of desired image properties. The scalar parameter α controls the trade-offbetween fit with the acquired data u and the regularisation. The choice of regulariser depends on theunderlying problem since there is no unique regulariser which performs best in all problems. The artof choosing the right regulariser and finding a balance between the two terms is a challenging task, asregularisation inevitably introduces bias into the solution. The bias is a price one pays for solving aproblem which is not solvable by other means [38].TV is the most commonly used regulariser. TV encourages a sparse image gradient and conse-quently favours piece-wise constant images with sharp boundaries [19, 20]. The TV model describeswell the spatial dimension of images of the type to be reconstructed in this study, but it is known tointroduce “staircase” artifacts for piece-wise affine signals, i.e. a ramp or features similar to the onesobserved in wavelength dependent attenuation coefficient (fig. 1) might be reconstructed as piece-wiseconstant, reminiscent of a staircase. TV is also known to suffer from intensity reductions and, as aresult, low contrast regions might be lost. Multiple modifications of TV have been proposed to curethe above shortcomings. 8he TNV approach proposed in [24] explicitly promotes reconstructions with common edges acrossall channels. The idea is based on the fact that having shared gradient directions is equivalent to havea rank-one Jacobian of a multi-channel image. Therefore, TNV penalises the singular values of theJacobian. In the spatial dimension, TNV has similar properties to TV regularisation, i.e. it also favoursa sparse image gradient. Consequently, TNV correlates channels and improves reconstruction qualityby promoting common structures in multichannel images. Application of TNV for reconstruction ofmulti-channel images has been successfully demonstrated in [24, 25, 39]. The reconstruction problemis then formulated as, F ( u ) = (cid:107) Au − b (cid:107) + α TNV( u ) . (6)Similar to TV, TNV suffers from a loss of contrast. Secondly, TNV does not allow the decoupling ofregularisation parameters for the spatial and spectral dimensions which makes it impossible to balancethe level of regularisation between dimensions.Here we propose a novel tailored regulariser which treats the spatial and spectral dimensions sep-arately. As the TV model captures piece-wise image properties in the spatial dimension, we rely onanother regulariser to support reconstruction in the spectral dimension. A natural remedy for the stair-case artifact is to use higher order derivatives in the regularisation term. Here, we use TGV [21] whichallows for both sharp changes in spectral signal and gradual intensity changes. TGV also effectivelyexploits inter-channel correlation. In this case the reconstruction problem is formulated as, F ( u ) = (cid:107) Au − b (cid:107) + β TV x,y ( u ) + γ TGV c ( u ) . (7)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 channelsindividually and then summed over all channels. TGV c is a one-dimensional operator applied in eachindividual 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-rayspectroscopy (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 caseBragg edge neutron CT data. The decoupling of regularisation into the spatial and spectral dimensionallows us to balance the two regularisation terms appropriately and promote the desired properties inboth 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 forparallelisation and software acceleration. A common way to solve the large-scale optimisation problem in eq. (5) is by means of the FastIterative 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 ImagingLibrary (CIL) [26, 27]. CIL provides a highly modular Python library for prototyping of reconstructionmethods for multi-channel data such as spectral and dynamic (time-resolved) CT. CIL wraps a numberof 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 andprovides a set of various regularisers through the CCPi Regularisation Toolkit [50]. TNV was solvedusing FISTA employing the CIL plugin for the CCPi Regularisation Toolkit [50] for the TNV proximaloperator. TV-TGV was implemented directly in CIL and solved using PDHG.Regularisation parameters were carefully chosen to achieve both noise suppression and featurepreservation in both spatial and spectral dimensions. For TNV reconstruction, the regularisationparameter α was set to 0.01, while TV-TGV parameters were set to β = 0 . and γ = 0 . . In bothcases, we ran 2000 iterations of the reconstruction algorithm.9 .3 Post-processing and analysis Quantitative analysis of the reconstructed spectra commonly includes detection and characterisationof the Bragg edges. The Bragg edge positions (in terms of d -spacing) allows compositional mapping, aseach crystalline structure will have unique set of lattice spacings and hence fingerprint in the neutrontransmitted spectrum. The shape of the detected Bragg edges, i.e. deviation from abrupt step-liketheoretical edge, supports characterisation of crystallographic properties.In [37] 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 fromBEAn 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 recon-structed 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 edgelocations. Note, smoothing is performed only to support the peak-finding routine and all further stepsare performed on unsmoothed spectra. (3) For each detected Bragg edge, a model function derivedin [51] is fitted using the non-linear least squares fitting method. Since the model function is highlynon-linear, the fitting procedure is very sensitive to initialisation parameters. Therefore we use bruteforce approach to help the fitting procedure to achieve acceptable fit. That is, initialisation parametersare drawn from uniform random distributions chosen to cover a realistic range for each parameter inthe model function. The best fit is chosen based on minimum root mean squared error between thefitted function and experimental data.Additionally, the spatial distribution of individual materials in a sample can be obtained by de-composing the reconstructed spectral images into individual material maps. Material decompositionin spectral CT is an ill-posed problem by its own and various approaches have been proposed in theliterature, mainly for spectral X-ray CT. The first group of methods performs material decompositionin the sinogram domain, following by reconstruction of material specific sinograms [52, 53]. The sec-ond group of methods employs simultaneous reconstruction and material decomposition [54]. A thirdapproach is to perform material decomposition in the image domain, i.e. after the reconstructionstep [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 volumeconservation principle [57, 58], where each voxel in the obtained material maps corresponds to a di-mensionless volume fraction occupied by the corresponding material with the full voxel correspondingto a unit volume. This is a valid assumption here as the materials employed in the current study donot mix. Under the volume conservation assumption, the voxel-wise sum of all material maps has tobe equal to 1 in each voxel.Let ˆ u be a K × D matrix version of u and let M be the material basis K × S matrix built up fromthe predicted wavelength-dependent neutron spectra (fig. 1), with S = 6 being the number of materialsemployed in this study (5 metals + air). Then, volume fractions of every material v , consisting of S row vectors with D elements, can be obtained by solving: ˆ u ≈ M v. (8)To solve eq. (8), we formulate the following optimisation problem: min v (cid:107) ˆ u − M v (cid:107) F s.t. v ≥ TS v = TD (9)where (cid:107) · (cid:107) F is the Frobenius norm and T is a row vector of all ones of appropriate dimension. Thefirst constraint is an element-wise non-negativity constraint enforcing the volume fractions v to begreater than or equal to 0; the second constraint expresses the volume preservation principle. We useCVX [59, 60] with the MOSEK solver to solve this optimisation problem.10 Results
In order to give a visual reference for the spatial dimension, we first perform conventional FBP recon-struction (implemented as an eponymous processor in CIL [26]) of the white beam sinogram (sum of allwavelength channels), as shown in fig. 5. All cylinders are seen clearly defined and with different con-trast for each element. The copper powder used in this study is known to have an average particle sizecomparable 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 occupiedby 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 makingthe phantom. There are also small “bumps” where a small portion of powder penetrated overlappingfoil layers.
Voxel V o x e l CuAl Fe NiZn μ , [ c m − ] Figure 5:
White beam (sum of all wavelength channels) reconstruction using the conventional FBP method.
Figure 6 shows two-dimensional slices for selected (individual) wavelength channels reconstructed usingthe 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 recon-struction methods. We only perform two-dimensional reconstructions in the spatial dimension in thisstudy but results can be generalised to the third dimension. We have chosen a slice roughly in themiddle of the upper half of the detector plate in order to avoid detector rows near to the top of thedetector or the gap between the detector read-out chips (section 2.3).As expected, at these signal levels the FBP reconstruction is barely interpretable (fig. 6, toprow). Both TNV and TV-TGV demonstrate drastic improvement in reconstruction quality and noisesuppression (figs. 6 and 7). A region-of-interest between two cylinders highlights significant smearingof features in the TNV reconstruction (fig. 6, middle row); the same features appear sharper in theTV-TGV reconstruction (fig. 6, bottom row). TNV suffers from the same limitations as conventionalTV-based regularisation: it might oversmooth images. Therefore, TNV suppresses noise and ringartifacts visible in FBP reconstruction but produces blurred and enlarged rings especially prominent11 B P T N V T V - T G V C Al Fe NiZn Σ t o t ( λ ) , [ c m − ] Figure 6:
Two dimensional reconstructions of selected (individual) wavelength channels. Top to bottom:channel-wise FBP reconstruction, iterative reconstruction with TNV regularisation and iterative reconstructionwith TV-TGV regularisation. A magnified region-of-interest marked with a white rectangle is shown inset.White lines mark profile lines chosen for examination in the next figure, whereas white dots mark individualvoxels chosen for spectral comparison in the next section. All slices are visualised using a common colourrange. The colour range in the magnified insets is scaled individually to cover only the voxel values in thecorresponding inset. in shorter and longer wavelength channels, where counts are much lower. Aluminium has very lowneutron attenuation and is invisible in the TNV reconstruction due to contrast loss; another knowndrawback 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 isvisible in the reconstructed slices and profile lines (fig. 7, bottom row). Fine features inside the coppercylinder are also partially preserved in the TV-TGV reconstruction (fig. 7, top row).
The ground truth for the spectral dimension is taken to be the predicted neutron cross-section for theselected materials (fig. 1). Individual spectra reconstructed for one . mm voxel in each of the5 materials are plotted in fig. 8 alongside the theoretical predictions. The voxel locations inside thecylinders were chosen arbitrarily (marked in fig. 6).The limited counting statistics acquired in our experiment make the FBP reconstructed spectrauninterpretable. By contrast, both the TNV and TV-TGV reconstructions show drastic improvementsin reconstruction quality. The reconstructed spectra for Fe, Ni and Cu closely follow the Bragg edge12 FBPTNVTV-TGV
Al Ni Σ t t ( λ ) , [ c m − ] Position along profile, [mm]
Figure 7:
Linescans corresponding to the vertical (top row) and horizontal (bottom row) white lines in fig. 6(bottom second left) passing through the Cu and Fe cylinders and the Al and Ni cylinders, respectively 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 between2.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 per-formance. While the TV-TGV produces a much smoother spectra, the Bragg edges appear to be lessprominent due to smearing (for instance, small edges around 2 Å in Fe). In the case of the TNVregularisation, the noise dominates over smaller Bragg edges. For the low-attenuation materials (Aland Zn), TV-TGV shows better performance as some Bragg edges are visible in the reconstructedspectrum, but are lost in the noise in the case of TNV.
Mapping of crystallographic information was performed on individual spectra reconstructed for one . mm 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 basedon fitting of a dedicated model function (fig. 9). Fitting was performed individually around eachdetected Bragg edge. The fitted function is shown as the orange solid line. Since the locations of theBragg edges are directly related to interplanar spacing d hkl , we use the reference λ hkl = 2 d hkl for thecorresponding materials to assess the mapping of the crystalline structural information (table 2). Wealso show an absolute error ε = (cid:107) λ hkl − λ (cid:48) hkl (cid:107) in the estimated d -spacing λ (cid:48) hkl . For some Bragg edges,the edge detection routine identified edges, but the least squares fitting procedure did not converge toa 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 recon-structions 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 theTNV results. Thus, in the TNV reconstruction we can estimate only some interplanar distances; noBragg edges can be detected and characterised for Al. Remarkably, the TV-TGV reconstruction isable to locate some Bragg edges even for the low-attenuation materials such as Zn and Al. Accuracyin 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-13 .51.01.5 F e FBP λ, Å
TNV TV-TGV N i λ, Å C u λ, Å Z n λ, Å A l λ, Å Σ t o t ( λ ) , [ c m − ] Figure 8:
Individual spectra (solid green line) reconstructed by FBP, TNV and TV-TGV for one representative . mm voxel located within in each material alongside the predicted signal (dotted black line). 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 themodel function (section 3.3). A more careful manual procedure might improve the Bragg edge fittingresults, especially for the noisier TNV reconstruction. Finally, in fig. 10 we show spatial material maps calculated on the basis of the material decompositionmethod described in section 3.3. The difference between the material maps calculated from TNV andTV-TGV reconstructed images is quite subtle, especially for high-attenuation materials. Advantagesof 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 moreuniform. The material maps for both reconstructions exhibit some artifacts, for instance, the faintFe cylinder is visible in the Ni map and similarly Zn is visible in the Cu map. The artifacts arecaused by the fact that the most prominent Bragg edges almost coincides in the Fe and Ni spectra, aswell as in the Cu and Zn spectra (fig. 1). In this study we used a fairly generic method for materialdecomposition. As we have already mentioned in section 3.3, the material decomposition problem is anill-posed problem and various regularisation methods can be used to improve material decomposition14 .51.01.5 F e ( )( )( ) TNV ( )( )( )( ) TV-TGV N i ( )( )( )( ) ( )( )( )( )( ) C u ( )( ) ( )( )( )( )( ) Z n ( )( )( )( ) A l ( )( )( )( ) Σ t o t ( λ ) , [ c m − ] λ, [Å] Figure 9:
Bragg edge fitting results. Reconstructed spectra (solid green line) are superimposed upon thepredicted signal (dotted black line). The model function was fitted around each Bragg edge separately. Thefitted model is shown as the solid orange line. Crystalline structural information is included in the reconstructedspectra. depending on underlying image properties in the spatial and spectral dimensions.
The unique features of the pulsed time of flight neutron source allow quantification of material prop-erties in bulk samples. Low flux and slow acquisition result in long scan times being required torecord sufficient counts across hundreds of wavelength channels, making material characterisation atthe full attainable resolution infeasible in practice. Joint reconstruction which exploits inter-channelcorrelations in multi-channel images is an effective way to shift feature recovery from exposure to re-construction. In this work we have studied the performance of advanced reconstruction methods withthe dedicated multi-channel regularisation techniques for Bragg edge neutron CT. We have demon-strated that the tailored TV-TGV regularisation technique, which favours specific image properties inthe spatial and spectral dimensions, allows retrieval of crystallographic properties at a resolution pre-viously unattainable through conventional reconstruction methods using the same exposure time. Theproposed technique was compared with the TNV method - a recent regularisation technique developedfor spectral CT. While both methods clearly outperform FBP, TV-TGV yields better reconstruction15 hkl
TNV TV-TGV λ (cid:48) hkl ε λ (cid:48) hkl ε Fe(110) 4.054 4.055 0.001 4.077 0.023Fe(200) 2.866 2.819 0.047 2.853 0.013Fe(211) 2.340 2.366 0.026 2.332 0.008Fe(220) 2.026 - (-)Fe(310) 1.812 (-) -Fe(222) 1.654 - -Fe(321) 1.532 - 1.492 0.04Ni(111) 4.064 4.073 0.009 4.076 0.012Ni(200) 3.520 3.533 0.013 3.545 0.025Ni(220) 2.490 2.522 0.032 2.512 0.022Ni(311) 2.122 2.122 0.000 2.136 0.014Ni(222) 2.032 - -Ni(400) 1.760 - -Ni(331) 1.616 (-) 1.609 0.007Cu(111) 4.168 4.216 0.048 4.151 0.017Cu(200) 3.610 (-) 3.662 0.052Cu(220) 2.552 (-) 2.627 0.075Cu(311) 2.176 2.163 0.013 2.215 0.039Cu(222) 2.084 - -Cu(400) 1.804 - -Cu(331) 1.656 - 1.629 0.027Al(111) 4.676 - 4.729 0.053Al(200) 4.050 - 4.040 0.01Al(220) 2.864 - (-)Al(311) 2.442 - 2.435 0.007Al(222) 2.338 - -Al(400) 2.024 - -Al(331) 1.858 - 1.982 0.124Zn(002) 4.950 - -Zn(100) 4.608 (-) 4.648 0.040Zn(101) 4.178 (-) 4.209 0.031Zn(102) 3.372 (-) 3.391 0.019Zn(103) 2.682 (-) 2.671 0.011
Table 2:
Results of mapping of crystallographic information. Estimated λ (cid:48) hkl is compared vs. reference λ (cid:48) hkl = 2 d hkl , where d hkl were extracted from NXS code [30]. For some Bragg edges, the edge detection routineidentified edges but the least squares fitting procedure did not converge to a good solution. We use “(-)” todesignate 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 per-formance of both methods. We extracted crystallographic properties from the reconstructed spectrabased on detection and characterisation of Bragg edges. TV-TGV facilitated extraction of interpla-nar spacings for all materials employed in this study; no Bragg edges could be characterised for thelow-attenuation materials in the TNV reconstruction.Our study is closely related to, and builds on, the proof-of-concept study presented in [17] usingthe IMAT beamline, a similar physical phantom (containing Ni, Fe, Cu, Ti, Pb, Al) and channel-wiseFBP reconstruction. The authors demonstrated Bragg edge fitting with an absolute error in Bragg16 N V Fe Ni Cu Al Zn Air T V - T G V Figure 10:
Material maps. Intensity value in each voxel is equal to the voxel volume fraction occupied bythe corresponding material (calculated based on the material decomposition method described in section 3.3). edge location within 0.04 Å for Ni ((111), (200), (220), (311)) and Fe ((110), (200), (211)). In thecurrent study however the dataset was acquired with three times shorter exposure time, two timeshigher spectral resolution and 25 times smaller scattering volume ( . × . × .
055 mm region-of-interest vs. . mm voxel) for Bragg edge characterisation. This represents counting statisticssome 150 times worse than in the proof-of-concept study yet comparable image quality and edgedetection is observed. This is what is achievable by replacing the conventional FBP reconstructionmethods with TNV or our proposed TV-TGV method. For the Bragg edges characterised in [17], theabsolute error in Bragg edge location in the present study is within 0.03 Å; for other less prominentBragg edges, which were not characterised in [17], error is within 0.1 Å.There are two well-known limitations of iterative reconstruction techniques. First, they presenta significant computational burden. Depending on implementation and available hardware, recon-struction of a four-dimensional dataset can take a few days. However for ToF neutron CT, wheretypical exposure time per projection is between 30 minutes and an hour, and where the beamtimecost for every measurement hour is extremely high, slow reconstruction is an acceptable price to payfor increased sample throughput. Secondly, choice of a regularisation parameter, which balances thefitting and regulatisation terms in the reconstruction procedure, is a topic of extensive research inthe inverse problems community. Here, we manually tuned regularisation parameters based on visualinspection. Automated procedures to define the regularisations parameters are needed to run theproposed methods on an everyday basis.The method proposed in this paper is by no means limited to ToF neutron CT but is applicableto other spectral CT modalities. We expect that further quality improvements can be achieved if thefull four-dimensional spatio-spectral volume is reconstructed as regularisation along the axial directionwill be used to penalise image variations and suppress noise. Furthermore, we plan to join recon-struction and material decomposition, such that the spatial volume fraction maps of basis materialsare reconstructed directly, instead of reconstruction of hundreds of channels followed by the materialdecomposition in the image domain. This approach allows a significant reduction of required computa-tions and storage and also have the potential to further improve image quality through the additionalregularisation.There is an ongoing work on providing the reconstruction methods described in this paper toIMAT scientists and users through MANTID Imaging [61] – a graphical toolkit for processing neutronimaging and tomography data. Relevant CIL modules [26, 27] will be integrated as MANTID Imagingplugins, supporting not only high-resolution ToF neutron CT investigations and a more efficient usageof valuable neutron beamtime but also bridging a gap between theory and applications of advancedreconstruction methods. 17 unding 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 IMATBeamline of the ISIS Neutron and Muon Source, Harwell, UK. EA was partially funded by the Fed-eral Ministry of Education and Research (BMBF) and the Baden-Württemberg Ministry of Scienceas part of the Excellence Strategy of the German Federal and State Governments. JSJ was partiallysupported by The Villum Foundation (grant no. 25893). WRBL acknowledges support from a RoyalSociety Wolfson Research Merit Award. PJW acknowledges support from the European ResearchCouncil grant No. 695638 CORREL-CT.
Acknowledgments
Authors would like to thank Dr. Winfried Kockelmann (ISIS STFC), Dr. Anton Tremsin (UC Berke-ley), Dr. Joe Kelleher (ISIS STFC), Dr. Søren Schmidt (ESS) and Alexander Liptak (University ofLiverpool) for the fruitful discussions on challenges in neutron CT and various approaches to overcomethem. This work made use of computational support by CoSeC, the Computational Science Centrefor Research Communities, through the Collaborative Computational Project in Tomographic Imaging(CCPi).
References [1] Maier-Leibnitz H, Springer T. The use of neutron optical devices on beam-hole experiments onbeam-hole experiments. Journal of Nuclear Energy Parts A/B Reactor Science and Technology.1963;17(4-5):217–225. DOI: 10.1016/0368-3230(63)90022-3.[2] Oien C, Bailey K, Barton C, Guenther R, Barton J. Resonance energy neutron radiography forcomputerized axial tomography. In: Materials Evaluation. vol. 35; 1977. p. 23–23.[3] Schlapper G, Brugger R, Seydel J, Larsen G. Neutron Tomography Investigations at the Mis-souri University Research Reactor. Transactions of the American Nuclear Society Supplement.1977;26(1):39.[4] Maire E, Withers PJ. Quantitative X-ray tomography. International materials reviews.2014;59(1):1–43. DOI: 10.1179/1743280413Y.0000000023.[5] Strobl M, Manke I, Kardjilov N, Hilger A, Dawson M, Banhart J. Advances in neutron radiographyand tomography. Journal of Physics D: Applied Physics. 2009;42(24):243001. DOI: 10.1088/0022-3727/42/24/243001.[6] Santisteban JR, Edwards L, Fitzpatrick ME, Steuwer A, Withers PJ, Daymond MR, et al. Strainimaging by Bragg edge neutron transmission. Nuclear Instruments and Methods in Physics Re-search Section A: Accelerators, Spectrometers, Detectors and Associated Equipment. 2002;481(1-3):765–768. DOI: 10.1016/S0168-9002(01)01256-6.[7] Santisteban JR, Edwards L, Fizpatrick ME, Steuwer A, Withers PJ. Engineering applica-tions of Bragg-edge neutron transmission. Applied Physics A. 2002;74(1):s1433–s1436. DOI:10.1007/s003390101241.[8] Kockelmann W, Frei G, Lehmann EH, Vontobel P, Santisteban JR. Energy-selective neutrontransmission imaging at a pulsed source. Nuclear Instruments and Methods in Physics ResearchSection A: Accelerators, Spectrometers, Detectors and Associated Equipment. 2007;578(2):421–434. DOI: 10.1016/j.nima.2007.05.207. 189] Kardjilov N, Manke I, Hilger A, Strobl M, Banhart J. Neutron imaging in materials science.Materials Today. 2011;14(6):248–256. DOI: 10.1016/S1369-7021(11)70139-0.[10] Steuwer A, Withers P, Santisteban J, Edwards L. Using pulsed neutron transmission forcrystalline phase imaging and analysis. Journal of applied physics. 2005;97(7):074903. DOI:10.1063/1.1861144.[11] Song G, Lin JY, Bilheux JC, Xie Q, Santodonato LJ, Molaison JJ, et al. Characterization ofCrystallographic Structures Using Bragg-Edge Neutron Imaging at the Spallation Neutron Source.Journal of Imaging. 2017;3(4):65. DOI: 10.3390/jimaging3040065.[12] Woracek R, Penumadu D, Kardjilov N, Hilger A, Strobl M, Wimpory R, et al. Neutron Bragg-edge-imaging for strain mapping under in situ tensile loading. Journal of Applied Physics.2011;109(9):093506. DOI: 10.1063/1.3582138.[13] Wensrich CM, Hendriks JN, Gregg A, Meylan MH, Luzin V, Tremsin AS. Bragg-edge neu-tron transmission strain tomography for in situ loadings. Nuclear Instruments and Methods inPhysics Research Section B: Beam Interactions with Materials and Atoms. 2016;383:52–58. DOI:10.1016/j.nimb.2016.06.012.[14] Reid A, Marshall M, Kabra S, Minniti T, Kockelmann W, Connolley T, et al. Application ofneutron imaging to detect and quantify fatigue cracking. International Journal of MechanicalSciences. 2019;159:182–194. DOI: 10.1016/j.ijmecsci.2019.05.037.[15] Woracek R, Penumadu D, Kardjilov N, Hilger A, Boin M, Banhart J, et al. 3D Mapping of Crystal-lographic Phase Distribution using Energy-Selective Neutron Tomography. Advanced Materials.2014;26(24):4069–4073. DOI: 10.1002/adma.201400192.[16] Watanabe K, Minniti T, Sato H, Tremsin AS, Kockelmann W, Dalgliesh R, et al. Cross-sectionalimaging of quenched region in a steel rod using energy-resolved neutron tomography. NuclearInstruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectorsand Associated Equipment. 2019;944:162532. DOI: 10.1016/j.nima.2019.162532.[17] Carminati C, Strobl M, Minniti T, Boillat P, Hovind J, Morgano M, et al. Bragg-edge attenuationspectra at voxel level from 4D wavelength-resolved neutron tomography. Journal of AppliedCrystallography. 2020 Feb;53(1):188–196. DOI: 10.1107/S1600576720000151.[18] Tremsin AS, Vallerga JV, McPhate JB, Siegmund OHW, Raffanti R. High Resolution PhotonCounting With MCP-Timepix Quad Parallel Readout Operating at > Frame Rates. IEEETransactions on Nuclear Science. 2013 April;60(2):578–585. DOI: 10.1109/TNS.2012.2223714.[19] Rudin LI, Osher S, Fatemi E. Nonlinear total variation based noise removal algorithms. PhysicaD: nonlinear phenomena. 1992;60(1-4):259–268. DOI: 10.1016/0167-2789(92)90242-F.[20] Sidky EY, Kao CM, Pan X. Accurate image reconstruction from few-views and limited-angledata in divergent-beam CT. Journal of X-ray Science and Technology. 2006;14(2):119–139.arXiv:0904.4495.[21] Bredies K, Kunisch K, Pock T. Total generalized variation. SIAM Journal on Imaging Sciences.2010;3(3):492–526. DOI: 10.1137/090769521.[22] Kazantsev D, Ovtchinnikov E, Withers PJ, Lionheart WR, Lee PD. Sparsity seeking to-tal generalized variation for undersampled tomographic reconstruction. In: 2016 IEEE 13thInternational Symposium on Biomedical Imaging (ISBI). IEEE; 2016. p. 731–734. DOI:10.1109/ISBI.2016.7493370. 1923] Holt KM. Total nuclear variation and jacobian extensions of total variation for vector fields. IEEETransactions on Image Processing. 2014;23(9):3975–3989. DOI: 10.1109/TIP.2014.2332397.[24] Rigie DS, La Rivière PJ. Joint reconstruction of multi-channel, spectral CT data via constrainedtotal nuclear variation minimization. Physics in Medicine & Biology. 2015;60(5):1741. DOI:10.1088/0031-9155/60/5/1741.[25] Kazantsev D, Jørgensen JS, Andersen MS, Lionheart WR, Lee PD, Withers PJ. Joint image re-construction method with correlative multi-channel prior for x-ray spectral computed tomography.Inverse Problems. 2018;34(6):064001. DOI: 10.1088/1361-6420/aaba86.[26] Jørgensen JS, Ametova E, Burca G, Fardell G, Papoutsellis E, Pasca E, et al. Core Imaging Library– Part I: a versatile Python framework for tomographic imaging. Philosophical Transactions ofthe Royal Society A. 2021; preprint. arXiv:2102.04560.[27] Papoutsellis E, Ametova E, Delplancke C, Fardell G, Jørgensen JS, Pasca E, et al. Core ImagingLibrary – Part II: Multichannel reconstruction for dynamic and spectral tomography. Philosoph-ical Transactions of the Royal Society A. 2021; preprint. arXiv:2102.06126.[28] Fundamentals D. Nuclear Physics and Reactor Theory. DOE-HDBK-1019/1-93; 1993.[29] Boin M, Hilger A, Kardjilov N, Zhang S, Oliver E, James J, et al. Validation of Bragg edgeexperiments by Monte Carlo simulations for quantitative texture analysis. Journal of AppliedCrystallography. 2011;44(5):1040–1046. DOI: 10.1107/S0021889811025970.[30] Boin M. nxs : a program library for neutron cross section calculations. Journal of AppliedCrystallography. 2012 Jun;45(3):603–607. DOI: 10.1107/S0021889812016056.[31] Xie Q, Song G, Gorti S, Stoica AD, Radhakrishnan B, Bilheux JC, et al. Applying neutrontransmission physics and 3D statistical full-field model to understand 2D Bragg-edge imaging.Journal of Applied Physics. 2018;123(7):074901. DOI: 10.1063/1.5013676.[32] Kockelmann W, Minniti T, Pooley DE, Burca G, Ramadhan R, Akeroyd FA, et al. Time-of-FlightNeutron Imaging on IMAT@ISIS: A New User Facility for Materials Science. Journal of Imaging.2018;4(3). DOI: 10.3390/jimaging4030047.[33] Burca G, Kockelmann W, James JA, Fitzpatrick ME. Modelling of an imaging beamline at theISIS pulsed neutron source. Journal of Instrumentation. 2013 oct;8(10):P10001–P10001. DOI:10.1088/1748-0221/8/10/P10001.[34] Tremsin AS, Vallerga JV, McPhate JB, Siegmund OHW. Optimization of Timepix count ratecapabilities for the applications with a periodic input signal. Journal of Instrumentation. 2014may;9(05):C05026–C05026. DOI: 10.1088/1748-0221/9/05/C05026.[35] Watanabe K, Minniti T, Kockelmann W, Dalgliesh R, Burca G, Tremsin AS. Characterizationof a neutron sensitive MCP/Timepix detector for quantitative image analysis at a pulsed neutronsource. Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrom-eters, Detectors and Associated Equipment. 2017;861:55 – 63. DOI: 10.1016/j.nima.2017.04.034.[36] Jørgensen JS, Ametova E, Burca G, Fardell G, Papoutsellis E, Pasca E, et al. Neutron TOFimaging phantom data to quantify hyperspectral reconstruction algorithms. STFC ISIS Neutronand MuonSource. 2019;DOI: 10.5286/ISIS.E.100529645.[37] Liptak A, Burca G, Kelleher J, Ovtchinnikov E, Maresca J, Horner A. Developments towardsBragg edge imaging on the IMAT beamline at the ISIS Pulsed Neutron and Muon Source: BEAnsoftware. Journal of Physics Communications. 2019;3(113002). DOI: 10.1088/2399-6528/ab5575.2038] Jørgensen JS, Hansen P, Schmidt S. Sparse image reconstruction in computed tomography. DTUCompute; 2013.[39] Zhong Z, Palenstijn WJ, Adler J, Batenburg KJ. EDS tomographic reconstruction regularized bytotal nuclear variation joined with HAADF-STEM tomography. Ultramicroscopy. 2018;191:34–43.DOI: 10.1016/j.ultramic.2018.04.011.[40] Knoll F, Holler M, Koesters T, Otazo R, Bredies K, Sodickson DK. Joint MR-PET reconstructionusing a multi-channel image regularizer. IEEE transactions on medical imaging. 2016;36(1):1–16.DOI: 10.1109/tmi.2016.2564989.[41] Schloegl M, Holler M, Schwarzl A, Bredies K, Stollberger R. Infimal convolution of total general-ized variation functionals for dynamic MRI. Magnetic resonance in medicine. 2017;78(1):142–155.DOI: 10.1002/mrm.26352.[42] Holler M, Huber R, Knoll F. Coupled regularization with multiple data discrepancies. Inverseproblems. 2018;34(8):084003. DOI: 10.1088/1361-6420/aac539.[43] Huber R, Haberfehlner G, Holler M, Kothleitner G, Bredies K. Total generalized variationregularization for multi-modal electron tomography. nanoscale. 2019;11(12):5617–5632. DOI:10.1039/C8NR09058K.[44] Chatnuntawech I, Martin A, Bilgic B, Setsompop K, Adalsteinsson E, Schiavi E. Vectorial to-tal generalized variation for accelerated multi-channel multi-contrast MRI. Magnetic resonanceimaging. 2016;34(8):1161–1170. DOI: 10.1016/j.mri.2016.05.014.[45] Beck A, Teboulle M. A fast iterative shrinkage-thresholding algorithm for linear inverse problems.SIAM journal on imaging sciences. 2009;2(1):183–202. DOI: 10.1137/080716542.[46] Chambolle A, Pock T. A first-order primal-dual algorithm for convex problems with applications toimaging. Journal of mathematical imaging and vision. 2011;40(1):120–145. DOI: 10.1007/s10851-010-0251-1.[47] Van Aarle W, Palenstijn WJ, Cant J, Janssens E, Bleichrodt F, Dabravolski A, et al. Fast andflexible X-ray tomography using the ASTRA toolbox. Optics express. 2016;24(22):25129–25147.DOI: 10.1364/OE.24.025129.[48] van Aarle W, Palenstijn WJ, De Beenhouwer J, Altantzis T, Bals S, Batenburg KJ, et al. TheASTRA Toolbox: A platform for advanced algorithm development in electron tomography. Ul-tramicroscopy. 2015;157:35–47. DOI: 10.1016/j.ultramic.2015.05.002.[49] Palenstijn WJ, Batenburg KJ, Sijbers J. Performance improvements for iterative electron to-mography reconstruction using graphics processing units (GPUs). Journal of Structural Biology.2011;176(2):250–253. DOI: 10.1016/j.jsb.2011.07.017.[50] Kazantsev D, Pasca E, Turner MJ, Withers PJ. CCPi-regularisation toolkit for computed to-mographic image reconstruction with proximal splitting algorithms. SoftwareX. 2019;9:317–323.DOI: 10.1016/j.softx.2019.04.003.[51] Tremsin AS, McPhate JB, Kockelmann WA, Vallerga JV, Siegmund OH, Feller WB. Energy-resolving neutron transmission radiography at the ISIS pulsed spallation source with a high-resolution neutron counting detector. In: 2008 IEEE Nuclear Science Symposium ConferenceRecord. IEEE; 2008. p. 2902–2908. DOI: 10.1109/NSSMIC.2008.4774973.[52] Roessl E, Proksa R. K-edge imaging in x-ray computed tomography using multi-bin photoncounting detectors. Physics in Medicine & Biology. 2007;52(15):4679. DOI: 10.1088/0031-9155/52/15/020. 2153] Schlomka J, Roessl E, Dorscheid R, Dill S, Martens G, Istel T, et al. Experimental feasibility ofmulti-energy photon-counting K-edge imaging in pre-clinical computed tomography. Physics inMedicine & Biology. 2008;53(15):4031. DOI: 10.1088/0031-9155/53/15/002.[54] Barber RF, Sidky EY, Schmidt TG, Pan X. An algorithm for constrained one-step inversionof spectral CT data. Physics in Medicine & Biology. 2016;61(10):3784. DOI: 10.1088/0031-9155/61/10/3784.[55] Firsching M, Niederlohner D, Michel T, Anton G. Quantitative Material Reconstruction inCT with Spectroscopic X-ray Pixel Detectors–a Simulation Study. In: 2006 IEEE NuclearScience Symposium Conference Record. vol. 4. IEEE; 2006. p. 2257–2259. DOI: 10.1109/NSS-MIC.2006.354363.[56] Xie B, Su T, Kaftandjian V, Niu P, Yang F, Robini M, et al. Material decomposition in X-rayspectral CT using multiple constraints in image domain. Journal of Nondestructive Evaluation.2019;38(1):16. DOI: 10.1007/s10921-018-0551-8.[57] Liu X, Yu L, Primak AN, McCollough CH. Quantitative imaging of element composition and massfraction using dual-energy CT: Three-material decomposition. Medical physics. 2009;36(5):1602–1609. DOI: 10.1118/1.3097632.[58] Ronaldson JP, Zainon R, Scott NJA, Gieseg SP, Butler AP, Butler PH, et al. Toward quantifyingthe composition of soft tissues by spectral CT with Medipix3. Medical physics. 2012;39(11):6847–6857. DOI: 10.1118/1.4760773.[59] Grant M, Boyd S, Ye Y. CVX: Matlab software for disciplined convex programming, version 2.0beta; 2008. http://cvxr.com/cvx .[60] Grant M, Boyd S. Graph implementations for nonsmooth convex programs. In: Blondel V, BoydS, Kimura H, editors. Recent Advances in Learning and Control. Lecture Notes in Control andInformation Sciences. Springer-Verlag Limited; 2008. p. 95–110. http://stanford.edu/~boyd/graph_dcp.htmlhttp://stanford.edu/~boyd/graph_dcp.html