Materials Separation via the Matrix Method Employing Energy-Discriminating X-ray Detection

Abstract

The majority of lab-based X-ray sources are polychromatic and are not easily tunable, which can make the 3D quantitative analysis of multi-component samples challenging. The lack of effective materials separation when using conventional X-ray tube sources has motivated the development of a number of potential solutions including the application of dual-energy X-ray computed tomography (CT) as well as the use of X-ray filters. Here, we demonstrate the simultaneous decomposition of two low-density materials via inversion of the linear attenuation matrices using data from the energy-discriminating PiXirad detector. A key application for this method is soft-tissue differentiation which is widely used in biological and medical imaging. We assess the effectiveness of this approach using both simulation and experiment noting that none of the materials investigated here incorporate any contrast enhancing agents. By exploiting the energy discriminating properties of the detector, narrow energy bands are created resulting in multiple quasi-monochromatic images being formed using a broadband polychromatic source. Optimization of the key parameters for materials separation is first demonstrated in simulation followed by experimental validation using a phantom test sample in 2D and a small-animal model in 3D.

1. Introduction

Polychromatic laboratory-based X-ray sources are used for imaging experiments on a wide range of life and materials science samples [1,2]. Since its initial discovery, X-ray imaging has undergone a process of continuous improvement in terms of data collection speed, as well as spatial resolution [3]. A key area of current research is the separation of composite materials into their discrete components based on the acquisition of a single image at each projection angle [4]. Differentiating materials into their constituent components is particularly challenging for weakly interacting biological samples [5,6], which typically do not produce strong transmission features such as absorption edges. Techniques such as X-ray phase contrast microscopy have been employed to help differentiate biological components, but these are often only applicable with coherent, highly monochromatic sources [7,8], particularly when very high spatial resolution is required [9].

However, a pathway towards the clinical application of phase contrast imaging using more compact laboratory-based X-ray sources has emerged in recent years [10]. The availability of high-resolution, energy-discriminating, X-ray detectors that enable discretization of the incident polychromatic signal, offer new opportunities for developing approaches that allow for the effective separation of materials [11,12]. By energy thresholding the polychromatic transmission image at the detector, multiple quasi-monochromatic images can be simultaneously collected [13]. Each image corresponds to a different narrow bandwidth without suffering from any spectral overlap [14,15,16]. Here, we combine energy discriminating X-ray detection with the application of a matrix-based approach to materials separation using the linear attenuation coefficient. We evaluate the performance of this approach for samples composed of two constituent materials, which together form a composite sample using X-ray CT data collected using an energy-discriminating detector.

Elemental contrast or mapping the spatial distribution of an element within samples can be readily achieved using tunable monochromatic X-ray sources, which are typically provided by synchrotrons. For example, by measuring the characteristic peaks associated with X-ray fluorescence spectra, specific elements can be mapped, and their spatial location and concentration determined [17,18,19]. This approach typically involves scanning a small X-ray probe across the sample to construct a 2D or 3D elemental distribution map. The requirements of this method include the need for tightly focused X-ray beams in order to achieve high spatial resolution and a sufficient amount of X-ray flux that data can be collected within a practical time frame. Another method for obtaining the distribution of the specific elements which can be used with a full-field imaging setup is Absorption-Edge Contrast Imaging (AECI). AECI exploits the discrete jumps (absorption edges), which are present in the transmission spectra associated with different elements as a function of the incident X-ray energy [20,21]. The AECI method is based on the subtraction of images collected either side of the absorption edge of the element of interest [21]. Hence, for AECI to work effectively, the tunability of the monochromatic source is critical.

For laboratory-based X-ray sources, which typically have a polychromatic spectrum, elemental mapping is more challenging; however, a number of solutions have been proposed to address this. Methods for elemental mapping that have been successfully demonstrated using lab-based X-rays include dual energy CT [22,23], Ross filter pairs [2] and energy-thresholding using energy-discriminating detectors [13,24]. Dual-energy CT has previously been used to separate bone and soft tissue by exploiting differences in their attenuation properties as a function of energy [25]. The basic principle of dual-energy CT involves collecting data using two distinct energy spectra peaked around different incident X-ray energies [25,26]. In the majority of cases, both sets of spectra have some degree of overlap which can impact the effectiveness of materials separation. Dual-energy CT has previously been used for materials separation in the absence of an absorption edge. By creating an in vitro ‘stone library’ comprising 37 stones of known chemical composition and attenuation ratios, separation and chemical characterisation of kidney stones in twenty seven patients was achieved via the use of a multi-layer X-ray detector [27]. Furthermore, dual energy CT has also been used to measure mass density [28] and to extract materials specific density for two different types of tissues [26]. Whilst dual-energy CT has been successfully applied to soft tissue differentiation, some of the known disadvantages include limitations in terms of flexibility and the potential need for beam-hardening correction [29].

