Enhanced hyperspectral tomography for bioimaging by spatiospectral reconstruction

Here we apply hyperspectral bright field imaging to collect computed tomographic images with excellent energy resolution (~ 1 keV), applying it for the first time to map the distribution of stain in a fixed biological sample through its characteristic K-edge. Conventionally, because the photons detected at each pixel are distributed across as many as 200 energy channels, energy-selective images are characterised by low count-rates and poor signal-to-noise ratio. This means high X-ray exposures, long scan times and high doses are required to image unique spectral markers. Here, we achieve high quality energy-dispersive tomograms from low dose, noisy datasets using a dedicated iterative reconstruction algorithm. This exploits the spatial smoothness and inter-channel structural correlation in the spectral domain using two carefully chosen regularisation terms. For a multi-phase phantom, a reduction in scan time of 36 times is demonstrated. Spectral analysis methods including K-edge subtraction and absorption step-size fitting are evaluated for an ex vivo, single (iodine)-stained biological sample, where low chemical concentration and inhomogeneous distribution can affect soft tissue segmentation and visualisation. The reconstruction algorithms are available through the open-source Core Imaging Library. Taken together, these tools offer new capabilities for visualisation and elemental mapping, with promising applications for multiply-stained biological specimens.

With the increasing development in X-ray detector technology, the interest in energy-selective tomography has grown in recent years, particularly for medical imaging. In conventional X-ray absorption computed tomography (CT), each detector pixel records the total number of detected photons, irrespective of their energy, building up a single radiograph for each projection angle. Contrast is hence solely provided by differences in attenuation based on a sample's local material composition. Issues may then arise in post-processing segmentation, particularly for phases or structures of similar electron density, such as different types of soft tissue. In addition, due to the use of polychromatic radiation in laboratory X-ray imaging, beam-hardening artefacts are common, as the nonlinear nature of attenuation as a function of energy is ignored without energy discrimination 1 . The introduction of spectroscopic detectors have enabled an additional dimension of information to be acquired, by measuring both the energy and position of each incident photon. Given that every element has a unique attenuation profile, these may be used as spectral 'fingerprints' in energy-sensitive CT imaging, such that full 3D elemental mapping may be performed from a single CT scan.
Recent spectral detectors take two main forms, differing mainly in their measurement processes, and number of energy 'channels' used for photon binning. Multispectral detectors use a set of 4-8 threshold (reference) channels, to which each incident photon is allocated based on the electrical signal generated upon detection. Such detectors have found use in the field of medical imaging for soft tissue differentiation [2][3][4][5][6] . While the low channel number enables high count-rates similar to conventional CT imaging 7 , multispectral imaging provides coarse energy resolution (~ 5-10 keV), thus for spectrally similar species, the threshold positions require specific pre-selection, as well as a priori knowledge of sample composition. In contrast, hyperspectral detectors can achieve very fine energy resolution (~ 1 keV), with the ability to store photons into hundreds of narrow energy channels. The result is that, for every pixel, we acquire a 'pixel spectrum' , containing a full absorption profile, for

