Polyquant CT: direct electron and mass density reconstruction from a single polyenergetic source

TL;DR

This paper introduces a novel segmentation-free, data-driven method for directly reconstructing electron and mass density from single polyenergetic CT scans, improving accuracy without prior segmentation or restrictive models.

Contribution

It presents a new statistical reconstruction algorithm that directly fits attenuation to density, bypassing the need for prior segmentation or physical model assumptions.

Findings

Higher accuracy in density quantification demonstrated

Superior quantitative imaging shown with numerical and real data

Method maintains similar complexity to existing iterative algorithms

Abstract

Quantifying material mass and electron density from computed tomography (CT) reconstructions can be highly valuable in certain medical practices, such as radiation therapy planning. However, uniquely parameterising the X-ray attenuation in terms of mass or electron density is an ill-posed problem when a single polyenergetic source is used with a spectrally indiscriminate detector. Existing approaches to single source polyenergetic modelling often impose consistency with a physical model, such as water--bone or photoelectric--Compton decompositions, which will either require detailed prior segmentation or restrictive energy dependencies, and may require further calibration to the quantity of interest. In this work, we introduce a data centric approach to fitting the attenuation with piecewise-linear functions directly to mass or electron density, and present a segmentation-free statistical reconstruction algorithm for exploiting it, with the same order of complexity as other iterative methods. We show how this allows both higher accuracy in attenuation modelling, and demonstrate its superior quantitative imaging, with numerical chest and metal implant data, and validate it with real cone-beam CT measurements.

Figures

Equations

$$I_{\mathrm{out}}=I_{\mathrm{in}}\exp \left(-{\int }_{\ell }\mu (\ell )d\ell \right),$$

$$I_{\mathrm{out}}={\int }_{\xi }{I}_{\mathrm{in}}(\xi )\exp \left(-{\int }_{\ell }\mu (\ell ,\xi )d\ell \right)d\xi .$$

$$y_i\sim \mathrm{Poisson}\left\{\sum _{j=1}^{N_{\xi }}{b}_i({\xi }_j)\exp \left(-[\mathbf{\Phi }\mathbf{\mu }({\xi }_j){]}_i\right)+s_i\right\}\text{mbox}fori=1,...,N_{\mathrm{ray}},$$

$$\mu (\xi )=\rho N_g({\sigma }_{\mathrm{pe}}(\xi )+{\sigma }_{\mathrm{incoh}}(\xi )+{\sigma }_{\mathrm{coh}}(\xi )),$$

$$N_g=N_A\sum _i\frac{{\omega }_i{Z}_i}{A_i},$$

$${\rho }_e=\frac{\rho Ng}{{\rho }_{\mathrm{water}}{N}_{g,\mathrm{water}}},$$

$$\text{HU}=1000\frac{\mu -{\mu }_{\mathrm{water}}}{{\mu }_{\mathrm{water}}-{\mu }_{\mathrm{air}}},$$

$$\mu (\xi )=\underset{\text{photoelectric}}{\underbrace{K_1{\rho }_e{Z}_{\mathrm{eff}}^n{\xi }^{-3}}}+\underset{\text{Compton scatter}}{\underbrace{K_2{\rho }_e{f}_{\mathrm{KN}}(\xi )}},$$

$$\mu (\xi )=a_1{\mu }_{\mathrm{water}}(\xi )+a_2{\mu }_{\mathrm{bone}}(\xi )=\rho \left(a_1{m}_{\mathrm{water}}(\xi )+a_2{m}_{\mathrm{bone}}(\xi )\right),$$

$$\hat{\mathbf{\mu }}(\mathbf{x},\xi )=\sum _{i=1}^{N_f}{\mathbf{f}}_i\odot ({\alpha }_i(\xi )\mathbf{x}+{\beta }_i(\xi )),$$