Another approach that was recently demonstrated for both two and three dimensional X-ray imaging uses Ross filter pairs to obtain elemental maps of the sample with a polychromatic source [2]. The basic principle of this method relies on the selection of pairs of filters which have an absorption edge located above and below the absorption edge of the elemental of interest. The Ross filter pair method for elemental contrast provides a very well-defined and sharp energy window which enables elemental enhancement for materials with an absorption edge within the defined energy range. In addition, X-ray fluorescence computed tomography (XFCT) was also used with laboratory-based X-ray sources to characterize the different elements within the sample or tissue [30]. For biological materials this approach generally provides excellent image contrast by incorporating contrast-enhancing elements (often in the form of nanoparticles) such as Mo and Au. The main limitation associated with XFCT in the context of the present work is the need to use these additional agents and ensure they perfuse homogeneously throughout the sample. Other issues include the contamination of the CT images due to Compton scattering and the longer data acquisition times required to raster scan a small pencil beam across the sample, which is particularly challenging for X-ray CT. Some of these issues have been addressed for example a recent study conducted by I. O. Romero et al. with the pencil-beam geometry showed improvements in the quality and speed of XFCT using a ring of spectrometers to detect the fluorescent signal [31]. Recent studies of elemental mapping using polychromatic laboratory sources include the use of polycapillary lenses in a confocal geometry [32,33]. The authors of this work were able to achieve both 2D and 3D elemental mapping via colour tomography. Their method enabled depth profile analysis for a variety of sample types and had excellent sensitivity due to the simultaneous use of two separate detectors optimized to measure the X-ray spectrum independently at both high and low energy.

Another approach to improving elemental selectivity using polychromatic lab-based X-ray sources is to use energy-discriminating X-ray detectors. These types of detectors have undergone rapid development and provide a high degree of spectral filtering via energy thresholds which can be easily defined in software. Detectors capable of full energy resolution imaging with a limited degree of pixelation first became commercially available well over a decade ago, however, the limited number of pixels and large pitch currently associated with these types of detectors [34] limits their application to high-resolution X-ray µCT. A large number of different types of energy-discriminating detector are now commercially available and could potentially be used with the matrix method presented here, including the Medipix [35], Pilatus [36], XPAD [37], PiXirad [15], HEXITEC detector [38] and more recently, the hyperspectral X-ray detector [39]. Hyperspectral X-ray detectors are able to provide a full X-ray spectrum at each individual detector pixel, typically by separating and classifying the incident X-ray photon energy into a very large number (up to several thousand) of discrete energy bins. A recent example of a hyperspectral detector used for X-ray imaging is the SLcam [40].

For a comparison of the performance and characteristics of these types of detector the reader is referred to the recent review of this field by K. Taguchi et al. [41]. Early work using energy-discriminating detectors to generate quasi-monochromatic X-ray images employed sensors which comprised a relatively small number of pixels with a large pixel pitch [42,43]. Through continuous development higher resolution energy discriminating detectors have now become available commercially. These detectors permit the specification of multiple energy thresholds resulting in the simultaneous production of multiple quasi-monochromatic images, each corresponding to a different characteristic energy range with a pixel resolution in the micron range [44]. One of the first demonstrations of material separation using energy-discriminating X-ray detectors was published by X. Wang et al. [11]. In this paper, the authors used specific X-ray contrast enhancing agents in the sample and based on their K-edge signature were able to obtain separate images for each of the different elements. Similarly, P.M. Shikhaliev [45] also used contrast enhancing agents including iodine and CaCO₃ to facilitate soft tissue enhancement in Computed Tomography (CT) using energy-resolving X-ray detectors in a fan-beam geometry. In this work, the author used a total of 5 discrete energy bands to enable separation of the various elements within a phantom test sample. Another study using energy discriminating detectors was also conducted by A.S. Gogolev et al. using a semiconductor matrix direct conversion detector and optimising the energy range in CT to enable differentiation between aluminum and a potassium iodide solution within a phantom test sample [46].