Results
Firstly we examine a simple aluminium phantom test sample, measuring 0.5 cm in diameter, and 1 cm in height. Three cylindrical holes (diameter 0.7 mm) were drilled to a sample depth of 75% (0.75 cm), and filled with a different powder (CeO 2 , ZnO and Fe). Comprising different elements, each powder has different absorption properties for analysis. The phantom was imaged using a micro-focus source at 60 kV with beam power 6 W. Two scans were acquired using different exposure times, with an effective voxel size of 98 μm. The full dimensions of the sample were captured within the 80 × 80 detector field of view (FOV). The number of channels used in hyperspectral imaging are determined by two system characteristics-the tube voltage and the energy-channel calibration. For the 60 kV tube voltage applied, the full spectral profile was segmented across 200 energy channels for each pixel, with a measured single energy channel width of 1.2 keV. Given we have a priori knowledge of sample composition and, thus, absorption edge position in this case, a reduced subset of channels were selected for both datasets (corresponding to channel numbers 100-200, or energy range 28-56 keV), in order to increase the speed of reconstruction. To evaluate the capabilities of the iterative method in handling low count and under-sampled datasets, a comparison was made between two acquisition schemes: Scan A was taken with 180 projections and 180 s exposure time, and Scan B was taken with 30 projections and 30 s exposure time. Through prior testing, it was found that the conditions used for Scan A produced sufficiently low noise and strong feature definition, such that it may be treated as a 'ground truth' , optimum state, when using conventional channel-wise www.nature.com/scientificreports/ FDK reconstruction. All results obtained under the conditions of Scan B were then compared directly to this reference case. Given the nature of spectral datasets, each channel may be analysed on an individual basis, akin to a set of multiple monochromatic scans. Figure 1a illustrates a set of sinograms for Scan A, corresponding to four equally spaced energy channels, combined with the corresponding central slice FDK reconstructions showing changes in attenuation value, μ as a function of energy. An important factor to note, however, is that we cannot directly identify chemicals through matching the absolute value of μ in each voxel to theoretical values for different elements. This is because non-linear effects, such as ' charge-sharing' between pixels, can lead to erroneous photon energies (and number of events) during scanning. For instances of lower X-ray flux, chargesharing may be corrected for, as detailed elsewhere 16 . In this study, we show that, even without correction, we may still perform chemical identification and mapping through the identification of absorption edges, which remain clearly visible at their theoretical energy positions. Through examination of the monochromatic reconstructions in Fig. 1b (covering the energy range 28-49 keV), we observe a significant rise in attenuation for one such phase, indicative of a spectrally significant elemental marker. By analysing voxel absorption spectra, chemical insight on sample composition can be obtained by examining changes in attenuation as a function of energy. Three regions of interest (ROIs) were selected in Fig. 1b covering each powder phase. The spectra shown in Fig. 1c corresponds to phase 1 where, as expected, an absorption edge is observed at 40.4 keV. Owing to the high energy resolution of the system, one may conclusively Figure 1. Attenuation variation in the spectral dimension for a multi-phase powder phantom. (a) Set of four sinograms taken between channels 100-175 at equal spacing of 25 channels (corresponding average energies are shown-channel width 1.2 keV) for Scan A. A single discontinuity in the sinograms appears due to an interrupted scan. (b) The corresponding FDK reconstructions for each energy channel of the sinograms in (a). Three distinct regions are observed, corresponding to the three metal-based powders. The colour scale measures the attenuation, and is consistent across both images. Three ROIs are highlighted (white/blue squares, marked by numbers for each respective powder phase). (c) Average voxel spectra of powder phase ROI 1. A line signifying the theoretical position of the cerium K-edge is overlaid for comparison. (d) Comparison of measured absorption spectra (top) for ROIs 2 and 3 located in the zinc oxide and iron phases respectively, and the theoretical values (bottom) over the same spectral range. www.nature.com/scientificreports/ confirm that such an edge corresponds to that of cerium (40.443 keV). For the two remaining material phases in Fig. 1d, little variation in attenuation is observed. In addition, no absorption edges belonging to these materials (ZnO, Fe) were seen, as they fall below the sensitivity of the hyperspectral system (< 10 keV). Other methods may be employed to segment such materials, for instance through measuring the relative change in attenuation as a function of energy, as has been shown elsewhere 16 . Here we note that the rate of change of μ with increasing energy does not precisely match the theoretical predictions. This may be due to effects such as pulse pile-up, which can distort the spectral profile due to the miscounting of coincident photon events 28 . A full investigation of the detector system electronics and signal acquisition properties, as a function of flux and energy range, is required to identify the precise cause. However, such detailed analysis is beyond the scope of this paper, with the observations not affecting results in terms of K-edge imaging. Unsurprisingly, the conditions of Scan A provide high quality reconstructions, with sharp feature definition, and little noise fluctuations across the full spectral range. Next we explore differences when reconstructing the low count and undersampled dataset of Scan B. As previously discussed, analytic reconstruction methods, such as FDK, often fail to adequately reconstruct features for low SNR data. Here we explore the extent to which noisy, few projection data can be reconstructed, through the use of a spatiospectral reconstruction algorithm. The reconstruction problem is typically formulated as a combination of a data fitting term between measured and reconstructed data, and one or more regularisation terms which encode desirable image characteristics. The problem is then solved using an iterative optimisation algorithm. CIL provides a number of building blocks to formulate the optimisation problems, and solvers to find a numerical solution. In this case we used a combination of two regularisation terms, known as the Total Variation (TV) and Total Generalised Variation (TGV), implemented along the spatial and spectral dimensions respectively. Further details of the regularisation terms, and how they were applied to the dataset, can be found in the methods section. The joint regularisation method, referred to here onwards as TV-TGV, is a novel method, chosen based on the prior knowledge that we expect noisy images, combined with the presence of an absorption edge. The CIL software enables such a method to be constructed and applied, and its advantages over other state-of-the-art methods has been demonstrated elsewhere in the case of Bragg-edge neutron tomography 29 .
For the reconstruction of Scan B, TV-TGV was applied simultaneously across the full set of energy channels. By doing so, we exploit the correlations between neighbouring channels, compared to FDK, which is applied channel-by-channel and therefore cannot make use of such structural similarities. As described in Eq. 1 (see methods section), three regularisation parameters ( α for TV, β 1,2 for TGV) were optimised. Final parameters were determined to be 0.002, 0.18 and 0.25 for α , β 1 and β 2 respectively. A range of iteration numbers were tested, with no discernible improvements in image quality observed beyond 1000 iterations. Final reconstruction of the full volume was achieved with a runtime of 25 minutes for 1000 iterations of the algorithm. Figure 2a shows a comparison of FDK and TV-TGV, through direct comparisons of reconstructed slices in the transverse and frontal planes. From the FDK reconstruction for Scan A, a distinct variation in attenuation is observed across the CeO 2 powder region, suggesting inhomogeneity in the powder. The frontal view confirms the non-uniform distribution, illustrating a consistently higher attenuation in the outer perimeter throughout the full depth of the sample. The slight ' cupping' feature is also reminiscent of beam hardening encountered for white beam imaging, which was not expected because of the high (1.2 keV) energy discrimination in this case, and is discussed further in the Supplementary Information (Fig. S3). Both reconstructed views illustrate the sufficiently high count and sampling of the dataset, with minimal noise and strong feature edge definition. In contrast, the undersampled reconstructions of Scan B highlight the limitations of FDK for low projection number and exposure time. Significant streak artefacts, projecting outwards from the CeO 2 powder phase, are observed, attributed to photon starvation due to the reduced photon counts stored in each energy channel, under the faster scan conditions. These artefacts are in addition to any further streaks caused by undersampling, combined with increased noise across the reconstructed slices.
Upon application of the TV-TGV method for Scan B, however, the data is recovered well, with the use of TV in the spatial domain providing enough smoothing to remove all streak artefacts, while maintaining edge preservation, restoring an image quality comparable to that of Scan A. We quantify this noise suppression through calculation of the contrast-to-noise ratio (CNR). Here, we apply the method employed by Bian et al., where the standard deviations of both the signal and background ROIs are taken into account 30 . In addition, we average out random variations both in the spatial and spectral domains through two steps: larger ROIs covering a greater range of pixels (5x5) average out spatial variations, while spectral fluctuations are reduced by computing the average standard deviation value for each ROI across all channels, and applying this to the CNR calculation for each individual channel. We use the ROIs marked in Fig. 2a (ROI S, ROI Bg) and calculate the CNR between the ZnO and Al phases in every channel. The result, shown in Fig. 3a, emphasises the successful feature restoration of the TV-TGV algorithm, with higher CNR measured in all channels, compared to both the Scan A and Scan B FDK reconstructions. As expected, with increasing energy, CNR gradually declines, with the difference in CNR values between TV-TGV and FDK Scan B narrowing, but still maintaining an average improvement of approximately 390% over the energy range. Additional quantitative analysis is shown in Fig. 3b in the form of channelwise root-mean-square-error (RMSE) measurements. Calculations are made by quantifying the difference in pixel values between ROI 1 of the cerium phase in Scan A, and the same ROI in each respective Scan B reconstruction. The residuals were then squared, averaged and square-rooted to determine an RMSE value for the energy channel. The method was then applied across all channels. The significant reduction in image error following TV-TGV reconstruction is highlighted, with RMSE values consistently lower than the FDK Scan B equivalent, across the full energy range. An additional single channel example of the RMSE calculation is shown in Supplementary Fig. S4 to map spatial errors across different image planes. The application of TV-TGV can be analysed further by comparing spatial and spectral profiles for each method, as shown in Figure 2b and 2c respectively. For the spatial profile, measured across two of the powder phases (Fig. 2b), the advantages of the TV-TGV method for low count data are evident. Noise levels are significantly reduced for Scan B, compared to www.nature.com/scientificreports/ the equivalent FDK profile, with the TV-TGV method almost completely replicating the profile acquired for Scan A. For the spectral profile ( Fig. 2c), TGV regularisation improves the linearity of the regions either side of the cerium K-edge for Scan B. It should be noted that, for such a simple phantom, identification of the absorption edge is not an issue here, regardless of reconstruction algorithm. However, the results demonstrate the capability of handling noisy, few projection data, and achieving high image quality in fast scan acquisitions. In this case, a reduction in scan time by a factor of 36 has been shown. We next explore a more realistic sample in the field of biological imaging, where low SNR may pose problems in material identification and feature segmentation. An iodine-stained lizard head (Anolis sp.), measuring approximately 17 mm in length, and 10 mm in width, was scanned at a beam voltage of 50 kV, at a maximum power of 0.7 W, reconstructed with a voxel size of 137 μm. A spectral subset of channels from 60 to 160 (energy range 17.3-45.0 keV) was chosen given the prior knowledge of the single iodine stain, eliminating the need to www.nature.com/scientificreports/ reconstruct the full spectral range for identification of key elemental signals. A reduced projection dataset was used, such that only 60 projections (at 120 s exposures) were reconstructed over the full 360° rotation, to test the capabilities of the reconstruction algorithms. The need to limit X-ray dose is more prevalent in the life sciences field and, as such, the ability to extract high quality information from low count and/or undersampled data is crucial. Once more, TV-TGV reconstruction was applied, with the results compared to the conventional FDK method. Here, initial parameter estimates were based on the optimal conditions achieved in the reconstruction of the phantom sample, allowing faster identification of the optimal parameters. Final reconstruction was performed with parameters of 0.002, 0.25 and 0.35 for α , β 1 and β 2 respectively, with a total runtime of 40 minutes following 1000 iterations. Figure 4 demonstrates the smoothing effects in both the spatial and spectral domains, with the lens and jaw adductor muscles highlighted as key regions of interest. For both the axial and sagittal reconstructed slices (Fig. 4a), spatial smoothing provided improved edge definition of the exterior head shape, as well as the outer eye and jaw regions. The presence of the iodine contrast agent is easily confirmed with the precise matching of the theoretical edge position (33.169 keV) in each spectral profile (Fig. 4b). The benefit of noise suppression in the energy domain is demonstrated, allowing the presence of absorption edges to be more clearly defined, where they can often be lost within noisier reconstructions like FDK. Areas of lower iodine uptake, such as that of the jaw adductor muscles, may therefore be confidently identified as iodine-containing structures. As we have no ground truth reconstruction available for RMSE evaluation, we use only CNR calculations for a quantitative image metric, validated by the strong alignment between the RMSE and CNR metrics for the powder phantom. Figure 4c shows the channelwise CNR calculation for ROIs marked in Fig. 4a (ROI S, ROI Bg), which cover the jaw muscle and background respectively. Increased contrast using TV-TGV is found in every channel over the reconstructed range, compared to FDK. The impact of the iodine K-edge is also clear in its ability to sharply improve image contrast, and as such CNR, relative to the background material. The impact of TV-TGV is significant, producing an average CNR improvement of approximately 300% compared to the standard FDK algorithm.
Given the presence of an absorption edge, spectral analyses may be performed to provide information on iodine distribution throughout the biological specimen. Firstly we utilise the availability of a spectral profile in each voxel to measure relative iodine concentration, by virtue of spectral profile fitting. As shown in Fig. 5, www.nature.com/scientificreports/ linear least squares fitting was applied to regions before and after the absorption edge step for both the FDK and TV-TGV reconstructed volumes. By extrapolating and evaluating these fits at the known position of the K-edge (33.169 keV), we can precisely measure the size of the step change, �µ 0 16 . Repeating this process at every voxel provides us with a map of �µ 0 , which is directly proportional to the concentration of iodine present in the sample. Calculated values of �µ 0 are shown for both reconstructed volumes in Fig. 5a. Significant noise distortions in the FDK spectra lead to erroneous linear fitting, and consequently inaccurate measurements of �µ 0 , as shown in Fig. 5b. Spectral smoothing due to TV-TGV, however, ensures improved precision in calculation of relative iodine concentration across the volume. Results indicate the diffusion of iodine fully into the lens, with high concentrations at the interior, and slightly lower levels on the exterior surface. Further, the increased reliability of �µ 0 measurements allows us to confidently identify 'hot spots' of higher iodine uptake, appearing on the brain and sections of the jaw muscles. The results are in good agreement with expected uptake regions of iodine contrast agent 11,12 . Our second analysis uses the absorption edge as a means of segmenting iodine-containing material from the remaining structures. For this, we used K-edge subtraction (KES). That is, spectral information is extracted from energy channels before the edge of interest, and subtracted from an equivalent set after the edge. The result is a dataset containing only the contrasting material, eliminating other structures where attenuation is slight across this energy range. The method of KES has previously been applied for monochromatic imaging either side of absorption edges for segmentation 31,32 , as well as in hyperspectral cases, highlighting its potential for segmenting materials where more than one K-edge is present 16 . A detailed description of the method is provided in Supplementary Information (see Fig. S5). Here, KES also offers an opportunity to evaluate reconstruction quality by direct comparison of tissue segmentation for both the FDK and TV-TGV methods. www.nature.com/scientificreports/ In order to measure the success of correct tissue segmentation, we match our hyperspectral reconstructed volumes against a DECT scan acquired of the same sample, reduced to the same spatial resolution (137 μm) as that of the hyperspectral data. DECT has long been regarded as the 'gold standard' of biological stain imaging, and thus works well both as a measure of where hyperspectral X-ray CT stands in comparison, as well as an evaluation tool for each case of our spectral KES method. Segmented views of the sample are shown in Fig. 6 for both the TV-TGV regularised method, as well as the FDK reconstructed volume. The resulting 3D visualisations are shown upon hyperspectral KES around the iodine edge. Using the DECT segmentation as reference for identifying key soft tissue features, the advantages of TV-TGV over FDK become clear. The increased level of noise due to FDK leads to reduced visibility, particularly in regions of lower iodine concentration, such as the tongue and jaw adductor muscles. In contrast, following KES of the TV-TGV volume, clear separation is observed for regions to which the iodine has diffused. Structures including the brain, lens, tongue and jaw muscle all show strong X-ray signal enhancement due to sufficient staining by elemental iodine. Segmentation of the remaining material offers the ability to observe 'non-contrast-enhanced' structures. In this case, the external skull structure, consisting mostly of hydroxyapatite (HA), remains. As confirmed by the DECT results, visualisation of certain bone structures, including the skull roof and mandible region, are achieved. Full definition of the skull structure is lost however, and this is attributed to the long-term storage of the sample prior to imaging. As such, bone mineral has accumulated in some regions, while having dissipated in others. Therefore, precise segmentation of HA material was not expected. Nevertheless, the advantage of combined hyperspectral imaging and advanced reconstruction algorithms is provided through the successful segmentation of iodine in the biological structure. Moreover, while 60 projections with 120 s exposures were taken, results from our phantom sample suggest a further reduction in exposure time is possible with minimal loss in reconstruction quality, owing to the application of regularised algorithms.