$$f_i=\{\begin{array}{ll}1& \text{if }{k}_{i-1}\le x<k_i\\ 0& \text{otherwise}\end{array}\text{mbox}fork=1,\dots ,N_f,$$

$${\alpha }_i{k}_i+{\beta }_i={\alpha }_{i+1}{k}_i+{\beta }_{i+1}.$$

$$\hat{\mathbf{\mu }}({\mathbf{\rho }}_{\mathbf{e}},\xi )={\mathbf{f}}_1({\mathbf{\rho }}_{\mathbf{e}})\odot {\alpha }_1(\xi ){\mathbf{\rho }}_{\mathbf{e}}+{\mathbf{f}}_2({\mathbf{\rho }}_{\mathbf{e}})\odot [{\alpha }_2(\xi ){\mathbf{\rho }}_{\mathbf{e}}+{\beta }_2(\xi )].$$

$$\hat{\mathbf{\mu }}(\mathbf{x},\xi )=\sum _{i=1}^{N_f}{\mathbf{f}}_i\odot \left(\left[\frac{{\alpha }_i}{{\xi }^3}+{\gamma }_i{f}_{\mathrm{KN}}(\xi )\right]\mathbf{x}+\frac{{\beta }_i}{{\xi }^3}+{\delta }_i{f}_{\mathrm{KN}}(\xi )\right),$$

$$\sum _{j=1}^{N_{\xi }}{b}_i({\xi }_j)\exp \left(-[\mathbf{\Phi }\mathbf{\mu }({\xi }_j){]}_i\right)=\sum _{j=1}^{N_{\xi }}{b}_i({\xi }_j)\exp \left(-[\mathbf{\Phi }\hat{\mathbf{\mu }}(\mathbf{x},{\xi }_j){]}_i\right)\text{mbox}fori=1,...,N_{\mathrm{ray}}.$$

$${\psi }_i(\mathbf{x},\xi )\equiv b_i(\xi )\exp \left(-[\mathbf{\Phi }\hat{\mathbf{\mu }}(\mathbf{x},\xi ){]}_i\right)\text{mbox}fori=1,...,N_{\mathrm{ray}},$$

$$\mathrm{NLL}(\mathbf{x};\mathbf{y})=\sum _{i=1}^{N_{\mathrm{ray}}}\sum _{j=1}^{N_{\xi }}{\psi }_i(\mathbf{x},{\xi }_j)+s_i-y_i\log \left(\sum _{j=1}^{N_{\xi }}{\psi }_i(\mathbf{x},{\xi }_j)+s_i\right),$$

$$\mathbf{d}(\mathbf{x})=\mathbf{y}\oslash \left(\sum _{j=1}^{N_{\xi }}\mathbf{\psi }(\mathbf{x},{\xi }_j)+\mathbf{s}\right)-1,$$

$$\frac{\partial \mathrm{NLL}(\mathbf{x};\mathbf{y})}{\partial \mathbf{x}}\approx \sum _{i=1}^{N_f}{\mathbf{f}}_i(\mathbf{x})\odot {\mathbf{\Phi }}^T\left[\sum _{j=1}^{N_{\xi }}{\alpha }_i({\xi }_j)\mathbf{\psi }(\mathbf{x},{\xi }_j)\odot \mathbf{d}(\mathbf{x})\right],$$

$$\hat{\mathbf{x}}=\underset{\mathbf{x}\in \mathcal{C}}{\mathrm{argmin}}\mathrm{NLL}(\mathbf{x};\mathbf{y})+\lambda R(\mathbf{x}),$$

$$L_0={\Vert {\mathbf{\Phi }}^T\left[\left(\sum _{j=1}^{N_{\xi }}{\alpha }_1({\xi }_j)\mathbf{b}({\xi }_j)\right)\odot \Phi 1\right]\Vert }_{\infty },$$

$${\mathbf{prox}}_{\delta R}(\mathbf{z})=\underset{\mathbf{\rho }\in \mathcal{C}}{\mathrm{argmin}}\frac{1}{2}\Vert \mathbf{z}-\mathbf{\rho }{\Vert }_2^2+\delta \lambda R(\mathbf{\rho }).$$

$${\mathbf{\mu }}_{\mathrm{mono}}=\underset{\mathbf{\mu }}{\mathrm{argmin}}(\mathbf{\Phi }\mathbf{\mu }-\mathbf{l}{)}^T\mathbf{W}(\mathbf{\Phi }\mathbf{\mu }-\mathbf{l})+\lambda R(\mathbf{\mu }),$$

$$\sqrt{\frac{1}{N_{\mathrm{vox}}}\sum _{i=1}^{N_{\mathrm{vox}}}({\hat{\rho }}_i-{\rho }_i}{)}^2,$$

Full text

1 Introduction

Physically quantifying reconstructions from CT is highly valuable in some medical practices. In radiation therapy for example, the electron density inferred from CT images of the patient allows the dose deposition from the treatment beam to be accurately modelled (Battista et al.,, 1980). Additionally, quantifying the bone mineral density allows osteoporosis to be characterised and the risk of bone fractures to be assessed (Marshall et al.,, 1996). Due to the non-linear energy-dependence of X-ray attenuation however, mapping from a set of raw measurements to a consistent physically quantitative reconstruction is not straight forward, and requires both actively accounting for the polyenergetic effects such as beam hardening, and establishing a relation to the quantity of interest.

Mapping from CT to physical density is usually treated in a three step process: linearising the measurements by correcting for scatter, taking the logarithm (Chang et al.,, 2014) and non-linearly calibrating from a polyenergetic to approximate monoenergetic source (Joseph and Spital,, 1978); inverting the linearised projections with analytic or iterative reconstruction algorithms (Fessler,, 2014); then finally applying a non-linear calibration to mass or electron density (Schneider et al.,, 1996).

In general, the mapping from a single polyenergetic to a monoenergetic measurement is an ill-posed problem, since the spectral information of the irradiated material is lost with the standard energy integrating detectors (Curry et al.,, 1990). Whilst imaging the human body however, most tissues may be reasonably modelled with few distinct materials: namely water and bone as in (Joseph and Spital,, 1978). If the amount of each material in a projection is known a priori from a water–bone segmentation, then it is possible to approximately reverse the polyenergetic effects through polynomial fitting or with a look up table. For more accuracy, this may also be brought into the reconstruction model, where (Elbakri and Fessler,, 2002) demonstrate statistical quantitative imaging of mass density, with the water–bone assumption. Given that the composition is itself highly correlated with density, then the prior segmentation may not be necessary, but estimated during reconstruction (Elbakri and Fessler,, 2003).

Another approach is to model the attenuation explicitly in terms of physical processes, given quantitative physical parameters of interest. One such choice is to model the photoelectric and Compton effects in terms of relative atomic number and electron density, which is reasonably accurate for relevant elements and energies (Weber and van den Berge,, 1969; Rutherford et al.,, 1976; Jackson and Hawkes,, 1981). Given measurements from two sufficiently different spectra—a technique known as dual-energy CT (DECT)—a projection of Compton attenuation can be uniquely determined (Alvarez and Macovski,, 1976), of which electron density is an analytic function (Klein and Nishina,, 1929). This DECT technique effectively bypasses the need for assumptions such as water–bone compositions, so should be applicable to a wider range in materials, although one faces practical difficulties in generating the two spectra (Johnson,, 2012). Additionally, one will expect a loss in accuracy between significantly different elemental compositions, such as soft tissues, bone and metallic implants, since the two parameter model is not consistent over a wide range in atomic species Jackson and Hawkes, (1981).

In (De Man et al.,, 2001), the authors introduce an iterative maximum-likelihood polychromatic algorithm for CT (IMPACT), which models the energy independent factors in the photoelectric–Compton model from (Alvarez and Macovski,, 1976) as piecewise-linear functions of monoenergetic attenuation, allowing reconstruction from a single source. This method does not require any prior segmentation, and is reported on a wide range of materials including metallic implants. By imposing the energy dependence of the two parameter model however, the physical consistency will also degrade throughout diverse material types, due to the inconsistency of these parameterisations in effective atomic number and electron density Jackson and Hawkes, (1981).

The second conversion from either reconstructed attenuation in Hounsfield units (HU) in (Joseph and Spital,, 1978; De Man et al.,, 2001) or mass density in (Elbakri and Fessler,, 2002, 2003) to electron density is also an ill-posed problem. This is because X-ray interaction depends on the environment of the electrons as well as their density, and will vary considerably with atomic number for imaging energy ranges (Jackson and Hawkes,, 1981; Curry et al.,, 1990). Although, again given the fact that most human tissues have similar properties, a single piecewise linear fitting is reasonably accurate in practice (Constantinou et al.,, 1992), though it will not be consistent with synthetic materials such as some plastics (Schneider et al.,, 1996).

In both the photoelectric–Compton and a material decomposition such as the water–bone model, there is some degree of fitting to materials and energies of interest. Instead, one could use a purely data-centric approach. Here, given a representative set of substances, one could use a model that accurately parameterises the energy dependent attenuation in terms of the quantity of interest, without necessarily any physical justification. One such method is to model the energy dependent attenuation as a piecewise-linear function of quantitative density, which may be fit to a set of materials of interest. When the transitions between linear sections are independent of energy we also get the nice property that the computation order in an iterative method is independent of the number of energies considered.

In this article, we will study in detail the piecewise linear quantitative model for CT, and will show how it may be incorporated in a regularised iterative reconstruction algorithm. Specifically, this provides a generalised method for directly quantifying the electron or mass density of a heterogeneous specimen, and is also able to model hard metallic structures without any prior segmentation. In (Mason et al.,, 2017) we presented a preliminary study using this idea for the specific case of calculating electron density for radiotherapy planning.

In preparation of this manuscript, we became aware of the commercial method DirectDensity™ from Siemens Healthineers® (Ritter and Mistry,, 2016), which reports direct reconstruction into relative electron density from a single polyenergetic source. This is a preprocessing technique combining bone detection with a projection-based material decomposition similar to (Joseph and Spital,, 1978). In this article we compare against the bone–model of (Elbakri and Fessler,, 2002), which has been shown to itself provide superior performance to that of (Joseph and Spital,, 1978).

1.1 Contributions

We establish and analyse a general modelling technique to allow direct quantitative reconstruction from a single polyenergetic source, where we study the cases of electron and mass density imaging. Unlike existing approaches of fitting to physical parameterisations of attenuation, such as bone–water or photoelectric–Compton, we fit directly to the data, which we demonstrate is more accurate over a wide range of biological tissues. We show how this model may be incorporated into statistical reconstruction, and propose an algorithm for performing this that allows further convex spatial regularisation to be used. By design, the complexity of using our model does not scale with the number of discrete energies, and will have an order of computational cost $2.5\times$ that of standard monoenergetic iterative algorithms. As another consequence of the fitting constraints, one could also quantify the equivalent attenuation from a mono-energetic source—known as ‘quasi-monoenergetic’ in DECT Johnson, (2012)—though this is not evaluated in this work. We also demonstrate how this model may also directly mitigate metal artefacts, without any need for segmentation.

1.2 Article Structure

We begin this article with background material in Section 2 on X-ray attenuation, existing polyenergetic parameterisations, and the probabilistic measurement model we will invoke for reconstruction. We then propose our generalised data centric model in Section 3.1, and outline the cases of electron and mass density quantisation for biological tissues, as well as synthetic materials and metal implants. In Section 3.2 we demonstrate how to utilise this model in statistical reconstruction, and detail one such algorithm in Section 3.2.1. The experimentation in Section 4 evaluates our method with: a model accuracy test in Section 4.1 against other physical parameterisations; a numerical reconstruction test in Section 4.2 with simulated fan-beam CT of a chest and pelvis with metallic hip implants; and reconstruction validation on a real physical phantom scanned with cone-beam CT. We then discuss important considerations in Section 5 leading to conclusions in Section 6.

2 Background

2.1 The CT measurement model

In CT, one is able to observe a specimen’s attenuation through the radiation intensity after transmission. The magnitude of this is found from the Beer–Lambert law, given for a monoenergetic beam as

[TABLE]

where $I_{\mathrm{in}}$ is the incident intensity, $\ell$ is the line-of-sight path of the beam through specimen, $I_{\mathrm{out}}$ is the output intensity one is able to measure. Since in practice, $\mu$ is energy dependent and typically the source is polyenergetic, the output intensity becomes

[TABLE]

For a finite number of photons, the measured intensity will be probabilistic with an approximate Poisson distribution (Chang et al.,, 2014). If we also move the attenuation, measurements and energy spectrum into a discretised setting, we can write the measurement process as

[TABLE]

where $N_{\mathrm{ray}}$ is the number of CT measurements, $N_{\mathrm{\xi}}$ is the number of energy bins, $\boldsymbol{b}(\xi)\in\mathbb{R}^{N_{\mathrm{ray}}}$ is a vector of incident intensities, $\boldsymbol{\mu}(\xi)\in\mathbb{R}^{N_{\mathrm{vox}}}$ is the vector of attenuation coefficients with $N_{\mathrm{vox}}$ the number of voxels, $\boldsymbol{\Phi}\in\mathbb{R}^{N_{\mathrm{ray}}\times N_{\mathrm{vox}}}$ is the system matrix describing the summation along the paths from source through specimen onto each detector, and $\boldsymbol{s}\in\mathbb{R}^{N_{\mathrm{ray}}}$ is the expectation of the scatter or other background noise reaching the detector.

2.2 Material dependent X-ray attenuation

The mechanism that allows various regions in a heterogeneous specimen to be differentiated is their degree of X-ray attenuation. For biological tissues irradiated with a diagnostic X-ray source, the significant phenomena contributing to the attenuation of incident radiation are photoelectric and scattering effects—consisting of Compton, Rayleigh and Thompson scatter (Curry et al.,, 1990). The combined attenuation strength of a given material can be quantified as

[TABLE]

where $\xi$ is the energy of the incident radiation, $\sigma_{\mathrm{pe}}$, $\sigma_{\mathrm{incoh}}$ and $\sigma_{\mathrm{coh}}$ represent the interactive cross sections—quantifying the probability of interaction—of photoelectric, incoherent (Compton) and coherent (Rayleigh and Thompson) effects, $\rho$ is the mass density, and $N_{g}$ is the number of electrons per unit volume defined as

[TABLE]

where $N_{A}$ is Avagadro’s number, and $Z_{i}$, $A_{i}$, $\omega_{i}$ are atomic number, atomic weight and relative fraction by mass of a material’s constituent elements (Schneider et al.,, 1996). A convenient parameter to use is the relative electron density, which is

[TABLE]

where $N_{g,\mathrm{water}}$ is the absolute electron density of water and $\rho_{\mathrm{water}}$ its mass density, having a value of almost exactly $\rho_{\mathrm{water}}=1\,\text{g/cm}^{3}$ at room temperature.

From (4), (5) and (6), one may quantify the X-ray attenuation in terms of the energy independent mass density or relative electron density—$\rho$ and $\rho_{e}$ respectively—given knowledge of the energy dependent interactive cross sections, which may be found from existing parameterisations in Section 2.3 or our proposed model in Section 3.1. A more common approach however, is simply to use calibration curves on the reconstructed images. Using the Hounsfield scale defined as

[TABLE]

examples of calibration curves are shown in Figure 1.

The calibration curves in Figure 1 are generated from 52 human tissues defined in (ICRP Publication 89,, 2002) and used in (ICRP Publication 110,, 2009), according to a monoenergetic X-ray source at 60 keV. As in (Schneider et al.,, 1996), the curves are piecewise linear, with 4 breakpoints at lung tissue, adipose tissue and around soft tissues, although these are some of these are difficult to discern in the figure. It can observed that most tissues have very similar relative electron density and mass density, which deviates more significantly with higher attenuating materials containing bone.

In order to use the calibration curve in Figure 1, the attenuation should be converted to a monoenergetic equivalent, and this is often approximated prior to reconstruction as in (Joseph and Spital,, 1978), which means performing nonlinear calibration both before and after reconstruction, to map into mass or electron density. The alternative that we consider here is to use explicit parameterisations of the attenuation.

2.3 Existing Physical Parameterisations

Due to each cross-section in (4) being a non-linear function of energy and material, the total attenuation of a tissue is complicated and difficult to quantify exactly. One approach is to parameterise (4) as a linear combination of basis functions. For example, in (Alvarez and Macovski,, 1976) a convenient choice is

[TABLE]

where $f_{\mathrm{KN}}(\cdot)$ is the Klein–Nishina function (Klein and Nishina,, 1929) describing the probability of Compton scatter, $Z_{\mathrm{eff}}$ is the effective atomic number for a composite material (Weber and van den Berge,, 1969), and $K_{1}$, $K_{2}$ and $n$ are scalar parameters to fit the model to data. It should be noted that for unbound electrons, one would have $\sigma_{\mathrm{incoh}}=f_{\mathrm{KN}}(\cdot)$. Additionally, the coherent scattering events are not explicitly modelled, though their contribution is low at the energies of interest (Jackson and Hawkes,, 1981).

According to (8), the energy dependent attenuation of any material may be parameterised by its effective atomic number and electron density. Several examples of these parameters for a range of differing material types are plotted in Figure 2(a).

Although there is clearly a strong correlation between the relative electron density $\rho_{e}$ and effective atomic number $Z_{\mathrm{eff}}$, in order to unambiguously model all classes of materials, at least two spectral measurements must be taken to separate the contribution from the two terms in (8)—since there is no one-to-one mapping from $\rho_{e}$ to $Z_{\mathrm{eff}}$. This is the basis of a DECT technique, where from the model in (8) and measurements from two distinct X-ray spectra, the attenuation can be decomposed into photoelectric and Compton, from which $\rho_{e}$ could be unambiguously calculated (Alvarez and Macovski,, 1976).

In (De Man et al.,, 2001), the authors use a piecewise-linear fit from a monoenergetic equivalent attenuation to both of the energy independent factors in (8): $K_{1}\rho_{e}Z_{\mathrm{eff}}^{n}$ for photoelectric attenuation; and $K_{2}\rho_{e}$ for Compton scatter. This fitting allows polyenergetic reconstruction from a single source. Due to the degradation of the energy dependent modelling in (8) at higher effective atomic number or higher energy however (Jackson and Hawkes,, 1981), this will not be quantitatively consistent in $\rho_{e}$ between hard and soft materials; the model is also shown not to be consistent between synthetic and biological materials (De Man et al.,, 2001).

Another idea is to use physical materials as basis functions. For example, for biological specimens, water and bone may be considered (Joseph and Spital,, 1978). The attenuation is then

[TABLE]

where $\rho$ is the mass density, $m_{\mathrm{water}}(\xi)$ and $m_{\mathrm{bone}}(\xi)$ are the energy dependent mass attenuation coefficients, and $a_{1}$ and $a_{2}$ can be binary (Elbakri and Fessler,, 2002) or water–bone fractions (Elbakri and Fessler,, 2003).

One may calculate $\rho$ and map into $\rho_{e}$ through a non-linear calibration if desired. The relation between mass and electron density is shown in Figure 2(b). Although higher density materials deviate from that of water, as long as one can generate an estimate of the mass density, then the trend shown in Figure 2(b) can be approximated as piecewise-linear and conversion to electron density is possible, which is similar to the HU against density plot shown in Figure 1. A possible weakness of this model is the inaccuracies that will occur when tissues have a dissimilar mass-attenuation profile to both water or hard bone, such as adipose tissue (Schneider et al.,, 1996).

3 Methodology

3.1 Polyquant Attenuation Model

We propose to take a data centric rather than a physical approach to parameterising the X-ray attenuation. To motivate our choice for this, we have plotted the relative attenuation coefficient against electron and mass density for 52 biological tissues in (ICRP Publication 89,, 2002) at a number of energies in Figure 3. We have normalised the attenuation to the maximum for an energy—the tooth in each case—simply to allow visualisation on a single graph.

It appears from Figure 3 that the chosen biological tissues approximately follow an energy dependent two part linear fit in both electron and mass density. We have superimposed such a fit onto each plot in Figure 3, where the transition point $k$ is fixed to be the same for each energy and set to minimise the norm of the residual through all the data. Imposing a constant $k$ throughout will motivate computational efficiency outlined in the following section. One interpretation of this parameterisation is taking the inverse of HU to density shown in Figure 1 and incorporating a dependency on energy.

Generalising the piecewise linear fit motivated from the data in Figure 2(b) may be written as

[TABLE]

where $N_{f}$ is the number of linear fits, $f_{i}$ is a class indicator function for materials belonging to the $i^{th}$ class, and $\boldsymbol{x}\in\mathcal{R}^{N_{\mathrm{vox}}}$ is the physical parameter one wishes to fit to—electron density or mass density for example. The class indicator functions are defined as

[TABLE]

where $k_{i}$ is the ‘knee’ between the $i^{th}$ and $(i+1)^{th}$ linear fits, with $k_{0}=0$. To ensure connected fits, we may also enforce the constraints

[TABLE]

We also have $k_{N_{f}}=0$ and $\beta_{1}=0$, which ensures that a vacuum has no attenuation.

Following on from the connection to the decomposed inverse of calibration in Figure 1, this general model would indeed allow several breakpoints around different tissue classes as in (Schneider et al.,, 1996), but here we focus on two or three for computational simplicity.

3.1.1 Special Cases and Connections to Existing Models

We now give several interesting special cases and connections to existing models:

• •

An example of (10) may be used for parameterising the attenuation of biological tissues as a function of relative electron density, which was the singular case studied in (Mason et al.,, 2017), and allows the fitting in Figure 3. For this we have $N_{f}=2$, $x=\rho_{e}$ and $k=1.0064$, which may be written as

[TABLE]

This is the instance we study, along with an equivalent fit for mass density.

• •

The model above may also include highly dense materials, such as titanium metallic implants with a mass density of $\rho=4.5$ g/cm3. To include this, a second knee point may be inserted between bone and metal, and one would have $N_{f}=3$. We also study this in our numerical experiment with hip implants.

• •

Another interesting special case with (10) is when all $b_{i}$ are set to 0, $x=\rho$, $N_{f}=2$, $\alpha_{1}(\xi)$ and $\alpha_{2}(\xi)$ are mass attenuation coefficients for water and bone, and the continuity constraints are ignored. This yields (9), with a density dependent threshold, and is equivalent to the formulation in (Elbakri and Fessler,, 2003) without smoothing. The consequence of ignoring the continuity in (12) will be a step at the ‘knee’ points.

• •

In a similar manner to the water–bone model utilised in (Elbakri and Fessler,, 2003), the IMPACT model in (De Man et al.,, 2001) may also be expressed as a special case of our general fitting model in (10). We note that adopting the same notation, the model may be written as

[TABLE]

where $\alpha_{i}$, $\beta_{i}$, $\gamma_{i}$ and $\delta_{i}$ are the coefficients or photoelectric and Compton linear fits, and $\boldsymbol{x}$ is a monoenergetic equivalent attenuation. A key difference between (14) and (10) is the energy dependence of coefficients and number of parameters. Whilst our model has $2N_{f}N_{\xi}-1$ free parameters, (14) has $4N_{f}$, which is likely to be less, thereby offering increased model flexibility but, as we will see in the next section, at no additional computational cost.

• •

We could also use (10) to predict the reconstruction of a monoenergetic equivalent image—in DECT this is known as ‘quasi monoenergetic’ reconstruction (Johnson,, 2012). This is possible since the imposition of a constant ‘knee’ position across energies allows any monoenergetic attenuation to be predicted using an equivalent piecewise linear function. The relation to the IMPACT model in this case is that the energy dependence is fitted from the data instead of implied from the photoelectric–Compton model in (8).

• •

Although not evaluated in this article, due to the same shape in HU to proton stopping power in (Schneider et al.,, 1996) as relative electron density, we suggest our model should be very applicable for proton interaction modelling also.

3.2 Direct Quantitative Density Reconstruction

The Polyquant model introduced in Section 3.1 describes the forward mapping from physical quantity to attenuation. We now show how this may be combined with the statistical CT measurement model in (3), to allow direct statistical inference of mass or electron density.

Combining (10) with (3) results in the relation

[TABLE]

If we introduce a function $\psi(\cdot,\cdot)$ to simplify notation as

[TABLE]

we can write the negative log-likelihood (NLL) for the Poisson model as

[TABLE]

where we note that this function is non-convex as with similar CT NLL functions in (Erdogan and Fessler, 1999a, ; Chang et al.,, 2014).

Reconstruction of the quantitative density map can be performed by finding an $\boldsymbol{x}$ that minimises (17). We will look at gradient descent methods, for which we require an expression for the derivative of NLL. If we simplify notation with the following

[TABLE]

where $\oslash$ represents component-wise division. An expression for the derivative is then

[TABLE]

where $\boldsymbol{\Phi}^{T}$ represents a transpose of the system matrix or ‘back-projection’, and $\odot$ is component-wise multiplication. We have shown this derivative as an approximation “$\approx$” in testament to the fact that there are discontinuities in the gradient at the ‘knee’ positions, where the gradient is not defined. Although in (Mason et al.,, 2017), we invoked the logistic function to mitigate this effect, we have found better empirical performance by instead using (19) directly.

Before we proceed, we note that (19) only has a single backprojection operation $\boldsymbol{\Phi}^{T}$ per linear fit, so the number of these operations is independent of the number of energies $N_{\xi}$. This is only possible due to constant ‘knee’ positions $k$ for each energy. Calculating $\boldsymbol{\Phi}\boldsymbol{\hat{\mu}}(\boldsymbol{x},\xi)$ is also independent on $N_{\xi}$, and can be evaluated with $2N_{f}-1$ applications of $\boldsymbol{\Phi}$, where one fewer is a consequence of the assumption that $\beta_{1}=0$. For example, applying the method for electron density reconstruction of tissues using (13), one would expect a computation cost of three forward- and two backprojections.

Although one may obtain a quantitative density reconstruction through maximum likelihood estimation, by iteratively minimising (17) through gradient descent with (19), incorporation of prior regularisation will typically improve statistical performance, especially as the noise increases or when few measurements are taken. Incorporating regularisation now gives the penalised log-likelihood or maximum a posteriori estimate as

[TABLE]

where $R(\cdot)$ is some regularisation function, and $\mathcal{C}$ is a set of box constraints on $\boldsymbol{x}$ so that $0\leq x_{i}\leq\zeta\mbox{ for }i=1,...,N_{\mathrm{vox.}}$, where $\zeta$ is the maximum allowable density value, and the constraint set ensures non-negative density values.

The choice of $R(\cdot)$ in (20) will vary on the imaging application, but some possibilities are generalized Gaussian Markov random field (GGMRF) (Bouman and Sauer,, 1993), total variation (TV) (Rudin et al.,, 1992) or wavelet sparsity (Daubechies et al.,, 2008). Without loss of generality, we opt in our experimental section for TV, since it promotes piecewise flat images, which we expect from homogeneous slabs of tissue.

3.2.1 Algorithm Design

Although there are many approaches for evaluating (20), we will detail here the algorithm we have used in our experiments, which we give in Algorithm 1.

We note that Algorithm 1 is an instance of the iPiano (Ochs et al.,, 2014), although we have made some slight changes, where we chose this method due to its analysis with non-convex objective functions as we have in (17). The parameter $\delta$ is the step size, which scales with a factor of $L_{0}$, which we define as

[TABLE]

where $\|\cdot\|_{\infty}$ is the maximum norm of the vector. $L_{0}$ represents the maximum in the diagonal of the Hessian of (19) at the point $\boldsymbol{0}$, and may be considered as a crude estimate of the global Lipschitz constant, $L$. We note that this step size is likely to be very conservative in practice. The parameter $\gamma$ sets the ‘inertia’ of the method, where we used $\gamma=0.8$ as resulted in the fastest performance in (Ochs et al.,, 2014). Finally, $\mathbf{prox}_{\delta R}$ is the proximity operator defined as

[TABLE]

In our experimental section, we use the TV (Rudin et al.,, 1992) as a regularisation function $R$, in which case (22) may be evaluated as in (Beck and Teboulle,, 2009).

To illustrate the convergence properties of our method, and investigate its robustness to more aggressive step sizes, we applied Algorithm 1 to the chest data detailed in Section 4.2, and plotted the evolution of the NLL through iterates in Figure 4, for different step size multiplication factors.

It follows that using $\delta$ is indeed very pessimistic, and one observes convergence in this case even with $128\delta$, where convergence is achieved around 100 iterations. For factors larger than 128, we have found the method fails to converge, and it can be seen than the objective does dot decrease monotonically for more aggressive step sizes. This analysis suggests that the local Lipschitz constant of the gradient through iterations is normally significantly less than the global maximum, or indeed the approximation with $L_{0}$. Adaptive step size schemes could exploit this fact, such as the backtracking line-search in (Ochs et al.,, 2014), but one should be aware that evaluating the cost function in (20) will be more than half the cost of calculating a new gradient, so each backtrack evaluation will increase the iteration cost significantly.

4 Experimentation

For our experiments, we investigate the model accuracy in Section 4.1 followed by a numerical fan-beam CT quantitative reconstruction test in Section 4.2, and finally validation with real data in Section 4.3.

4.1 Polyquant Model Evaluation

To investigate the accuracy of our proposed attenuation model, and its comparison to the other parameterisations given in Section 2.3, we calculated the predicted linear attenuation coefficient for a number of representative tissues from ICRP 89 (ICRP Publication 89,, 2002): adipose tissue, muscle tissue, spongy tissue (upper femoral spongiosa), and hard bone.

For the DECT model in (8), we optimised the three scalar parameters $K_{1}$, $K_{2}$ and $n$ for best fit to all materials in the ICRP 89—shown in Figure 3. We also fitted our piecewise linear Polyquant model in Section 3.1 for both relative electron density and monoenergetic attenuation at 60 keV to the same data. Similarly, we evaluated the accuracy of the IMPACT parameterisation in (De Man et al.,, 2001) to monoenergetic attenuation, again using the same fitting data from ICRP 89. These are plotted along with water and bone attenuation models in Figure 5 and the residual norms are tabulated in Table 1.

An observation that can be made from the plots in Figure 5 and results in Table 1 are that both IMPACT and our Polyquant models are very accurate over all materials. Fitting from monoenergetic attenuation as in IMPACT or ‘Polyquant-tissue-$\mu_{\mathrm{mono}}$’ is understandably more accurate than $\rho_{e}$, due to its closer similarity to other monoenergetic attenuations, and one will inevitably lose this accuracy if one later calibrates to $\rho_{e}$ or $\rho$ using the trend in Figure 1. We note that using our energy dependent fitting strategy in the case of monoenergetic attenuation consistently outperforms IMPACT, and this difference is an order of magnitude for the case of hard bone.

Another feature of the numerical results in Table 1 is that explicitly using the photoelectric–Compton relation in (8) as in DECT is less accurate than fitting it to biological materials as in IMPACT. This highlights that the physical photoelectric–Compton model is not very accurate over a wide range of material types, but also that one should not expect the fitted parameters from IMPACT to necessarily be physically meaningful in terms of $\rho_{e}$ and $Z_{\mathrm{eff}}$.

Finally, the water model is reasonably accurate for the soft tissues and spongy bone and not hard bone, and the bone model conversely so, which suggests the component-wise model in (9) is sensible. However, even if one selected the best cases from bone or water, as is the essence of (Joseph and Spital,, 1978; Elbakri and Fessler,, 2002, 2003), then the error would still be considerably higher than either IMPACT or our proposed approach.

To demonstrate that these fitted parameters are not universal, we also ran the same models with the same fitted parameters on two synthetic plastic materials: Plexiglass®/acrylic (polymethyl methacrylate) and Teflon® (polytetrafluoroethylene). These are interesting materials since they have similar electron densities to muscle and bone respectively but significantly different attenuation. In order to extend these models to include metallic implants, we also looked at how well they may model the attenuation of solid Titanium, where both IMPACT and Polyquant included a second ‘knee’ point to incorporate its attenuation also.

From the plots in Figure 6 and results in Table 1, it is apparent that our proposed model is not universal across material types. Indeed, a water model is the best performing upon the specific plastic case, and the DECT is the best universal model, though this is also to be expected given the data in Figure 2 since it uses a two-dimensional parameterisation.

For the metallic implants, we note how the DECT, water and bone models are very inaccurate. IMPACT also shows significant errors, despite fitting the photoelectric–Compton model directly to Titanium, which provides further evidence that the energy dependence model is not universal. On the other hand, by imposing no physical model, our Polyquant approaches are able to capture the attenuating profile of the implant material very closely.

The implications from the model experiment is that although our model is able to fit very closely to both tissue materials and metallic implants, once these parameters are fitted, they are inaccurate in synthetic plastic materials. An explanation for this mismatch is due to the fact that the correlation between electron density and attenuation coefficient is significantly different with tissues and synthetic materials (Schneider et al.,, 1996). With this, although parameterising the mass attenuation coefficient with water and bone gives more consistent predictions across material types, if one wishes to map into electron density for radiation therapy applications, then one must also take care to fit to the appropriate class of materials.

We finally note that our model is still able to account for plastic materials, and the ‘Polyquant-plastic-$\rho_{e}$’ in Table 2 is a result of fitting to a family of polymers—we used Teflon®, Delrin®, Plexiglass®, polystyrene, LDPE (low-density polyethylene) and PMP (polymethylpentene)—and the resulting accuracy is notably very high. We also show the relationship to attenuation for these materials in the graphs in Figure 7, which confirms that plastics interestingly follow a similar piecewise trend of their own, though this clearly does not extend to Titanium. If one wishes to quantify the attenuation of a mixture of both synthetic and biological tissues, then it seems that no single energy parameterisation would be consistent. One may opt for DECT measurements, but at the cost of requiring two diverse spectral sources.

4.2 Low Dose Numerical Reconstruction Test

The data we used for our numerical reconstruction test were derived from the Adult Reference Computational Phantom (ICRP Publication 110,, 2009), which is a segmented image of defined density and chemical composition to represent real tissues. To investigate the ability to image metal implants, we inserted a pair of prosthetic hip joints (marked in green) with a solid titanium pin and shell. The slices through the chest and pelvis—which are the images we selected for testing—are shown in Figure 8. The resolution is $299\times 137$.

The system geometry used for testing was a flat detector fan-beam CT, which was modelled in the Monte-Carlo package Gate (Jan et al.,, 2011). We included a focused lead strip collimator to mitigate scatter, a 0.6 mm thick CsI scintillator, a shaped source to put more photons into the centre of the specimen, and generated spectra using SpekCalc (Poludniowski et al.,, 2009): 6 mm aluminium filtration at 120 kVp. We simulated a total of $3\times 10^{9}$ photons into 512 detection elements over 360 projections in $1^{\circ}$ increments, which represented a low dose acquisition. Although the spectrum used to generate the measurements was sampled at $1\,\text{keV}$ increments, for reconstruction we used $N_{\xi}=21$, which was deemed a sufficient number in (De Man et al.,, 2001), and corresponded to our sampling of mass attenuation coefficients.

4.2.1 Reconstruction Methods Under Test

Most of the methods we tested are iterative reconstructions, with some degree of polyenergetic modelling. For fairness, we used TV regularisation for each (Rudin et al.,, 1992), and selected the regularisation parameter giving the highest performance in each case — curves showing the influence of these parameters are shown in Figure 10. All iterative methods were run for 500 iterations, which resulted in empirical convergence in each case, after which the figure of merit did not changed by more than the reported precision.

Details of the methods under test are:

• •

FBP: Filtered backprojection (FBP) is a popular reconstruction algorithm, that ‘analytically’ approximates the inverse of linearised CT projections (Fessler,, 2014). In our implementation: we take the logarithm of raw measurements; correct for beam hardening with a water only polynomial fitting (Joseph and Spital,, 1978); apply a ramp filter with Hann windowing, with cut-off frequency optimised to give maximum performance; and finally apply a single backprojection operation $\boldsymbol{\Phi}^{T}$. To map to density from attenuation coefficient, we use the calibration curve in Figure 1 and detailed in Schneider et al., (1996).

We included FBP as a crude baseline and to indicate the level of noise in our system, but we expect it to perform significantly worse than other competitive approaches under test due to the low dose.

• •

PWLS: Penalised weighted-least-squares (PWLS) approximates the CT model as linear by taking the logarithm of the raw measurements (Chang et al.,, 2014), calibrating them to correct for beam-hardening artefacts for the polyenergetic source (Joseph and Spital,, 1978), and includes a statistical weighting to approximate the Poisson noise in (3). The objective function is then

[TABLE]

where $\boldsymbol{W}$ is a diagonal statistical weighting matrix with entries $w_{ii}=(y_{i}-s_{i})^{2}/y_{i}$, and $\boldsymbol{l}$ is the collection of linearised monoenergetic projections (Chang et al.,, 2014). Converting $\boldsymbol{\mu}_{\mathrm{mono}}$ to electron density is then done through a nonlinear calibration according to (Schneider et al.,, 1996), as with the FBP. It should be noted that we are not actively modelling the metal implant in this case.

• •

Poly-SIR: Polyenergetic statistical iterative reconstruction (Poly-SIR) is the segmented water–bone model from (9) (Elbakri and Fessler,, 2002). As this requires prior knowledge on material classes, we give it the oracle segmentation of the hard bone structures, and it treats everything else as water. In the pelvis case, we also pass the oracle segmentation of the metal implants, with corresponding mass attenuation coefficients. Since the physical model gives the mass density, we convert to relative electron density where appropriate using the curve shown in Figure 2(b).

• •

IMPACT: We implement the IMPACT model of (De Man et al.,, 2001) as a special case of our generalised fitting in (14). This allows us to use the same algorithm as presented in Section 3.2.1 for its minimisation. We use a three component piecewise-linear fitting from monoenergetic attenuation at 60 keV to energy dependent attenuation, which accounts for metal implants as well as biological tissues. We then use the same post-processing calibration technique as with FBP and PWLS to convert to either mass or relative electron densities.

• •

Polyquant: In our proposed model, we use the piece-wise linear fitting in Section 3.1 and reconstruction strategy in Section 3.2. We use Algorithm 1, but with $10\delta$ which from Figure 4 is still rather conservative, but still exhibited a monotonic decrease in the objective function. For the pelvis case, we extend the model with a second ‘knee’ point and linear section to include the attenuation of solid titanium according to the generalisation in (10). We used separate fittings to electron and mass density shown in top and middle graphs in Figure 3, to demonstrate the ability to reconstruct directly into either quantity, where each mapping was a least-squares fitting to the materials in ICRP 89 (ICRP Publication 89,, 2002).

4.2.2 Quality Assessment Metric

The metric for quantifying the accuracy of the various methods under test is the root-mean-squared-error (RMSE) in relative electron and mass density, calculated as

[TABLE]

where $\hat{\boldsymbol{\rho}}\in\mathcal{R}^{N_{\mathrm{vox}}}$ is the estimated mass or electron density of a tested method, and $\boldsymbol{\rho}$ is the ground-truth. As an error, the lower the score represents a higher quantitative accuracy.

4.2.3 Reconstruction Results

The results from our reconstruction test are illustrated in Figure 9 and Table 3.

In terms of both electron density and mass density accuracy, our proposed method is the best performing method under test. The lack in estimation accuracy in using the water–bone model with Poly-SIR is likely to be due to discrepancies in adipose and spongy bone tissues from water, which were illustrated in Figure 5. We note that Poly-SIR is second best performing in the pelvis case, but it was provided with oracle information on the implant and hard bone, which would be difficult to segment in practice from a preliminary reconstruction such as FBP due to its high noise. IMPACT is understandably very closely performing to our method in the chest case, as these performed similarly in the model test, but the advantage of fitting at each energy is clear in the case of the metal implant, where the IMPACT performance is similar to that of PWLS.

The relationship between regularisation strength and quantitative performance of the iterative methods are shown in Figure 10. In both cases, this demonstrates that the Polyquant not only reaches the best performance, especially in the pelvis case, but its numerical accuracy is reasonably robust to the setting of this parameter. The difference in the location of the minima between methods may be accounted for by the different scales in image parameter — for example IMPACT uses the monoenergetic attenuation at 60 keV, whereas Polyquant uses the physical density, and these have a relative difference in intensities. Another implication from these curves is that although the Polyquant chest image in Figure 9(i) exhibits a higher level of intratissue variation than IMPACT, this will be mitigated by increasing the regularisation parameter, which can increase by over 40% whilst maintaining the best numerical performance.

Another interesting feature of the results is the difference between electron and mass density scores. Although it may seem counter-intuitive that the Poly-SIR for example would have a higher error before calibration than after, this is due to the RMSE being absolute and not relative. Since from Figures 2(b) and 3 the mass densities of materials are on average higher than relative electron density, this will account for the difference.

4.3 Cone-Beam CT Validation

To validate our method, we also tested its ability to perform quantitative reconstruction from real X-ray cone-beam CT (CBCT) measurements. For this, we acquired a scan with a Varian® TrueBeam™ On-Board Imager® of a CIRS STEEV head phantom. The phantom consists of synthetic resins to mimic the attenuating properties of human tissues, allowing quantitative assessment of relative electron density accuracy. There was also a metal structure in the centre of the phantom, consisting of the plug section from a PTW PinPoint® ionisation chamber, allowing us to investigate the mitigation of metal induced artefacts.

4.3.1 CBCT Data Processing

Our CBCT acquisition consists of 499 projections at a 100 kVp tube potential and 20 mA current for 15 mSec on each, which were the default settings for a head acquisition. Compared to the numerical test, the relative X-ray flux was roughly $2.5\times$ higher, which coupled with the smaller specimen volume and larger number of projections, implies that this test was at a significantly higher dose.

The raw measurements were all pre-corrected for detector responses, and the effect of bow-tie shifting with gantry rotation, with default TrueBeam™ corrections to give $\boldsymbol{y}$ in (3). The scatter estimate $\boldsymbol{s}$ in (3) was also taken from the Varian system’s default scatter correction. For testing the FBP—realised with the Feldkamp–Davis–Kress method (Feldkamp et al.,, 1984)—and PWLS, we calculated the linearised projection vector $\boldsymbol{l}$ in (23).

In the case of the fully polyenergetic reconstruction methods Poly-SIR, IMPACT and Polyquant, we require explicit knowledge of the X-ray spectrum $b_{i}(\xi)$ in (3). Due to the variable thickness of Aluminium in the bow-tie filter, this will be spatially varying, and we calculated it analytically from the appropriate spectrum in the system’s calibration parameters, and the spectral response of the various metal filters and scintillator in the beam path. As in our other experimental sections, we discretised the spectrum into 21 energies, and used the same parameters as were fitted in Section 4.1 for the ICRP biological tissues, but supplemented by the mass attenuation of the metal implant according to information provided by the manufacturer.

4.3.2 CBCT Reconstruction

For reconstruction, we mapped into a resolution of $512\times 512\times 144$, and used each method as detailed in Section 4.2.1. We ran each iterative method for 500 iterations. For the regularisation parameter, we heuristically used $0.5\lambda_{\textrm{pelvis}}$, where $\lambda_{\textrm{pelvis}}$ were the same TV regularisation parameters from the digital pelvis experiment, and gave good empirical performance on the CBCT data. Finally, for the bone and metal segmentations required for Poly-SIR, we obtained these through applying thresholds on the FBP and PWLS separately. To illustrate the critical role of this step, we have shown both images in Figures 11(c) and 11(d).

Reconstructions of the 83rd slice from the CBCT data are shown in Figure 11, along with a region from the 91st slice containing a higher mass of metal in Figure 12. Although the FBP appears to suffer strongly from the presence of the metal structure, most of the iterative methods mitigate its effect considerably, with the PWLS showing a more pronounce dark region in the soft tissue. The Poly-SIR based on a segmentation from the FBP does suffer from its streaking as shown in Figure 11(d). Although a better performance is achieved through segmenting from the PWLS as in Figure 11(d), this will have a considerably higher total computational cost; the Poly-SIR region in Figure 12(c) is based on this PWLS initialisation. The preservation of bone structure appears to be similar between all the iterative methods. The Polyquant method appears to have the smoothest regions around the metal implant, whilst maintaining the same level of structure in the bone and objects surrounding the phantom.

To evaluate the quantitative accuracy of each reconstruction, we isolated regions of interest (ROI) in the 83rd slice located in soft tissue and spongy bone — these are shown in Figure 11(a). The RMSE of these regions are calculated relative to the electron density of the soft tissue and bone equivalent resins, and shown in Table 4.

From Table 4, we note that our proposed model is the most accurate method under test by at least 16% over other approaches in all regions.

5 Discussion

Aspects of our method we have not evaluated in this study are its robustness, and the practicality of its computational implementation, though these are both worth discussing. In the first case, we have noted that due to the discontinuous gradient from the piecewise linear fitting, there is no theoretical guarantees for convergence. Two approaches we have adopted in implementations are using a smooth function for $f(x)$ in (13) such as a generalised logistic, and connecting the two linear fits with a quadratic function for some interval around the ‘knee’. However both options increase the computational cost of the gradient term in (19) considerably, and we have found neither give any empirical advantage in convergence or accuracy over just using the non-smooth version. Since it is common in CT reconstruction to use empirically well performing methods that have no convergence guarantees such as pre-computed curvatures in separable paraboloid surrogates (Erdogan and Fessler, 1999a, ) or ordered subsets (Erdogan and Fessler, 1999b, ), we believe this is reasonable. We also highlight, that from Figure 4, our method does empirically appear to converge even for very aggressive step size multiplication factors.

Another potential robustness issue with any polyenergetic model are partial volume effects, where a discretised voxel contains different classes of material. Due to its linearity, given that the two materials belong to the same fit interval, such as lung and fat or muscle and bone as in Figure 3, then our model would correctly estimate the attenuation from these materials. If a voxel contains materials from different intervals however, such as fat and muscle or metal and bone, then our model will overestimate the attenuation, according to the trends in Figure 3. Although we have not evaluated the degree of this effect, we note it will be in common with other approaches (Elbakri and Fessler,, 2002, 2003; De Man et al.,, 2001), and may be mitigated by increasing the resolution of reconstruction.

The performance of iterative methods under a low dose acquisition is an important consideration for reducing the amount of ionising radiation delivered to the patient. By performing our numerical test at a low dose, we demonstrated the method is relatively robust in this setting. Between the iterative methods tested, PWLS is expected to suffer the most from a low dose, due to its approximation to the noise model, and linearisation of the projections that become unstable for very low photon fluxes (Chang et al.,, 2014), and has been shown to perform worse at low doses than ‘pre-log’ methods such as Polyquant, Poly-SIR and IMPACT (Fu2014). Between these three method, whether there exists a relative change of interplay between partial volume effect photon flux and accuracy is yet to be determined, although Polyquant is likely to benefit from its superior attenuation modelling.

In terms of computational cost, our method is comparable to other full gradient iterative methods. Given that the bottleneck is in calculating the forward- and back-projection operators $\boldsymbol{\Phi}$ and $\boldsymbol{\Phi}^{T}$, then we note that each gradient step in (19) may be implemented with 3 forward and 2 backward evaluations, which represents a $2.5\times$ larger cost than PWLS, which is the same as our implementation of IMPACT. Comparatively, Poly-SIR has a $2\times$ larger cost than PWLS, which given the consistent accuracy advantage of our model is unlikely to be worth this slight speed advantage. For faster implementation, we suggest that algorithmic acceleration such as ordered subsets (Erdogan and Fessler, 1999b, ; Wang et al.,, 2015) and using parallel hardware (Yan et al.,, 2008) are likely suitable, and will be investigated in future work.

6 Conclusions

We have introduced a general quantitative attenuation model, which allows direct inference of mass or electron density from raw CT measurements with a single polyenergetic source. Not only have we demonstrated this allows more accurate modelling than explicit physical models such as water–bone or photoelectric–Compton, but have shown how it may be exploited in a flexible reconstruction algorithm that allows accurate quantitative medical imaging, even with metal implants and real CBCT data. As with other single source methods, we have highlighted its inconsistency between synthetic and biological tissues, but this may not be of relevance for medical imaging, in which opting for the more general DECT model is significantly less accurate over materials of interest. Since our method has a similar computational cost to other iterative approaches, but offers markedly higher accuracy, it offers both a practical and beneficial approach to CT imaging.

Acknowledgements

The authors would like to sincerely thank Dr Adam Wang from Varian Medical Systems, for preprocessing our cone-beam CT data and providing invaluable information and advice about the imaging system. This work was supported by the Maxwell Advanced Technology Fund, EPSRC DTP studentship funds and ERC project: C-SENSE (ERC-ADG-2015-694888). MD is also supported by a Royal Society Wolfson Research Merit Award.