In addition, to enabling elemental separation based on the characteristic K-edge absorption profile of different elements [13] energy-discriminating X-ray detectors also provide a simple and effective means of eliminating or reducing beam-hardening artifacts [24,47]. In a recent study, a high-resolution energy-discriminating detector (PiXirad) was used in X-ray CT to separate soft-tissue from bone, without employing any contrast enhancing agents. The tomographic reconstructions were then quantitatively analyzed enabling bone densitometry using a polychromatic lab-based X-ray source [24]. In summary, the increasing availability of energy-resolving X-ray detection has driven a significant amount of interest in exploiting these detectors for materials differentiation.

In the present study, we have investigated the use of the energy-discriminating PiXirad detector for materials separation in the context where there are no K-edges or contrast enhancing agents. Here we present the first application of the matrix method using multiple quasi-monochromatic images generated with each single exposure to solve the linear attenuation matrix equation resulting in the three-dimensional separation of the sample components. We find that when combining this approach with energy-resolving X-ray detection the method can be generalized to apply to any polychromatic X-ray source and samples where no defined absorption edge is present in the transmission spectrum and no contrast enhancing agent is used. In order to demonstrate this, a sample consisting of a mouse hand is imaged and the 3D distributions of extracted bone and soft tissue.

2. Materials and Methods

2.1. The PiXirad Detector

The PiXirad detector system (PiXirad IC srl c/o INFN Pisa, Italy) was used for these measurements. The detector consists of a CdTe solid state sensor connected to a complementary metal-oxide semiconductor (CMOS) Application-Specific Integrated Circuit (ASIC) readout [16]. The CMOS ASIC has 512 $\times$ 476 pixels arranged on a hexagonal pitch of 60 µm $\times$ 52 μm. The sensor thickness for the PiXirad detector is 650 µm with a corresponding detection efficiency of 100% at 10 keV and 98% at 50 keV [48]. Each pixel includes an electrode, connected to a charge amplifier, which provides two discriminators and two 15-bit counters per register. This enables each pixel to discriminate photons in terms of their energy. The X-ray photons are converted to electric charges at the semiconductor, which are collected and transferred to the ASIC where the signal per pixel is electronically treated. This technology enables the PiXirad detector to generate a maximum of four colour images per single exposure, by setting four individual energy thresholds. For each energy threshold, the PiXirad detector counts photons with energies greater than the threshold energy up to the cut off X-ray energy of the spectrum. The four-energy thresholds produced by the PiXirad detector (shown in Figure 1 as $E_1$, $E_2$, $E_3$ and $E_4$) are used to produce two discrete energy bandwidths ($\Delta E_a$ and $\Delta E_b$), which only contain photon energies within the specified energy range. The discrete energy bandwidths, $\Delta E_a$ and $\Delta E_b$ are approximated here as two quasi-monochromatic beams each having a bandwidth (Δ𝐸/𝐸) less than 34%. This bandwidth can readily be defined using the PiXirad whose spectral energy resolution is shown in Figure S2 in the Supplementary Material.

A key consideration here is that the performance of the PiXirad detector can be affected by charge-sharing. Charge-sharing is an artefact which arises when a single photon is absorbed and produces a signal that is ‘shared’ across two or more detector pixels. The spread of the charge cloud during drift or cloud splitting arising from X-ray fluorescence generally leads to a reduction in the measured energy of the detected photons [14]. In pixelated energy-resolving X-ray detectors, charge-sharing can impact the majority of detection events and thus needed to be taken into account in this study. The characteristics of charge-sharing include the appearance of low energy tails in the observed photopeaks when using monochromatic X-rays which, particularly for CdTe and CdZnTe detectors, can be due to carrier trapping, which results in ‘hole tailing’ [49]. However, we note that in the PiXirad (and other high-resolution energy-resolving X-ray detectors) the effect of hole tailing is ameliorated due to the small pixel size of the detector. There are several interconnected phenomena that determine the degree and characteristics of charge sharing when using polychromatic X-ray sources [41]. The most obvious of these are the apparent increase in the number of detected counts within the low energy region of the spectrum and concomitant reduction in the number of counts within the higher energy region [12].

2.2. Matrix Equation for Materials Separation Using the PiXirad Detector