Discussion
The above case studies have highlighted the current state of the art for hyperspectral imaging in a lab-based setting and shown that, given an appropriate reconstruction algorithm, comparable levels of feature definition and characterisation of tissues and structures can be obtained to those in dual-energy and multispectral CT. However, our work is not without limitations. Spatial resolution is still far behind the standard set in conventional, and dual-energy, CT, emphasised by the downsampling of the aforementioned DECT results from their initial 9 μm voxel size. While currently limiting in the range of samples and features that may be observed, spatial resolution is only expected to improve with new detector technology. The use of an iodine contrast agent did not fully test the capability of chemical detection in hyperspectral imaging, given its high affinity for multiple soft tissue regions. A reasonable extension to the experiment is the use of additional chemical tracers, which bind to specific biological structures at varying concentrations. Given the high energy resolution of the hyperspectral (HEXITEC) detector, visualisation and segmentation of multiply-stained structures may be easily performed by virtue of measuring several spectral fingerprints. It has also previously been shown that absolute values of concentration may be determined in spectral imaging, given the use of an appropriate calibration phantom 33 , www.nature.com/scientificreports/ offering potential diagnostic insight on tracer concentration and distribution as a function of time. The use of an ex vivo biological sample here may also be used as a step towards in vivo structure analysis in the future. Biological scans performed in vivo are limited more severely in terms of allowable X-ray dose and chemical concentration, hence the analysis performed in this work offers a solution in the form of regularised reconstruction for noisy, short exposure cases. Alternatively, there is a case to be made for the applicability of hyperspectral imaging in the absence of distinct spectral markers. With the availability of full spectral profiles at each pixel, beam-hardening artefacts may be eliminated now that changes in attenuation may be discriminated as a function of energy. In addition, while methods like K-edge subtraction are no longer applicable, measurement of relative attenuation changes over the energy range enables segmentation of poorly contrasting materials, which otherwise cannot be differentiated in conventional X-ray CT 16 . Previous work on bone densitometry and soft tissue segmentation has been shown to be possible with few channel spectroscopic imaging, without the need for contrast agents 34 . The use of advanced iterative reconstructions, in particular with application of regularisation terms, is also in its early stage development. The chosen method here of TV-TGV in the spatial-spectral domains is applicable given the composite materials producing sharp absorption edges, however optimisation of the regularised parameters for each dimension can be a slow process, due to the decoupled nature of the reconstruction algorithm. While every individual dataset will require specific tuning of each parameter, further development in the reconstruction protocols, combined with software packages like CIL to enable them, continues to allow these processes to be performed more quickly, while also laying a foundation upon which future methods may be based. In addition, the above studies have demonstrated the reduction in overall scan time enabled by the algorithms in CIL, offsetting the time taken to optimise reconstructed image quality. We predict that these methods will continue to allow further reductions in scan time, as well as limiting overall X-ray dose, without suffering significant losses in image quality, an advantage that may prove to be hugely beneficial to imaging in the biological field.