The method for separating two materials distributions requires two quasi-monochromatic intensity images to be defined using four energy thresholds ($E_1,E_2,E_3$ and $E_4$) without any spectral overlap. In order to separate a larger number of materials, additional energy thresholds are required which necessitates more than one exposure per projection [13]. The transmitted intensity image corresponding to the first energy bandwidth, $I(\Delta E_a)$, is created via the subtraction of the measured intensities corresponding to the two energy thresholds $I(E_1)$ and $I(E_2)$. This subtraction can be represented as

$$\begin{gathered} I(\Delta E_a)=I\left(E_1\right)-I\left(E_2\right) \\ =\underset{E_1}{\overset{E_c}{{\displaystyle \int }}}S\left(E\right) \exp \left({{\displaystyle \sum }}_{j=1}^l-{\mu }_j\left(E\right).s_j\right)dE-\underset{E_2}{\overset{E_c}{{\displaystyle \int }}}S\left(E\right) \exp \left({{\displaystyle \sum }}_{j=1}^l-{\mu }_j\left(E\right).s_j\right)dE \\ =\underset{E_1}{\overset{E_2}{{\displaystyle \int }}}S\left(E\right) \exp \left({{\displaystyle \sum }}_{j=1}^l-{\mu }_j\left(E\right).s_j\right)dE \\ =\exp \left({{\displaystyle \sum }}_{j=1}^l-{{\displaystyle \int }}_{E_1}^{E_2}S\left(E\right){\mu }_j\left(E\right).dE.{\mathrm{s}}_j\right) \\ \approx \exp \left({{\displaystyle \sum }}_{j=1}^l-{\mu }_j\left(\Delta E_a\right).{\mathrm{s}}_j\right) \end{gathered}$$

(1)

where $S\left(E\right)$ is the spectral density of the incident X-rays, $E_c$ is the cut-off energy of the X-ray spectrum, ${\mu }_j\left(E\right)$ is the linear attenuation coefficient of material $j$, and $s_j$ is the thickness of material $j$. We note that Equations (1) and (2) are not normalized by the total intensity and that $I(\Delta E_a)$ and $I(\Delta E_a)$ are a function of the integrated flux between the relevant energy thresholds.

Similarly, the second energy bandwidth, $I(\Delta E_b),$ defined by the subtraction between the two other measured energy threshold intensities, $I(E_3)$ and $I(E_4)$, is given by

$$\begin{gathered} I(\Delta E_b)={{\displaystyle \int }}_{E_3}^{E_4}S\left(E\right) \exp \left({{\displaystyle \sum }}_{j=1}^l-{\mu }_j\left(E\right).s_j\right)dE=\exp \left({{\displaystyle \sum }}_{j=1}^l-{{\displaystyle \int }}_{E_3}^{E_4}S\left(E\right){\mu }_j\left(E\right).dE.{\mathrm{s}}_j\right) \\ \approx \exp ({{\displaystyle \sum }}_{j=1}^l-{\mu }_j\left(\Delta E_b\right).s_j) \end{gathered}$$

(2)

By defining narrow energy bandwidths, we assume that we can approximate two quasi-monochromatic X-ray beams represented by $\Delta E_a$ and $\Delta E_b$ for the first and second energy bands respectively as shown in Equations (1) and (2). Note that the effective linear attenuation coefficient at $\Delta E_a$ and $\Delta E_b$ can be defined as

$${\mu }_j\left(\Delta E_a\right)={\int }_{E_1}^{E_2}{\mu }_j\left(E\right).S\left(E\right).dE , {\mu }_j\left(\Delta E_b\right)={\int }_{E_3}^{E_4}{\mu }_j\left(E\right).S\left(E\right).dE$$

(3)

where the integral of S(E) should be unity within the relevant energy range.

We consider two different materials present within the sample denoted as material 1 and material 2, with ${\mu }_1\left(E\right)$ and ${\mu }_2\left(E\right)$ denoting their respective linear attenuation coefficient, and with $s_1$ and $s_2$ respectively their sample thicknesses. Taking the intensity of the incident beam to be unity, we can express $I(\Delta E_a)$ and $I(\Delta E_b)$ as

$$I(\Delta E_a)=\exp \left(-{\mu }_1\left(\Delta E_a\right).s_1-{\mu }_2\left(\Delta E_a\right).s_2\right)$$

(4)

and

$$I(\Delta E_b)=\exp \left(-{\mu }_1\left(\Delta E_b\right).s_1-{\mu }_2\left(\Delta E_b\right).s_2\right)$$

(5)

These equations can be written in matrix form as:

$$\left[\begin{array}{c}\log (I(\Delta E_a))\\ \log \left(I(\Delta E_b\right))\end{array}\right]=\left[\begin{array}{cc}-{\mu }_1(\Delta E_a)& -{\mu }_2(\Delta E_a)\\ -{\mu }_1(\Delta E_b)& -{\mu }_2(\Delta E_b)\end{array}\right]\left[\begin{array}{c}{s}_1\\ s_2\end{array}\right]$$

(6)

Inversion of Equation (6) yields the final expression for the thickness of each of the components:

$$\left[\begin{array}{c}{s}_1\\ s_2\end{array}\right]=\frac{1}{det}\left[\begin{array}{cc}-{\mu }_2(\Delta E_b)& {\mu }_2(\Delta E_a)\\ {\mu }_1(\Delta E_b)& -{\mu }_1(\Delta E_a)\end{array}\right]\left[\begin{array}{c}\log (I(\Delta E_a))\\ \log \left(I(\Delta E_b\right))\end{array}\right]$$

(7)

where the determinant ($det$) is given by: $({\mu }_1(\Delta E_a)\times {\mu }_2(\Delta E_b))-\left({\mu }_2(\Delta E_a\right)\times {\mu }_1(\Delta E_b))$ and cannot be vanishingly small. Equation (7) shows with a 2 $\times$ 2 matrix, two intensities at two different energies are required to separate two materials.

2.3. Simulation Study

A simulation study was conducted using the measured incident polychromatic X-ray spectrum of the lab source shown in Figure 1. The simulation study assumed an energy sampling frequency of 0.2 keV and did not take into account the actual finite experimental energy resolution of the PiXirad detector (Figure S2). For the simulation, the MATLAB (The MathWorks Inc., Sydney, Australia) programming package was used to model the experiment. The simulation includes the effect of the energy threshold, the spectral distribution of the source, information about the linear attenuation as function of energy, as well as experimental geometry and sample size and thickness. However, we note that the simulation study did not take into account the effect of charge-sharing due to the complex nature of this phenomenon for polychromatic X-rays. The transmitted intensity is calculated via a summation over all energies present in the spectrum weighted by their spectral contribution. A schematic diagram providing an overview of the simulation is presented in Figure S3 in the Supplementary Material.

In the simulation, three different phantom samples were considered, shown in Figure 2. Each phantom consisted of a 10-mm diameter cylinder of epoxy resin (density 1.13 g/cm³) with three different inserts of the same resin enriched with varying amounts of hydroxyapatite (HA), Ca₅(PO₄)₃(OH), with densities of 1.39, 1.65 and 1.90 g/cm³, corresponding to HA400, HA800 and HA1200 respectively, see

Table 1. Composition and density for the components of the phantom samples (used for both simulation and experimental studies). The manufacturer for the experimental phantom was QRM GmbH, Germany [51].

Material	Density ± 0.02 (g/cm3)	Base (Epoxy Resin) (%)	HA Content (%)
epoxy	1.13	100	0
HA400	1.39	71.31	28.87
HA800	1.65	51.2	48.8
HA1200	1.9	36.51	63.49

. The diameter of all the inserts was 2 mm. The absorption coefficients for the elements comprising the phantom as a function of X-ray energy were obtained from the NIST materials database [50]. The polychromatic intensity in the simulations was calculated via a summation of the monochromatic intensities over all incident wavelengths, weighted by their relative contribution to the measured incident X-ray spectrum, S(E). $s_j$ is the one-dimensional projected thickness of the phantom (as viewed from the XZ plane).

Figure 3 summarizes the workflow for the simulations applying the matrix method to separate epoxy and HA400 (using the first phantom shown in Figure 2a). The top row of Figure 3 shows the four transmitted intensity projections corresponding to four energy thresholds. The second row shows the projection intensity difference at $I(\Delta E_a)$ which is obtained from Equation (1) and at $I(\Delta E_b)$ which is obtained from Equation (2). The third row shows the normalized intensity difference for $\overline{I}(\Delta E_a)$ and $\overline{I}(\Delta E_b)$ such that the value in the absence of a sample is equal to unity. These results form the input for Equation (7) (the matrix equation). The linear attenuation parameters for epoxy and HA400 in the matrix equation at energies $\Delta E_a$ and $\Delta E_b$ were calculated using Equation (3). The bottom row shows the retrieved thickness of the separated epoxy and HA400 components respectively based on the solution of the matrix equation.