Conclusion
This paper has applied hyperspectral imaging to a simple test phantom, as well as, for the first time, mapping the location of staining in a biological sample, using as many as two hundred energy channels (~ 1 keV resolution). We have highlighted the vulnerabilities of conventional reconstruction methods for lab-based hyperspectral X-ray imaging. Through the use of a novel, spatiospectral reconstruction algorithm, we have enabled precise chemical identification and mapping at the micrometer scale. Examination of a multi-phase phantom emphasised the significant reductions (36 times shorter) in scan time achievable by implementing regularised reconstruction to compensate for noisy datasets. In performing an ex vivo spectral CT scan of an iodine-stained lizard head sample, we have shown the capability of hyperspectral CT to have the elemental sensitivity to compete with existing techniques, such as DECT, in soft tissue segmentation and structural analysis, but with definitive identification of the iodine location through its characteristic K-edge. While here a single stain was measured and visualised, the exploration of multi-labelled biological samples is possible, given the high spectral resolution of www.nature.com/scientificreports/ the detector. Further, the weaknesses of analytic regimes such as FDK have been highlighted for spectral imaging, particularly for short scan acquisitions, reinforcing the need for standardised, iterative algorithms such as those provided in CIL 24 . Together with the reduced scan times they enable, the correlated reconstruction methods open up the potential for hyperspectral studies in fields including non-destructive testing, security scanning and chemical catalysis. With improving detector technology and multi-staining methods, we conclude that lab-based hyperspectral CT offers great future prospects for biological research, among a number of other fields, such as chemical engineering, geology, materials science and cultural heritage.