In order to investigate the optimum conditions required to separate the epoxy and HA, three different energy bandwidths for energy threshold positions ranging from 16 keV to 30 keV were investigated as shown in

Table 2. Simulation parameters used for the four energies thresholds $E_1$, $E_2$, $E_3$ and $E_4$ to obtain three pairs of different energy bandwidths.

Threshold Position (keV)	Bandwidth (keV)
	$\Delta E$ = 3	$\Delta E$ = 6	$\Delta E$ = 9
$E_1$	16	16	16
$E_2$	19	22	25
$E_3$	27	24	21
$E_4$	30	30	30

. Note that the two bandwidths at 9 keV are overlapping. All three simulated phantom samples were analyzed. The matrix method, employing four intensity images, enabled the separation of epoxy and HA400, epoxy and HA800, and epoxy and HA1200. An example result of the separation of epoxy from HA400 is shown in Figure 4, together with the actual thickness value based on the input sample parameters for comparison.

The difference between the resulting retrieved thickness, $s_{retrieved}$, and the actual thickness, $s_{actual}$, is evaluated using:

$${\chi }^2=\frac{{{\displaystyle \sum }}_{p=1}^n{\left(s_{p,actual}-s_{p,retrieved}\right)}^2}{{{\displaystyle \sum }}_{p=1}^n{\left(s_{p,actual}\right)}^2}$$

(8)

where $p$ denotes the pixel number $p=$1, 2,… $n$. The calculated ${\chi }^2$ error metric for all data sets is presented in Figure 5 in the absence of noise (a, c, e) and with 1% noise included (b, d, f). The noise was modelled by adding randomly generated values with a standard deviation of 1% about the mean intensity value. Figure 5 shows that in the absence of noise, the smallest bandwidth of 3 keV results in the most effective separation (smallest ${\chi }^2$ value) of HA from epoxy for all material densities. This is due to the fact that as the bandwidth increases the assumption that Equations (4) and (5) represent quasi-monochromatic intensity images which can be represented by a single effective energy becomes increasingly invalid decreasing the effectiveness of the materials separation. Moreover, the 9 keV bandwidth was found to negatively impact the retrieval owing to the larger bandwidth and subsequent breakdown of the assumption of quasi-monochromaticity as well as the spectral overlap.

Figure 5 also indicates that, under the conditions simulated here, the most effective materials separation occurs for the lowest density material (HA400) and epoxy with the highest density material (HA1200) yielding the least effective separation. However, when 1% noise is added (consistent with the typical noise level observed in the experiment) the lowest ${\chi }^2$ value occurs for the 6 keV bandwidth. This result likely reflects the trade-off between monochromaticity (which is best for the 3 keV bandwidth) and signal-to-noise (SNR), which is highest for the 9 keV bandwidth.

3. Experimental Results

Experimental X-ray imaging data were collected using the Xradia© micro XCT200 (Carl Zeiss X-ray Microscopy, Inc., Pleasanton, CA, USA) polychromatic X-ray laboratory source combined with the PiXirad detector system, as shown in Figure 6. X-ray spectra were recorded using the XR-100T-CdTe AMPTEK detector designed to characterize X-ray tubes and for X-ray spectroscopy applications. The source was a closed X-ray tube (Hamamatsu) consisting of a tungsten target and operated at a tube voltage of 40 kVp at a power of 10 W. The PiXirad detector was cooled down to an operating temperature of −20 °C using a combination of water cooling and nitrogen gas. A Galil motion controller was added to the system and was connected between the PiXirad detector and the rotation stage in order to synchronize their communication during the tomography scan. The effective image pixel size is defined by the physical detector pixel size divided by the geometric magnification factor, M, where $M=\frac{R_1+R_2}{R_1}$. Here, $R_1$ is the source-to-sample distance, which is equal to 125 cm, and $R_2$ is the sample-to-detector distance equal to 115 cm. In this experiment M = 1.9, giving an effective pixel size of 30 μm. The phantom sample (QRM GmbH) consisted of a 10 mm diameter cylinder of epoxy resin with three inserts (of 2 mm diameter each) of HA400, HA800 and HA1200. The phantom sample composition in the experiment was identical to that used in the simulation.

The PiXirad detector system simultaneously outputs 4 images, each corresponding to a different energy threshold. Whilst there is significant flexibility over the choice of energy threshold, due to the finite energy resolution of the detector they cannot be chosen arbitrarily and are restricted to certain ranges as determined by the PiXirad control software [16]. When setting energy thresholds in the control software for the PiXirad detector the software takes into account the finite energy resolution of the detector and automatically limits the minimum separation allowed between thresholds. The energy resolution values (and hence minimum separation) for each of the thresholds used in this paper are shown in Figure S2. Based on the possible energy thresholds, the following threshold values were used for the experiment: $E_1=16 \mathrm{keV}$ and $E_2=21.1 \mathrm{keV}$, and $E_3=$ 23.9 keV and $E_4=29.7 \mathrm{keV}$ resulting in energy bandwidths of 5.1 keV and 5.8 keV, respectively. Note that the small discrepancy between the experimental energy threshold values and the simulated threshold values is well below the energy resolution of the PiXirad. Experimental 2D projection data was first collected from a well-characterized phantom sample comprising a HA400 insert embedded in epoxy (identical to Figure 2a). The detector exposure time was set to 3.5 s which was determined by the dynamic range of the detector at the lowest energy threshold (where the flux was greatest). The total time for the X-ray µCT experiment was 2 h including sample preparation and alignment; with 45 min for actual data collection which was required in order to build up sufficient statistics for analysis. Correction of the measured intensity to account for the effect of charge-sharing was applied according to Equation (S1) (see Supplementary Material) and the procedure summarised in Figure 3 was used to separate the HA400 insert from the epoxy. The experimental results are shown in Figure 7a,b. Note that in Figure 7b transmission through half of the resin matrix was blocked due to the presence of additional inserts (HA800 and HA1200) within the phantom sample. The application of this mask results in the sudden drop in transmission but does not otherwise impact the separation of the insert. The tube voltage was then increased to 80 kVp to retrieve the denser (HA800, epoxy and HA1200, epoxy) material inserts as shown in Figure 7c–f. The thickness of the material pairs (HA400, epoxy), (HA800, epoxy) and (HA1200, epoxy) retrieved using the combined energy-resolved detector and matrix method are observed to be in very good agreement with the known component thicknesses (supplied by the manufacturer). The calculated χ² error inside the range indicated by the black arrows can be seen on each plot.

Following a successful proof-of-concept with the phantom, the matrix approach using the PiXirad detector was applied to 3D material separation of a ‘real-world’ biological sample comprising a mouse hand (with dimensions of 5 mm (H) × 14 mm (V)). The sample was provided by Dr. Philippe Bouillet (based at the Walter and Eliza Hall Institute of Medical Research) [52]. The sample was stored in formalin but removed from solution during the measurement and placed inside a plastic container to prevent dehydration. The same experimental conditions (e.g., energy thresholds and exposure time) were used for tomographic imaging of the mouse hand. In the present case the attenuation coefficient values used for bone were ${\mu }_{bone}\left(\Delta E_a\right)=$ 1109 m⁻¹ and ${\mu }_{bone}\left(\Delta E_b\right)=$ 392 m⁻¹ and the attenuation coefficient values for soft tissue are ${\mu }_{tiss}\left(\Delta E_a\right)=$ 102 m⁻¹ and ${\mu }_{tiss}\left(\Delta E_b\right)=$ 47 m⁻¹. These attenuation coefficients were derived from the X-ray properties of HA400, which has similar X-ray attenuation characteristics to mouse bone [53] and water, which is commonly used to represent hydrated soft tissue. The variation of these coefficients as a function of the X-ray energy were obtained from the National Institute of Standards and Technology (NIST) database [50]. In the case that the X-ray attenuation properties of the sample components are not known a priori, they can be measured in a separate experiment using a coupon of known thickness and measurements of the incident and transmitted X-ray intensity [54]. The attenuation coefficients within the energy ranges of $\Delta E_a$ and $\Delta E_b$ can then be calculated using Equation (3). The HA400 calibration data was used to correct for the effects of charge-sharing (see Equation (S1) and Figure S1 in Supplementary Material). For the tomography scan, 360 projections were collected by rotating the sample in 0.5 degrees steps between −90° and +90°. An example of the corrected projection showing the four energy threshold images is presented in Figure 8a–d. The matrix equation method using the PiXirad detector yielded the results shown in Figure 7f and Figure 8e (with 2 × 2 binning applied). This step was repeated for each projection resulting in two tomographic datasets, one for the soft tissue and the other for the bone. A filtered back projection (FBP) reconstruction algorithm together with ring artefact correction was used to obtain the final 3D reconstruction (Figure 9).