Materials and methods
Phantom sample preparation. An aluminium cylinder was used as the matrix for three internal powders, due to its low attenuation relative to other metals. The choice of CeO 2 offered a clear example of a characteristic spectral marker (absorption edge), while the remaining two materials were used to evaluate the ability to restore lower attenuating structures following severe noise distortion. In addition, they enable discussion to be made on the identification of phases, absent of spectral markers.
Biological sample preparation. For this experiment the head of a lizard (Anolis sp.) was scanned and analysed. No live animals were used as part of the study. The sample was purchased as a fixed/preserved specimen from Nasco Education (USA). For long-term stability, the sample was fixed in formalin and stored in 70% ethanol, prior to staining. The sample was then dehydrated to 100% ethanol, before staining with 1% elemental iodine in absolute ethanol ( I 2 E ). It has previously been shown that I 2 E offers strong contrast in soft tissue, allowing for discrimination from bone and teeth (hydroxyapatite) structures 10,11 . After staining, the sample was washed with 100% ethanol. Finally the biological specimen was mounted in 1.5% Agarose.
X-ray detector. The hyperspectral imaging was performed using an energy-sensitive HEXITEC detector 25 , consisting of a 1 mm thick CdTe single-crystal semiconductor, bump-bonded to an ASIC producing a 2 cm × 2 cm detection area. The system is split into an 80 × 80 pixel array, with a 250 μm pitch. The detector offers an energy resolution of up to 800 eV at 59.5 keV and 1.5 keV at 141 keV. All raw data was acquired on an event-byevent basis using the HEXITEC detector software.