The 3D reconstruction obtained without applying the matrix method to the energy resolved PiXirad data to separate the components is also shown for comparison. Figure 9 shows that the matrix approach is able to effectively separate bone from soft tissue using a lab based polychromatic X-ray source and just a single exposure per projection. However, the SNR in the retrieved images is degraded compared to the individual images corresponding to each of the four energy thresholds. The reduced SNR arises from the image subtraction implicit in solving Equation (7), i.e., noise will add in quadrature under following (weighted) image subtraction, with the net effect being a reduction in SNR [55]. This is a common artifact of matrix inversion and could be addressed using the ‘brute force’ approach of improving the statistics by collecting a larger number of images per projection (leading to longer data collection times). However, another common method to address this problem is matrix regularization which is able to reduce errors introduced by the inversion process by introducing modifications into the matrix operator that bias particular solutions with desirable properties. For example, Tikhonov regularization is commonly used to introduce high-pass operators to enforce smoothness if, for example, the underlying solution is known to be continuously varying [56].

In simulations of the experiment, we applied increasing amounts of experimental noise and where able to obtain an acceptable result up to a Poisson noise level of 6% (which is greater than expected for the experiment). Another limitation arises from energy-resolving X-ray detectors in general, which, to a greater or lesser degree depending on the sensor type, exhibit some charge-sharing effects. The quality of the final result for the retrieved sample components thus depends both on the SNR and the characterisation of the charge-sharing characteristics of the specific detector used. A number of factors contribute to the SNR of the data and the degree of charge-sharing artifacts. These include, the X-ray tube flux, the linear absorption coefficient for the sample, and the overall sample thickness. Here, we study a sample of just 5 mm × 14 mm in dimensions, extending this method to e.g., patient limbs or even whole animals would require a source capable of delivering higher energy X-rays and most likely the implementation of additional matrix regularisation to maintain an acceptable result whilst minimizing the absorbed dose which would lead to lower SNR data.

4. Conclusions

We have described and demonstrated an approach to material separation based on solution of the linear attenuation matrix equations which can be applied to polychromatic laboratory sources using commercially available energy resolving X-ray detectors. The approach presented here enables a more general application of this method by using the energy thresholding capabilities of an energy-discriminating (PiXirad) detector to separate multiple components. This concept was explored in both simulation and experiments for a phantom sample in 2D and a ‘real-world’ biological sample in 3D. In order to separate two materials two energy bandwidths, need to be defined based on four intensity images in order to obtain two quasi-monochromatic images. Using the PiXirad, this can be achieved in a single exposure. However, as we have shown previously, the number of elements to be separated can be increased be taking multiple images with a different set of energy thresholds defined per exposure [13]. The ability to carry out this type of experiment is now widely available via a variety of different commercial X-ray energy discriminating detectors. Based on the results of this study the matrix method for materials separation combined with commercial energy-discriminating X-ray detection appears to be a promising avenue of research with a wide range of potential applications particularly in the field of medical imaging.

An important question is how materials with very similar densities could be separated using the matrix method combined with energy-discriminating X-ray detection. An interesting application of this would be the differentiation of glandular and adipose regions in breast tissue or other samples that contain very similar density materials. The current state-of-the-art for research using energy discriminating detectors for soft tissue differentiation is to address the challenge of separating materials of comparable bulk density but differing μ(E) dependence which is not amenable to analysis by traditional methods. In the present case, the key limiting factor is likely be the signal to noise ratio in the data and whether data of sufficient quality can be collected in order to reliably separate materials with almost identical X-ray interaction characteristics. Moreover, in this current setup, options for improving the SNR include increasing the number of frames per projection, binning together pixels at the cost of reduced spatial resolution, or by using an alternate energy-resolving detector with a higher dynamic range (e.g., 32 bit). In order to explore this application, future work will focus on the further optimisation of the energy ranges and bandwidths used for these measurements which, in the case of energy-discriminating detectors, can be easily varied by the user.