Data acquisition routines.
For imaging of the phantom sample, a parallel source-sample-detector configuration was implemented in the custom-built Colour Bay, part of the Henry Moseley X-ray Imaging Facility (HMXIF) at The University of Manchester. The walk-in X-ray bay contains a 225 kV source, and full manipulator control is available via MATLAB scripts. Two scans of the phantom were acquired. Identical scanning conditions were implemented, at a geometrical magnification of 2.54, with the polychromatic X-ray source operating at a tube voltage of 60 kV with beam power 6 W. The chosen parameters ensured sufficient counts for energies close to measured absorption edges, optimised image contrast, while remaining within the count-rate limit of the detector. Further, X-ray flux was kept as low as possible to minimise the issue of photon saturation for the detector (X-ray source profile shown in Supplementary Information, Fig. S1). Exposure time, however, was varied, from 30 s up to 180 s per projection. In total, 180 projections were acquired per scan, with 2° step size over a full 360° rotation, resulting in total scan times of 2.5 and 11 hours, accounting for buffer times and bias voltage refreshing between projections in single photon detection 25 . Prior to reconstruction, one dataset was downsampled to 30 projections with 30 s exposure time to evaluate reconstruction quality on undersampled data. As a further demonstration of the flexibility in lab-based hyperspectral imaging, the HEXITEC was next combined with the Nikon High Flux Bay system at the HMXIF. Once more a 225 kV source was utilised, however improved contrast was possible due to lower flux capabilities, necessary for imaging of soft tissue within the biological sample. Source and sample manipulation were controlled by Nikon's proprietary software Inspect-X. A software module was created using the IPC interface to Inspect-X 35,36 , enabling communication between Nikon software and the spectral detector. The biological sample was secured to the rotation stage, held in place such that the specimen was suspended vertically during the acquisition. A geometric magnification of 1.81 was determined to project the sample fully in the detection area. The polychromatic X-ray beam was operated at a peak voltage of 50 kV, at a maximum power of 0.7 W (X-ray source profile shown in Supplementary Information, Fig. S1). Projection images were recorded at an angular step size of 2° over a full rotation, with exposure times of 120 s for each of the 180 projections, for an 8 hour total scan time. A reduced subset of the dataset was later taken for reconstruction, as described in the main text. For all scans throughout the study, sets of four flat-field projections were acquired both before and after scanning for fixed-pattern noise subtraction, while a further dark current correction was applied during set-up, to minimise reduced spectral response due to increased leakage current, typically found at the edges of the detector 37 . Detected events were binned into a set of spectral channels, with the number determined by the maximum X-ray energy.
Prior to scanning, a preliminary energy calibration procedure was performed through the measurement of characteristic fluorescence signals from a series of metals 38 . The calibration procedure provides a direct transition between energy channels and their corresponding energies. Further, an inter-pixel gain correction was applied through the use of a correlative optimised warping algorithm using the same data 39,40 . A calibration dataset was acquired before scanning of each sample, with the energy resolution of the system determined by measuring the distribution of FWHM values for characteristic peaks (see Supplementary Information, Fig. S2). Calibration values were consistent across each sample dataset. In each case, the energy resolution of the detector at 59.5 keV was found to be 1.21 ± 0.40 keV.
The accompanying DECT scan was conducted using a Zeiss Xradia XRM-400. The dual-energy acquisition performed scans at 40 kV and 80 kV, with a 0. 17  Spectral CT data reconstruction. Initial processing of the acquired datasets was handled via MATLAB routines. In the case of the biological sample, a combined wavelet-based Fourier filter was applied in every channel for the removal of ring artefacts across the dataset 41 . For each dataset, the resulting 4D matrix (3 spatial dimensions, 1 spectral dimension) was then reconstructed using the CIL software. Full functionality and operation of the Python-based framework has been discussed elsewhere 23,24 . Given its traditional use in conventional cone-beam reconstruction, the FDK method formed the baseline from which all other results were compared and evaluated. For comparison, an iterative reconstruction algorithm was chosen, combining data fitting with spatial and spectral regularisation terms. The iterative algorithm took the form: where our first term concerns classic least-squares data fitting, with hyperspectral projection data, b, related to voxel value, u through the operator, A, based on the system geometry and properties. The latter terms concern the addition of spatial and spectral regularisation. Since 4D datasets exhibit different image properties in the spatial and spectral dimensions, two different regularisation terms were applied for their noise reduction capabilities. The result is a ' decoupled' regularisation algorithm. We define TV x,y,z (u) as the Total Variation (TV), applied across each spatial dimension 42 for a single channel, before summing over all channels (channel-wise). TV is one of the most widely used regularisers, as it favours piece-wise constant images, with sharp edge boundaries. The application of TV therefore allows for noise suppression of flat signal regions, while maintaining the local discontinuities at structure edges. The TV model fits well for CT images and has been successfully used previously for undersampled CT data 43 , as well as noise suppression purposes 42,44 . The use of TV regularisation has also previously demonstrated its benefits in spectral image reconstruction over analytic methods like FDK 21,22 . However, TV is known for introducing 'staircasing' artefacts for piece-wise affine or smooth signals, for example ramp structures, which can result in patchy, unnatural reconstructed images 44 . Further, as the aim of TV regularisation is to reduce unwanted signal variations such as noise, the final reconstruction can suffer from a loss of contrast due to reductions in intensity 44 . For the spectral domain, we define TGV c (u) as the Total Generalised Variation (TGV), based on the method proposed by Bredies et al. 45 . Here we apply TGV regularisation along the channel direction, for each individual voxel in the image. TGV becomes applicable in the case of K-edge imaging, where we expect a smoothed step shape to our energy profiles across the edge position, spanning a few channels. TGV is able to reduce noise variations while maintaining the definition of any absorption edges present. Further, TGV avoids the staircasing effects experienced with TV. CIL provides a number of iterative algorithms to solve Equation 1, here we applied the primal dual hybrid gradient (PDHG) method 46 . Three regularisation parameters ( α for TV, β 1,2 for TGV) were used to control the strength of penalisation for each regulariser, and must be optimised for each dataset. Therefore parameter values were chosen to suppress noise while preserving edge definition of reconstructed images and spectral profiles. In the case of noise suppression, optimal parameters vary based on the noise level in the dataset. Values of α > 1 have previously been used for images of significantly lower SNR 47 , while for TGV regularization, selected values of the β 1,2 parameters are often based on a ratio between the two, with similar work on denoising suggesting this ratio, β 2 β 1 be either √ 2 or 2 [47][48][49] . These factors were therefore used to determine limits for initial parameter estimates. For the phantom sample, optimal values were determined based on three factors: visual comparison with reconstructed slices of Scan A, precision of absorption edge position, and minimum obtained 'primal-dual gap' . The latter is a quantitative measure of convergence for the PDHG method, and has previously been used as a stopping criterion for the iterative method 46,50 .

Data availability
The phantom sample datasets generated and analysed during the current study are available from https:// doi. www.nature.com/scientificreports/