A generalized simultaneous algebraic reconstruction technique (GSART) for dual-energy X-ray computed tomography

Abstract

BACKGROUND:

Dual-energy computed tomography (DECT) is a widely used and actively researched imaging modality that can estimate the physical properties of an object more accurately than single-energy CT (SECT). Recently, iterative reconstruction methods called one-step methods have received attention among various approaches since they can resolve the intermingled limitations of the conventional methods. However, the one-step methods typically have expensive computational costs, and their material decomposition performance is largely affected by the accuracy in the spectral coefficients estimation.

OBJECTIVE:

In this study, we aim to develop an efficient one-step algorithm that can effectively decompose into the basis material maps and is less sensitive to the accuracy of the spectral coefficients.

METHODS:

By use of a new loss function that employs the non-linear forward model and the weighted squared errors, we propose a one-step reconstruction algorithm named generalized simultaneous algebraic reconstruction technique (GSART). The proposed algorithm was compared with the image-domain material decomposition and other existing one-step reconstruction algorithm.

RESULTS:

In both simulation and experimental studies, we demonstrated that the proposed algorithm effectively reduced the beam-hardening artifacts thereby increasing the accuracy in the material decomposition.

CONCLUSIONS:

The proposed one-step reconstruction for material decomposition in dual-energy CT outperformed the image-domain approach and the existing one-step algorithm. We believe that the proposed method is a practically very useful addition to the material-selective image reconstruction field.

Keywords

X-ray imaging computed tomography (CT)Dual-energy CT material decomposition iterative reconstruction one-step method

1 Introduction

Dual-energy computed tomography (DECT) has been widely used in many clinical applications where its advantages over single-energy CT are exploited [1, 2]. DECT can provide material decomposition of the scanned object from the dual-energy data or images. In addition to material decomposition, image artifacts such as beam hardening and scattering can be corrected by DECT for the correction algorithms can effectively accommodate the nonlinear nature of the X-ray attenuation [3, 4]. Recently, along with the advanced particle therapy for cancer patients, the development of imaging techniques that characterize human tissues according to the particle interaction, for which DECT or spectral CT can play major roles, has also been active [5 –8].

Including the first demonstration of dual-energy X-ray imaging by Alvarez and Macovski [9], relatively early methods for material decomposition in DECT can be categorized into image-domain and projection-domain decomposition methods. The image-domain decomposition method performs material decomposition after image reconstruction at each energy [10 –15]. This method is generally applicable not only for the matched ray measurements but also for the unmatched ray measurements, where matching refers to the dual-energy projections at the same source-detector position [16 –19]. However, in the image domain, it is challenging to develop models that can correct for the physical factors such as the etic nature of an X-ray spectrum. On the other hand, the projection-domain decomposition methods perform material decomposition prior to image reconstruction [3 , 20]. Such methods have an advantage of the physical artifact correction over the image-based approach since the nonlinearity of X-ray attenuation can be directly handled in the projection domain. However, angularly well-matched ray sampling is a prerequisite and this scondition often requires additional scanning with the scanned object staying still or expensive technologies such as dual-layer detectors or photon-counting detectors.

In addition to the image-domain and projection-domain decomposition methods, a new category of approaches called ‘joint-reconstruction methods,’ ‘one-step inversion methods,’ or ‘direct reconstruction methods’ have been proposed. In this paper, we refer to them as one-step methods. These methods can overcome several intermingled limitations of the previous approaches by exploiting its non-linear forward model and iterative updating manner [21 –27]. The non-linear forward models, which have been designed to approximate Beer’s law, has a beam-hardening correction effect based on pre-calibrated spectral coefficients, and, unlike the conventional methods that directly perform material decomposition with angularly consistent spectral measurements, the one-step methods largely relax the requirement on the angular consistency by calculating estimations based on the corresponding spectral coefficients and minimizing errors between the estimations and measurements iteratively.

However, the conventional one-step methods have several side effects of using the non-linear forward model. First of all, they have expensive computational costs for convergence. In SECT, many fast-iterative reconstruction algorithms have been developed based on mathematical approximations mostly using convexity of loss functions, but it is not straightforward to demonstrate the convexity of loss functions having the non-linear forward model as like a linear case. Next, in the previous studies, the one-step methods mostly accompany a laborious calibration process. The spectral coefficients are pre-determined by calibration, and the performance of material decomposition largely depended on how accurately the actual spectrum is approximated to the spectral coefficients so an algorithm and a phantom for the calibration should be carefully determined [28]. In addition, for fast convergence of the algorithms, most of the one-step methods conduct initial material decomposition, which requires another calibration and calibration phantom.

In this paper, we propose a new one-step method named generalized simultaneous algebraic reconstruction technique (GSART). The algorithm employs the non-linear forward model and weighted squared errors as a loss function. We demonstrate the convexity of the loss function and derive the algorithm by collaborating with the optimization technique called separable quadratic surrogate (SQS), which was first proposed in SECT and has been successfully used with statistical reconstruction algorithms [29]. The GSART guarantees fast convergence and has the benefits of being less sensitive to noise of individual measurements due to its simultaneous updating scheme. This proposed method is named the GSART because the SART algorithm is one of the particular forms of the GSART [30, 31]. Moreover, the proposed method in this study simplifies the calibration process by performing two calibrations for the spectral coefficients and initial decomposition coefficients with the same spectral measurements of a calibration phantom. It would be laborious and inefficient if two calibrations are separately performed. The calibration phantom is composed of two cylindrical materials made of two representing basis materials.

This paper is organized as follows. In Section II, the derivation of the GSART, the calibration processes, and the details on numerical simulation are described. Reconstructed images and comparison results are presented in Section III and Section IV includes the discussion and conclusion.

2 Methods

2.1 Forward model for DECT

Let $y = (y_{1}, . . ., y_{M}) \in ℝ_{+}^{M}$ be our raw measurements, and $μ (E) = (μ_{1} (E), . . ., μ_{N} (E)) \in ℝ_{+}^{N}$ be the vectorized linear attenuation map in the unit of inverse length, where M and N denote the number of ray measurements and the number of image voxels. We represent the i^th ray measurement as following: ${\hat{y}}_{i} = \int s (E) e^{- {[A μ (E)]}_{i}} dE + r_{i}$ (1)

where s (E) is the spectral intensity of the X-ray as a function of the photon energy E, and A represents the M × N system matrix. r_i is the noise included in the i^th measurement. The noise includes both electrical noise and photon noise.

In DECT, the linear attenuation coefficient is generally divided into two factors, and each factor has both energy dependency and material dependency as follows: $μ (E) = φ (E) f + θ (E) g$ (2)

where f and g represent two material-selective images, or basis material maps. φ (E) and θ (E) are the energy dependencies of the corresponding materials. A functional form of the energy dependency can vary depending on how material-selective images are presented, e.g. density map or attenuation coefficient map [3 , 32]. In this study, we use the linear attenuation coefficient to represent the material-selective images at reference energy E′. In other words, the material-selective images represent monochromatic images at the reference energy. Then the energy dependencies can be written as following: $φ (E) = \frac{μ_{1} (E)}{μ_{1} (E^{'})}, θ (E) = \frac{μ_{2} (E)}{μ_{2} (E^{'})}$ (3)

With a discrete representation of the X-ray spectrum, Equation (1) can be further discretized as following: ${\hat{y}}_{i}^{s} = \sum_{j = 1}^{N_{e}^{s}} w_{j}^{s} e^{- φ_{j} {[Af]}_{i} - θ_{j} {[Ag]}_{i}} + r_{i}$ (4)

where s is the index for spectrum so $w_{j}^{1}$ and $w_{j}^{2}$ represent the discretized spectral intensities for the high-energy and low-energy scans, respectively. Accordingly, a total background signal would be $I^{s} = \sum_{j = 1}^{N_{e}^{s}} w_{j}^{s}$ . $N_{e}^{s}$ denotes the number of energy bins used.

2.2 Optimization algorithm

The loss function of the proposed method is expressed by the sum of the data fidelity term, L (f, g), and the regularization term, R (f, g): $min_{f, g} L (f, g) + λ \cdot R (f, g)$ (5)

where λ is a weighting vector to control the influence between the data fidelity term and the regularization term. In this study, we used the weighted l₂ norm square for the data fidelity term. From Equation (4), when the number of measurements for the dual-energy spectrum is the same, the weighted l₂ norm square for the two spectral logarithmic measurements can be represented by $L (f, g) = \sum_{s = 1}^{2} {(p^{s} - {\hat{p}}^{s})}^{T} W (p^{s} - {\hat{p}}^{s}) = \sum_{s = 1}^{2} {\sum_{i = 1}^{M} W_{i} [p_{i}^{s} + log (\sum_{j = 1}^{N_{e}^{s}} \frac{w_{j}^{s}}{I^{s}} e^{- φ_{j} {[Af]}_{i} - θ_{j} {[Ag]}_{i}})]}^{2},$ (6)

where $p^{s} = (- log (y_{1}^{s} / I^{s}), . . ., - log (y_{M}^{s} / I^{s})) \in ℝ_{+}^{M}$ represents the log-transformed measurement, and ${\hat{p}}^{s}$ represents the estimated attenuations. W is a diagonal matrix in which W_i represents the i^th diagonal element as written by $W_{i} = \frac{1}{\sum_{k = 1}^{N} a_{ik}} .$ (7)

In order to reduce the burden of computing the derivatives of the inseparable loss function with respect to each image voxel, the previous studies have typically used a separable quadratic surrogate function to minimize instead of the given loss function. The above loss function in Equation (6) is inseparable with respect to each voxel of the material-selective images. Thus, we derive our proposed algorithm by applying the SQS to the loss function. The detailed procedure is given in the following section.

2.2.1 Generalized Simultaneous Algebraic Reconstruction Technique (GSART)

To ensure that an algorithm derived from a surrogate function monotonically reduces the original function, the following three conditions should be satisfied [33]: $\begin{matrix} 1 . ψ (ρ; ρ^{(n)}) = Ψ (ρ) \\ 2 . {\nabla_{ρ} ψ (ρ; ρ^{(n)}) |}_{ρ = ρ^{(n)}} = {\nabla Ψ (ρ) |}_{ρ = ρ^{(n)}} \\ 3 . ψ (ρ; ρ^{(n)}) \geq Ψ (ρ), \forall ρ \geq 0 \end{matrix}$ (8) where Ψ (ρ) is an original function for the variable ρ and ψ (ρ ; ρ⁽ⁿ⁾) is its surrogate function determined by the n^th value ρ⁽ⁿ⁾.

We looked for a surrogate function corresponding to the inverse problem which satisfies the above three conditions. Note that we will be showing the derivation of the surrogate function only once for one of the material-selective images, e.g., f_k, since the derivation for the other one, g_k is quite straightforward. First of all, we define loss component as following $l_{i}^{s} ({[Af]}_{i}, {[Ag]}_{i}) = W_{i} {(p_{i}^{s} - {\hat{p}}_{i}^{s})}^{2}$ (9)

for the convenience of derivation. It is convex with respect to f_k and g_k due to the properties of the log-sum-exp function. The proof of the convexity is provided in Appendix A. Now, by using the second-order Taylor series, one can set the following inequality: $\begin{matrix} l_{i}^{s} ({[Af]}_{i}, {[Ag]}_{i}) \leq h_{i}^{s, (n)} ({[Af]}_{i}, {[Ag]}_{i}) \\ = l_{i}^{s} ([A α_{1}]_{i}^{(n)}, [A α_{2}]_{i}^{(n)}) + \sum_{l = 1}^{2} {\frac{\partial l_{i}^{s} ({[A α_{1}]}_{i}, {[A α_{2}]}_{i})}{\partial {[A α_{l}]}_{i}} |}_{α_{1} = α_{1}^{(n)}, α_{2} = α_{2}^{(n)}} ({[A α_{l}]}_{i} - {[A α_{l}^{(n)}]}_{i}) \\ + \frac{1}{2} \sum_{l = 1}^{2} \sum_{m = 1}^{2} C_{l, m, i}^{s, (n)} ({[A α_{l}]}_{i} - {[A α_{l}^{(n)}]}_{i}) ({[A α_{m}]}_{i} - {[A α_{m}^{(n)}]}_{i}), \end{matrix}$ (10)

where α ₁ = f and α ₂ = g are used for the ease of presentation. $C_{i}^{s, (n)}$ is referred to as a curvature. The curvature should be determined so that the above inequality is satisfied. With thus determined curvature, the quadratic surrogate function consequently satisfies the following inequality: $L (f, g) \leq \sum_{s = 1}^{2} \sum_{i = 1}^{M} h_{i}^{s, (n)} ({[Af]}_{i}, {[Ag]}_{i}) = Q_{1} (f, g) .$ (11)

To make it a separable surrogate function with respect to f and g, we apply the convexity method of De Pierro [34]. We rewrite [Af] _i and [Ag] _i as following: $\begin{matrix} {[Af]}_{i} = \sum_{k = 1}^{N} β_{ik} [\frac{a_{ik}}{β_{ik}} (f_{k} - f_{k}^{(n)}) + {[A f^{(n)}]}_{i}] \\ {[Ag]}_{i} = \sum_{k = 1}^{N} β_{ik} [\frac{a_{ik}}{β_{ik}} (g_{k} - g_{k}^{(n)}) + {[A g^{(n)}]}_{i}] \end{matrix},$ (12)

where $β_{ik} = \frac{a_{ik}}{\sum_{k = 1}^{N} a_{ik}}$ and $\sum_{k = 1}^{N} β_{ik} = 1$ . By Jensen’s inequality [35] and the convexity of the quadratic surrogate function, one can place $\sum_{k = 1}^{N} β_{ik}$ out of the function $h_{i}^{s, (n)}$ to set up an inequality relationship below: $Q_{1} (f, g) \leq \sum_{s = 1}^{2} \sum_{i = 1}^{M} \sum_{k = 1}^{N} β_{ik} h_{ik}^{s, (n)} (f_{k}, g_{k}; f^{(n)}, g^{(n)}),$ (13)

where $h_{ik}^{s, (n)} (f_{k}, g_{k}; f^{(n)}, g^{(n)}) = h_{i}^{s, (n)} (\frac{a_{ik}}{β_{ik}} (f_{k} - f_{k}^{(n)}) + {[A f^{(n)}]}_{i}, \frac{a_{ik}}{β_{ik}} (g_{k} - g_{k}^{(n)}) + {[A g^{(n)}]}_{i})$ . Then, the final separable quadratic surrogate function can be expressed as following: $Q_{2} (f, g; f^{(n)}, g^{(n)}) = \sum_{k = 1}^{N} q_{k} (f_{k}, g_{k}; f^{(n)}, g^{(n)}),$ (14)

where $q_{k} (f_{k}, g_{k}; f^{(n)}, g^{(n)}) = \sum_{s = 1}^{2} \sum_{i = 1}^{M} β_{ik} h_{ik}^{s, (n)} (f_{k}, g_{k}; f^{(n)}, g^{(n)})$ .

We are now ready to derive an update algorithm for solving the inverse problem. We first calculate the first and second derivative of q_k (f_k, g_k ; f⁽ⁿ⁾, g⁽ⁿ⁾). Note that the derivatives with respect to f_k of q_k (f_k, g_k ; f⁽ⁿ⁾, g⁽ⁿ⁾) are the same as the derivatives with respect to f_k of Q₂ (f, g ; f⁽ⁿ⁾, g⁽ⁿ⁾) because Q₂ (f, g ; f⁽ⁿ⁾, g⁽ⁿ⁾) is separable. The first and second derivatives with respect to f_k at n^th iteration are expressed as below: $\begin{matrix} {\frac{\partial q_{k} (f_{k}, g_{k}; f^{(n)}, g^{(n)})}{\partial f_{k}} |}_{f_{k} = f_{k}^{(n)}, g_{k} = g_{k}^{(n)}} = \sum_{s = 1}^{2} \sum_{i = 1}^{M} a_{ik} {\dot{h}}_{i}^{s, (n)} ({[A f^{(n)}]}_{i}, {[A g^{(n)}]}_{i}) \\ {\frac{\partial^{2} q_{k} (f_{k}, g_{k}; f^{(n)}, g^{(n)})}{\partial f_{k}^{2}} |}_{f_{k} = f_{k}^{(n)}, g_{k} = g_{k}^{(n)}} = \sum_{s = 1}^{2} \sum_{i = 1}^{M} a_{ik} γ_{i} {\ddot{h}}_{i}^{s, (n)} ({[A f^{(n)}]}_{i}, {[A g^{(n)}]}_{i}), \end{matrix}$ (15)

where $γ_{i} = \sum_{k = 1}^{N} a_{ik}$ . According to Equations (9) and (10), ${\dot{h}}_{i}^{s, (n)} ({[A f^{(n)}]}_{i}, {[A g^{(n)}]}_{i})$ and ${\ddot{h}}_{i}^{s, (n)} ({[A f^{(n)}]}_{i}, {[A g^{(n)}]}_{i})$ are $\begin{matrix} {\dot{h}}_{i}^{s, (n)} ({[A f^{(n)}]}_{i}, {[A g^{(n)}]}_{i}) = {\frac{\partial l_{i}^{s, (n)} ({[Af]}_{i}, {[Ag]}_{i})}{\partial {[Af]}_{i}} |}_{f = f^{(n)}, g = g^{(n)}} \\ {\ddot{h}}_{i}^{s, (n)} ({[A f^{(n)}]}_{i}, {[A g^{(n)}]}_{i}) = C_{1, 1, i}^{s, (n)} \end{matrix},$ (16) and ${\dot{h}}_{i}^{s, (n)} ({[A f^{(n)}]}_{i}, {[A g^{(n)}]}_{i})$ can be further expressed by ${\frac{\partial l_{i}^{s, (n)} ({[Af]}_{i}, {[Ag]}_{i})}{\partial {[Af]}_{i}} |}_{f = f^{(n)}, g = g^{(n)}} = \frac{W_{i} Φ_{i}^{s, (n)}}{{\hat{y}}_{i}^{s, (n)}} ({\hat{p}}_{i}^{s, (n)} - p_{i}^{s}),$ (17)

where $Φ_{i}^{s, (n)} = \sum_{j = 1}^{N_{e}} w_{j}^{s} φ_{j} e^{- φ_{j} {[A f^{(n)}]}_{i} - θ_{j} {[A g^{(n)}]}_{i}}$ . Then, Equation (15) can be rewritten by: $\begin{matrix} {\frac{\partial q_{k} (f_{k}, g_{k}; f_{k}^{(n)}, g_{k}^{(n)})}{\partial f_{k}} |}_{f_{k} = f_{k}^{(n)}, g_{k} = g_{k}^{(n)}} = \sum_{s = 1}^{2} \sum_{i = 1}^{M} a_{ik} \frac{W_{i} Φ_{i}^{s, (n)}}{{\hat{y}}_{i}^{s, (n)}} ({\hat{p}}_{i}^{s, (n)} - p_{i}^{s}) \\ {\frac{\partial^{2} q_{k} (f_{k}, g_{k}; f_{k}^{(n)}, g_{k}^{(n)})}{\partial f_{k}^{2}} |}_{f_{k} = f_{k}^{(n)}, g_{k} = g_{k}^{(n)}} = \sum_{s = 1}^{2} \sum_{i = 1}^{M} a_{ik} γ_{i} C_{1, 1, i}^{s, (n)} . \end{matrix}$ (18)

One of the simple curvatures, $C_{1, 1, i}^{s, (n)}$ , that satisfy the third condition in Equation (8) is the maximum second derivative of $l_{i}^{s, (n)} ({[Af]}_{i}, {[Ag]}_{i})$ in the feasible range. One can obtain the maximum second derivative as follows: $\begin{matrix} C_{1, 1, i}^{s, (n)} = W_{i} {[{(\frac{Φ_{i}^{s, (n)}}{{\hat{y}}_{i}^{s, (n)}})}^{2} + [{(\frac{Φ_{i}^{s, (n)}}{{\hat{y}}_{i}^{s, (n)}})}^{2} - \frac{χ_{i}^{s, (n)}}{{\hat{y}}_{i}^{s, (n)}}] ({\hat{p}}_{i}^{s, (n)} - p_{i}^{s})]}_{+} \\ \leq {\begin{matrix} W_{i} [{(\frac{Φ_{i}^{s, (n)}}{{\hat{y}}_{i}^{s, (n)}})}^{2} + {(\frac{Φ_{i}^{s, (n)}}{{\hat{y}}_{i}^{s, (n)}})}^{2} ({\hat{p}}_{i}^{s, (n)} - p_{i}^{s})], if {\hat{p}}_{i}^{s, (n)} \geq p_{i}^{s} \\ W_{i} [{(\frac{Φ_{i}^{s, (n)}}{{\hat{y}}_{i}^{s, (n)}})}^{2} - \frac{χ_{i}^{s, (n)}}{{\hat{y}}_{i}^{s, (n)}} ({\hat{p}}_{i}^{s, (n)} - p_{i}^{s})], if {\hat{p}}_{i}^{s, (n)} < p_{i}^{s} \end{matrix} \end{matrix},$ (19)

where $χ_{i}^{s, (n)} = \sum_{j = 1}^{N_{e}^{s}} w_{j}^{s} φ_{j}^{2} e^{- φ_{j} {[A f^{(n)}]}_{i} - θ_{j} {[A g^{(n)}]}_{i}}$ . The above inequality holds because all the terms in the equation are non-negative.

2.2.2 Predefined curvature

The maximum second derivative of $l_{i}^{s, (n)} ({[Af]}_{i}, {[Ag]}_{i})$ in Equation (19) is calculated at each iteration, and the process to calculate the second derivative includes an additional back-projection step so it is computationally expensive. Furthermore, the use of maximum second derivative results in a reduced step size of the algorithm and slows down the convergence. One effective way to increase convergence is to employ a predefined curvature. The predefined curvature can be used on a condition that the estimated attenuation, ${\hat{p}}^{s, (n)}$ , does not change substantially as the algorithm goes on. Also, we assumed that initial material decomposed maps derived by the image-domain material decomposition provide the initial estimation, ${\hat{p}}^{s, (0)}$ , that is similar to the real measurement, i.e. ${\hat{p}}^{s, (0)} \approx p^{s}$ . This helps to come up with a predefined curvature that is close to the real curvatures at each iteration. It is noteworthy that the monotonicity of the algorithm is not guaranteed anymore due to the violation of the third condition in Equation (8) when using a predefined curvature. However, since the initial image estimate that is fed into the iterative optimization procedure is reasonably close to the desired solution image in this work, the use of predefined curvature turns out to be empirically acceptable.

Then, we can rewrite the curvature and the second derivative of q_k (f_k, g_k ; f⁽ⁿ⁾, g⁽ⁿ⁾) by applying the assumption as following: $\begin{matrix} C_{1, 1, i}^{s, (n)} \approx W_{i} {(\frac{Φ_{i}^{s, (0)}}{{\hat{y}}_{i}^{s, (0)}})}^{2} \\ {\frac{\partial^{2} q_{k} (f_{k}, g_{k}; f_{k}^{(n)}, g_{k}^{(n)})}{\partial f_{k}^{2}} |}_{f_{k} = f_{k}^{(n)}, g_{k} = g_{k}^{(n)}} = \sum_{s = 1}^{2} \sum_{i = 1}^{M} a_{ik} {(\frac{Φ_{i}^{s, (0)}}{{\hat{y}}_{i}^{s, (0)}})}^{2} \end{matrix},$ (20) since W_iγ_i = 1. The final algorithm for both material-selective images with the predefined curvatures is expressed as following: $\begin{matrix} f_{k}^{(n + 1)} = {[f_{k}^{(n)} + \frac{\sum_{s = 1}^{2} \sum_{i = 1}^{M} a_{ik} \frac{Φ_{i}^{s, (n)}}{{\hat{y}}_{i}^{s, (n)}} \frac{(p_{i}^{s} - {\hat{p}}_{i}^{s, (n)})}{\sum_{l = 1}^{N} a_{il}}}{\sum_{s = 1}^{2} \sum_{i = 1}^{M} a_{ik} {(\frac{Φ_{i}^{s, (0)}}{{\hat{y}}_{i}^{s, (0)}})}^{2}}]}_{+} \\ g_{k}^{(n + 1)} = {[g_{k}^{(n)} + \frac{\sum_{s = 1}^{N_{s}} \sum_{i = 1}^{M} a_{ik} \frac{Θ_{i}^{s, (n)}}{{\hat{y}}_{i}^{s, (n)}} \frac{(p_{i}^{s} - {\hat{p}}_{i}^{s, (n)})}{\sum_{l = 1}^{N} a_{il}}}{\sum_{s = 1}^{N_{s}} \sum_{i = 1}^{M} a_{ik} {(\frac{Θ_{i}^{s, (0)}}{{\hat{y}}_{i}^{s, (0)}})}^{2}}]}_{+} \end{matrix},$ (21)

where $Θ_{i}^{s, (n)} = \sum_{j = 1}^{N_{e}^{s}} w_{j}^{s} θ_{j} e^{- φ_{j} {[A f^{(n)}]}_{i} - θ_{j} {[A g^{(n)}]}_{i}}$ . The classical SART algorithms for single-energy CT reconstruction, which uses a linear system Af = p, is one of the particular forms of this algorithm having parameters as $w_{j}^{s} = {\begin{matrix} I^{s}, if j = 0 \\ 0, otherwise \end{matrix}$ , and φ₀ = 1. With the parameters, the algorithm for the material-selective image, f_k in Equation (21) becomes the classical SART.

2.2.3 Material-specific total-variation (TV) minimizations

An advantage of using the quasi-Newton method in Equation (21) is that one can easily combine the proposed algorithm with any second-order differentiable regularizers. In this study, we adopted the total-variation (TV) loss for regularization to reduce the aforementioned noise [36]. The TV minimization is well known to be effective against angular under-sampling [37]. In DECT, two basis materials having sufficiently different density and attenuation coefficients are typically chosen, and the signal-to-noise ratio (SNR) for two material-selective images can be substantially different. Thus, to enhance the effect of TV minimization, we applied the TV minimization to the two material-selective images with different weights. We empirically determined the weights, λ_f and λ_g, in the experiments. The regularization term used in this study is shown below: $λ \cdot R (f, g) = λ_{f} {∥ f ∥}_{TV} + λ_{g} {∥ g ∥}_{TV}$ (22)

where f and g stand for the acrylic-only image and the aluminum-only image in this paper, respectively. The algorithm with the TV minimization, named GSART-TV, is expressed as following: $\begin{matrix} f_{k}^{(n + 1)} = {[f_{k}^{(n)} + \frac{\sum_{s = 1}^{2} \sum_{i = 1}^{M} a_{ik} \frac{Φ_{i}^{s, (n)}}{{\hat{y}}_{i}^{s, (n)}} \frac{(p_{i}^{s} - {\hat{p}}_{i}^{s, (n)})}{\sum_{l = 1}^{N} a_{il}} + λ_{f} \frac{\partial {∥ f ∥}_{TV}}{\partial f_{k}}}{\sum_{s = 1}^{2} \sum_{i = 1}^{M} a_{ik} {(\frac{Φ_{i}^{s, (0)}}{{\hat{y}}_{i}^{s, (0)}})}^{2} + λ_{f} \frac{\partial^{2} {∥ f ∥}_{TV}}{\partial f_{k}^{2}}}]}_{+} \\ g_{k}^{(n + 1)} = {[g_{k}^{(n)} + \frac{\sum_{s = 1}^{2} \sum_{i = 1}^{M} a_{ik} \frac{Θ_{i}^{s, (n)}}{{\hat{y}}_{i}^{s, (n)}} \frac{(p_{i}^{s} - {\hat{p}}_{i}^{s, (n)})}{\sum_{l = 1}^{N} a_{il}} + λ_{g} \frac{\partial {∥ g ∥}_{TV}}{\partial g_{k}}}{\sum_{s = 1}^{2} \sum_{i = 1}^{M} a_{ik} {(\frac{Θ_{i}^{s, (0)}}{{\hat{y}}_{i}^{s, (0)}})}^{2} + λ_{g} \frac{\partial^{2} {∥ g ∥}_{TV}}{\partial g_{k}^{2}}}]}_{+} \end{matrix}$ (23)

2.3 Calibration

In this study, we start out with the material-decomposed images by use of the image-domain material decomposition (ID-MD) method as initial image estimates of the iterative algorithm. The ID-MD method that expresses the material-selective images as high-order polynomials of the dual-energy images is adopted [3, 38]. It accompanies a calibration process to determine the coefficients of such polynomials, and the calibration process therefore should be performed for each basis material. In this work, we used the calibration phantom which is composed of two separate cylinders made of aluminum and acrylic in both simulation and real experiment. The calibration phantom used in the real experiment is shown in Fig. 1(a). The calibration phantom consisting of two cylindrical phantoms was used for it can provide various thickness combinations of two materials through a circular CT scan with the two cylinders placed in the scanning trajectory plane [3]. We present an example sinogram of the calibration phantom in Fig. 1(b).

The initial material-selective images using the high-order polynomials are expressed as following: $\begin{matrix} f^{(0)} = c_{f} \cdot b (h_{low}, h_{high}) = \sum_{m = 0}^{4} \sum_{n = 0}^{m + n \leq 4} c_{f, m, n} b_{m, n} (h_{low}, h_{high}) \\ g^{(0)} = c_{g} \cdot b (h_{low}, h_{high}) = \sum_{m = 0}^{4} \sum_{n = 0}^{m + n \leq 4} c_{g, m, n} b_{m, n} (h_{low}, h_{high}) \end{matrix}$ (24)

where h_low and h_high represent the low- and high-energy images of the calibration phantom reconstructed by the filtered back projection (FBP) [39]. c_f and c_g are the coefficients for material decomposition. b (h_low, h_high) is a vectorized basis images acquired by use of the high-order polynomials, and its element is represented as $b_{m, n} (h_{low}, h_{high}) = h_{low}^{m} h_{high}^{n}$ with the exponents of m and n. We have empirically determined to use a fourth-order polynomial, meaning m + n ≤ 4, in this work.

In the calibration process, the coefficients, c_f and c_g, are determined by minimizing the error between the high-order polynomials of the dual-energy images and the ground-truth material-selective images, f^* and g^*, respectively. The ground-truth material-selective images are generated by thresholding the reconstructed images of the calibration phantom. In the reconstructed images of the calibration phantom, the region that corresponds to the target material is assigned by the linear attenuation coefficient value at the reference energy, and the rest regions are set to be zero. In this study, we used 50 keV as the reference energy. The minimization problem can be expressed as following: $\begin{matrix} c_{f}^{*} = \underset{c_{f}}{arg min} {∥ M \circ (c_{f} \cdot b (h_{low}, h_{high}) - f^{*}) ∥}_{2}^{2} \\ c_{g}^{*} = \underset{c_{g}}{arg min} {∥ M \circ (c_{g} \cdot b (h_{low}, h_{high}) - g^{*}) ∥}_{2}^{2} \end{matrix}$ (25)

where M is a binary map to exclude unreliable regions, such as edges of the imaged object. Here, ∘ represents the Hadamard product. The solution for each optimization problem is obtained by a well-known least-square solver that uses differentiation of the cost function with respect to each c.

On top of this material decomposition calibration, an additional calibration for estimating the discrete spectral intensities, w^s, for low- and high-energy needs to be performed in the proposed GSART. The discrete spectral intensities in the numerical simulation were assumed to be known. In contrast, for the real experiment, we performed the calibration process to obtain the spectral coefficients for each energy spectrum. For the calibration phantom, the same calibration phantom used in the ID-MD was used. We approximated the whole continuous spectrum to be binned at seven energy levels at 20, 35, 50, 70, 90, 115, and 140 keV. The corresponding energy-dependent coefficients in Equation (2) for aluminum and acrylic, φ (E) and θ (E), were obtained from the National Institute of Standard Technology database (NIST) [40]. The discrete spectral intensities w^*s were determined by minimizing the loss between the measurements and the estimation from each spectrum as following: $w^{* s} = \underset{w^{s}}{arg min} {\sum_{i = 1}^{M} [p_{i}^{s} + log (\sum_{j = 1}^{N_{e}^{s}} \frac{w_{j}^{s}}{I^{s}} e^{- φ_{j} {[A f^{*}]}_{i} - θ_{j} {[A g^{*}]}_{i}})]}^{2}$ (26)

where f^* and g^* are the same ground-truth material-selective images used in Equation (25). The above optimization problem can be minimized straightforwardly with respect to $w_{j}^{s}$ by any gradient-based methods due to the monotonicity of the logarithm. In this work, we used Newton’s method.

2.4 Experiment

We first conducted a numerical phantom simulation. We used a lab-made GPU-based CT forward projection simulator that performs forward projection at each energy bin for an input image having material indices and sums the projections over the energy bins based on the X-ray spectrum input. This simulator also adds Poisson random noise based on the input intensity. The simulated numerical phantom is shown in Fig. 2. It is composed of five pairs of aluminum cylinders and five pairs of low-density acrylic cylinders of different radii inserted in a large acrylic cylinder phantom. We simulated a fan-beam CT system in which the source to detector distance and the source to iso-center distance were 1300 mm and 900 mm. The array size of detectors is 512 × 1, and the pixel pitch is 0.8 mm. For the source, the X-ray spectra at tube voltages of 80 kVp and 140 kVp filtered by a 10 mm aluminum filter were obtained from an open-source X-ray spectrum simulator [41]. We acquired 180 projections for each energy spectrum in an interlacing fashion, similar to the data obtained from the kV-switching system; therefore, the projection angles for each spectral projection are inconsistent. For all spectral ray measurements, we performed simulation by setting the mean number of background counts to 3 × 10⁴ and added Poisson noise corresponding to it.

Additionally, we conducted an experiment with a bench-top cone-beam CT system (EBSCAN 3DCT, EBTECH, Daejeon, Republic of Korea). Two lab-made ellipsoidal phantoms, which are made of acrylic with aluminum cylinder inserts, were used and are shown in Fig. 3. The experimental geometry was set similar to the numerical simulation. The source to detector distance and the source to isocenter distance were set to be 1300 mm and 900 mm, respectively. The detector array size is 1024 × 1024, and the pixel pitch is 0.4 mm. To reduce the scattering effects, a collimator was placed in front of the source so that we acquire fan-beam-like data. For dual-energy spectra, tube voltages of 90 kVp and 140 kVp with 7 mm aluminum filter were used. As were used in the numerical simulation, 180 views for each spectrum were acquired.

In both the simulation and the real experiments, we compared the performance of the proposed method with the ID-MD method using high-order polynomials and the extended algebraic reconstruction technique (E-ART) [28]. The ID-MD method that is based on high-order polynomials is known to less suffer from the beam-hardening artifacts than the linear ID-MD, and therefore has been compared in this study as a representing method in the image-domain approaches. The E-ART is considered a representative one-step iterative image reconstruction approach which updates material-selective images for each ray-by-ray measurement. The E-ART is hard to be combined with a regularizer due to its update nature. We therefore compared it with the non-regularized version of the proposed method (GSART) as well as the regularized version from the data fidelity point of view.

Computations in the proposed method including back-projection and forward-projection were implemented using graphics processing units (GPUs). We used an Intel^® i7-3770 CPU and a GeForce^® GTX 1080 Ti as the computing resources.

3 Results

3.1 Numerical simulation results

We compared the performance of the ID-MD, the E-ART, the GSART, and the GSART-TV. Here, the GSART and GSART-TV refer to the algorithms in Equation (21) and (23). A total number of iterations was set to be 200 for all the iterative algorithms. Reconstructed images of the phantom are shown in Fig. 4. Fig. 4(1) and 4(2) show the material-selective images for acrylic and aluminum, respectively. As seen in the zoomed-in images of Fig. 4(1a), the acrylic-only image of the ID-MD suffers from beam-hardening artifacts and is rather noisy. On the contrary, no obvious beam-hardening artifacts are observable in other iterative algorithm results. The E-ART results show a severely noisy image. The GSART reconstructed image is less noisy than that of the E-ART. In the GSART-TV reconstruction results, the noise is much reduced compared to other images. The weights for TV, λ_f and λ_g, were empirically tuned to suppress noise while preserving the detailed structures: λ_f = 1.0 and λ_g = 4.0. A similar trend was found in the reconstructed aluminum-only images.

For a quantitative evaluation, we calculated root-mean-squared-errors (RMSEs) with respect to the ground truth image. The results consistently support that the GSART-TV performs best in terms of material-selective image reconstruction. The E-ART works noticeably poor particularly in the acrylic-only image, and it is thought to be due to the ray-by-ray updating scheme.

3.2 Experimental results

Figure 5 summarizes the reconstructed material-selective images of the beam-hardening phantom. For the TV weights of the GSART-TV, λ_f = 1.0 and λ_g = 5.0 were chosen. In Fig. 5(1a), the acrylic-only image reconstructed by the ID-MD has severe beam-hardening artifacts particularly in between the aluminum cylinder inserts. The beam-hardening artifacts seem to be rather overly compensated for in the results of the E-ART. Consistently to the numerical study, the GSART-TV results in the best material decomposition performance. The beam-hardening artifacts have not been suppressed as much as in the numerical study in GSART and GSART-TV, which is supposed to be due to the imperfect estimation of the x-ray spectrum in this work. In Table 2, we summarized RMSEs between the ground truth image and the reconstructed material-selective images. The E-ART results show the highest RMSEs than others.

Table 1

RMSEs between the ground truth and the material-selective images of ID-MD, E-ART, GSART, and GSART-TV for the numerical simulation in Fig. 3

	ID-MD	E-ART	GSART	GSART-TV
Acrylic	1.7393	6.2694	1.2713	0.2703
Aluminum	1.5597	0.6994	1.1008	0.2809

Table 2

RMSEs between the ground truth and the material-selective images for the experiment by the ID-MD, E-ART, GSART, and GSART-TV in Fig. 4

	ID-MD	E-ART	GSART	GSART-TV
Acrylic	0.4278	1.3459	0.3703	0.1325
Aluminum	0.4566	1.7396	0.3978	0.1070

We present the reconstructed material-selective image of the dental phantom in Fig. 6. The same TV weights with the previous experiment were used. In the zoomed-in image in Fig. 6(1a), strong beam-hardening artifacts between the aluminum cylinders are observed. In the GSART-reconstructed image, the beam-hardening artifacts are largely reduced but the noise mostly remains. In contrast, both the beam-hardening artifacts and the noise are effectively suppressed in the results of the GSART-TV.

4 Discussion

We developed a one-step algorithm named GSART for material-selective image reconstruction from the dual-energy projection data that are not necessarily view-matched. By exploiting the convexity of the loss function, we have successfully demonstrated the feasibility of the proposed algorithm for practical uses. In both the numerical simulation and the real experiments, the proposed GSART improved the accuracy of material decomposition and reduced the beam-hardening artifacts compared to the ID-MD. The E-ART seems to properly remove the beam-hardening artifacts in the simulation study but not in the real experiments. It is actually overly corrected for the beam hardening artifacts and produces rather aggravated image quality in the real experiment. This is thought to be due to the sensitivity of the E-ART algorithm to the spectral information which is supposed to be accurate in the numerical study but not in the experimental study. However, the GSART is relatively insensitive to the accuracy of the spectral estimation compared to the E-ART. In general, the simultaneous updating method, which is used in the GSART, is known to be less sensitive to data noise than the serial updating method [30]. Inaccuracy in the spectral estimation would lead to increased inconsistency between the real projection data and the forward projection of the image estimate, which may act as a harder challenge in the E-ART case compared to the GSART much like the aforementioned noise sensitivity.

The calibration for estimating the discrete spectral intensities is an additional requirement for using one-step methods. It is challenging in general to acquire an accurate spectrum estimation, and a more sophisticated calibration phantom consisting of multiple material combinations and its associated calibration processes may be accordingly required. Even though a simple phantom consisting of only two basis materials is used in this work, the accuracy of spectrum estimation may still depend on the number of energy bins in the imaging model. In order to check the performance dependence of the GSART algorithm on the number of energy bins, we have investigated RMSE of the reconstructed material-selective images across the number of energy bins. Figure. 7 shows the RMSEs between the ground truth images and the material-selective images as a function of the number of energy bins in the experimental study. The energy bins were decided to be approximately evenly distributed in consideration of the intensity of the X-ray spectrum. For example, we used a narrower bin-width for the lower-energy part of the spectrum and a wider one for the higher-energy part since the X-ray energy spectra typically have higher intensity in the lower-energy part. It is observed that at the number of energy bins from 7 and after the RMSEs are relatively stable and that those become larger at the smaller number of bins.

The one-step algorithm that employs the simultaneous updating was previously introduced by K. Mechlem [21]. However, the algorithm is derived from the Poisson likelihood, whereas the proposed one-step algorithm is derived from the l₂ norm. Compared to the l₂ norm, the Poisson-based likelihood tends to weight less the highly attenuated measurements in the backprojection and addresses data noise accordingly. However, it results in spatial resolution degradation particularly around the high-density materials. Thus, the denoising regularizers such as TV may not contribute to accurate material decomposition in the Poisson-based approaches as much as in the proposed method.

In this study, to focus on the performance of the proposed algorithm, the experiments were conducted in a nearly scatter-free fan-beam-like system. Cooperating with a proper scattering correction method, the proposed method is expected to be also applicable for the wide cone-angle system. Such an application in the cone-beam CT system with a more realistic body phantom will be conducted in our future study.

5 Conclusions

In this work, we have developed an iterative reconstruction method for one-step DECT image reconstruction. The proposed method effectively reduced the beam-hardening artifacts and substantially increased the material decomposition accuracy compared to the ID-MD. Moreover, the calibration processes for the proposed method are considerably simplified by using the same calibration phantom for both initial material decomposition and spectrum estimation. We believe that the proposed method is a unique and practically very useful addition to the existing material-selective image reconstruction techniques.

Footnotes

Acknowledgment

This work was supported in part by the National Research Foundation of Korea under Grant NRF-2019M2A2A4A05031487 and NRF-2020R1A2C2011959, in part by the institute of Civil Military Technology Cooperation funded by the Defense Acquisition Program Administration and Ministry of Trade, Industry and Energy of Korean government under Grant UM19207RD2, in part by the Korea Medical Device Development Fund under Grant 1711137888, and in part by the Ministry of Trade, Industry & Energy(MOTIE, Korea) under Grant 20014921.

Appendix A

This section describes the convexity of Equation (9). For simplicity, the function $l_{i}^{s} ({[Af]}_{i}, {[Ag]}_{i})$ can be summarized as $g (x) = {(f (x) - c)}^{2} (A 1)$

where $\forall x \in ℝ^{N}, \forall c \in ℝ$ and $f (x) = - log (\sum_{i}^{N} e^{- x_{i}})$ . The log-sum-exp function log(∑_ie^{x
_i}) is well known as a convex function and strictly monotonically increasing, so f (x) is concave due to its negative sign. The convexity of the function g (x) is verified if and only if the following inequality condition is satisfied [42]: $\forall x, y \in ℝ^{N}, \forall θ \in [0, 1] : {(f (θ x + (1 - θ) y) - c)}^{2} \leq θ {(f (x) - c)}^{2} + (1 - θ) {(f (y) - c)}^{2} (A 2)$

which is obtained by substituting g (·) for h (·) in the following convex-equivalent condition: $\forall x, y \in ℝ^{N}, \forall θ \in [0, 1] : h (θ x + (1 - θ) y) \leq θ h (x) + (1 - θ) h (y)$

We have the following inequality by arranging the Equation (A2): $f {(θ x + (1 - θ) y)}^{2} - 2 cf (θ x + (1 - θ) y) \leq θ f {(x)}^{2} + (1 - θ) f {(y)}^{2} - 2 c θ f (x) - 2 c (1 - θ) f (y) . (A 3)$

Due to the concavity of f (x), the second term of the left-hand side has the upper bound. It is written as $- 2 cf (θ x + (1 - θ) y) \leq - 2 c θ f (x) - 2 c (1 - θ) f (y) . (A 4)$

Then, we can make the sufficient condition for the inequality condition in Equation (A2) by use of Equations (A3) and (A4): $f {(θ x + (1 - θ) y)}^{2} \leq θ f {(x)}^{2} + (1 - θ) f {(y)}^{2} . (A 5)$

Substituting f (x) = - log(∑_ie^-x_i), it follows: ${(log (\sum_{i} e^{- (θ x_{i} + (1 - θ) y_{i})}))}^{2} \leq θ {(log (\sum_{i} e^{- x_{i}}))}^{2} + (1 - θ) {(log (\sum_{i} e^{- y_{i}}))}^{2} . (A 6)$

From Hölder’s inequality, the left-hand side in Equation (A6) has the upper bound: $\begin{matrix} {(log (\sum_{i} e^{- (θ x_{i} + (1 - θ) y_{i})}))}^{2} \leq {(log ({(\sum_{i} {| e^{- θ x_{i}} |}^{\frac{1}{θ}})}^{θ} \cdot {(\sum_{i} {| e^{- (1 - θ) y_{i}} |}^{\frac{1}{1 - θ}})}^{1 - θ}))}^{2} \\ = {(θ log (\sum_{i} e^{- x_{i}}) + (1 - θ) log (\sum_{i} e^{- y_{i}}))}^{2} \end{matrix} . (A 7)$

We further yield the sufficient condition for Equation (A2) with the right-hand side in Equations (A6) and (A7): ${(θ log (\sum_{i} e^{- x_{i}}) + (1 - θ) log (\sum_{i} e^{- y_{i}}))}^{2} \leq θ {(log (\sum_{i} e^{- x_{i}}))}^{2} + (1 - θ) {(log (\sum_{i} e^{- y_{i}}))}^{2},$

which can be rearranged as $θ (1 - θ) {(log (\sum_{i} e^{- x_{i}}) - log (\sum_{i} e^{- y_{i}}))}^{2} \geq 0 . (A 8)$

The inequality in Equation (A8) is obviously satisfied for $\forall x, y \in ℝ^{N}, \forall θ \in [0, 1]$ so the inequality condition in Equation (A2) is also satisfied. Therefore, the function g (x)in Equation (A1) is convex for $\forall x \in ℝ^{N}$ , $\forall c \in ℝ$ . In our case, we restrict all parameters to be non-negative, i.e. $\forall x, y \in ℝ_{+}^{N}$ and $\forall c \in ℝ_{+}$ , so the convexity is still guaranteed.

References

Landry

, Seco

, Gaudreault

, Verhaegen

, Deriving effective atomic numbers from DECT based on a parameterization of the ratio of high and low linear attenuation coefficients, Phys. Med. Biol. 58(19) (2013), 6851–6866.

Schmidt

, Flohr

, Principles and applications of dual source CT, Phys. Medica 79(October) (2020), 36–46.

Stenner

, Berkus

, Kachelriess

, Empirical dual energy calibration (EDEC) for cone-beam computed tomography, Med. Phys. 34(9) (2007), 3630–3641.

Lee

, Lee

, Min

, et al., Efficient material decomposition method for dual-energy X-ray cargo inspection system, Nucl. Instruments Methods Phys. Res. Sect. A Accel. Spectrometers, Detect. Assoc. Equip. 884(September) (2018), 105–112.

Hünemohr

, Paganetti

, Greilich

, et al., Tissue decomposition from dual energy CT data for MC based dose calculation in particle therapy, Med. Phys. 41(6) (2014), 1–14.

Tremblay

J.É.

, Bedwani

and Bouchard

, A theoretical comparison of tissue parameter extraction methods for dual energycomputed tomography, Med. Phys. 41(8) (2014).

Bär

, Lalonde

, Royle

, et al., The potential of dual-energy CT to reduce proton beam range uncertainties, Med. Phys. 44(6) (2017), 2332–2344.

Men

, Dai

, Chen

, et al., Dual-energy imaging method to improve the image quality and the accuracy of dose calculation for cone-beam computed tomography, Phys. Medica 36 (2017), 110–118.

Alvarez

, Macovski

, Energy selective reconstructions in X-ray Computerized Tomography, Phys. Med. Biol. 21(5) (1976), 733–744.

10.

Maaß

, Baer

and Kachelrieß

, A practical new approach of material decomposition in image domain, Med. Phys. 36(8) (2009), 3818–3829.

11.

Niu

, Dong

, Petrongolo

, Zhu

, Iterative image-domain decomposition for dual-energy CT, Med. Phys. 41(4) (2014), 041901.

12.

Zbijewski

, Gang

G.J.

, Xu

, et al., Dual-energy cone-beam CT with a flat-panel detector: Effect of reconstruction algorithm on material classification, Med. Phys. 41(2) (2014).

13.

Xue

, Ruan

, Hu

, et al., Statistical image-domain multimaterial decomposition for dual-energy CT, Med. Phys. 44(3) (2017), 886–901.

14.

Harms

, Wang

, Petrongolo

, et al., Noise suppression for dual-energy CT via penalized weighted least-square optimization with similarity-based regularization, Med. Phys. 43(5) (2016), 2676–2686.

15.

Tilley

, Zbijewski

, W.

, Stayman, Model-based material decomposition with a penalized nonlinear least-squares CT reconstruction algorithm, Phys. Med. Biol. 64(3) (2019).

16.

Maaß

, Grimmer

and Kachelrieß

, Dual energy CT material decomposition from inconsistent rays (MDIR), IEEE Nucl. Sci. Symp. Conf. Rec. 1(1) (2009), 3446–3452.

17.

Maaß

, Meyer

and Kachelrieß

, Exact dual energy material decomposition from inconsistent rays (MDIR), Med. Phys. 38(2) (2011), 691–700.

18.

Kim

, Ye

J.C.

, W. et al., Sparse-view spectral CT reconstruction using spectral patch-based low-rank penalty, IEEE Trans. Med. Imaging 34(3) (2015), 748–760.

19.

Zhang

, Wang

, Li

, et al., Reconstruction method for DECT with one half-scan plus a second limited-angle scan using prior knowledge of complementary support set (Pri-CSS), Phys. Med. Biol. 65(2) (2020).

20.

Min

, Lee

, Kim

K.W.

, et al., Low-dose dual-energy cone-beam CT using a total-variation minimization algorithm, Med. Imaging 2011 Phys. Med. Imaging 7961(March 2011) (2011), 796132.

21.

Mechlem

, Ehn

, Sellerer

, et al., Joint statistical iterative material image reconstruction for spectral computed tomography using a semi-empirical forward model, IEEE Trans. Med. Imaging 37(1) (2018), 68–80.

22.

Cai

, Rodet

, Legoupil

, Mohammad-Djafari

, A full-spectral Bayesian reconstruction approach based on the material decomposition model applied in dual-energy computed tomography, Med. Phys. 40(11) (2013), 1–16.

23.

Long

, Fessler

, Multi-material decomposition using statistical image reconstruction in X-ray CT, Second Int. Conf. image Form. X-ray Comput. Tomogr. (Mmd) (2012), 413–416.

24.

Weidinger

, Buzug

T.M.

, Flohr

, et al., Polychromatic iterative statistical material image reconstruction for photon-counting computed tomography, , Int. J. Biomed. Imaging 2016 (2016).

25.

Foygel Barber

, Sidky

E.Y.

, Gilat Schmidt

and Pan

, An algorithm for constrained one-step inversion of spectral CT data, Phys. Med. Biol. 61(10) (2016), 3784–3818.

26.

Chen

, Zhang

, Sidky

E.Y.

, et al., Image reconstruction and scan configurations enabled by optimization-based algorithms in multispectral CT, Phys. Med. Biol. 62(22) (2017), 8763–8793.

27.

Mory

, Sixou

, Si-Mohamed

, et al., Comparison of five one-step reconstruction algorithms for spectral CT, Phys. Med. Biol. 63(23) (2018).

28.

Zhao

, Zhao

, Zhang

, An extended algebraic reconstruction technique (E-ART) for dual spectral CT, IEEE Trans. Med. Imaging 34(3) (2015), 761–768.

29.

Erdogan

, Fessler

J.A.

, Ordered subsets algorithms for transmission tomography, Phys. Med. Biol. 44(11) (1999), 2835–2851.

30.

Andersen

, Kak

, Simultaneous algebraic reconstruction technique (SART): A superior implementation of ART, Ultrason. Imaging 6(1) (1984), 81–94.

31.

Jiang

, Wang

, Convergence of the simultaneous algebraic reconstruction technique (SART), IEEE Trans. Image Process. 12(8) (2003), 957–961.

32.

Elbakri

I.A.

, Fessier

J.A.

, Segmentation-free statistical image reconstruction for polyenergetic X-ray computed tomography, Proc. - Int. Symp. Biomed. Imaging 2002-Janua(2) (2002), 828–831.

33.

Meyer

R.R.

, Sufficient conditions for the convergence of monotonic mathematicalprogramming algorithms, J. Comput. Syst. Sci. 12(1) (1976), 108–121.

34.

De Pierro

A.R.

, On the relation between the isra and the em algorithm for positron emission tomography, IEEE Trans. Med. Imaging 12(2) (1993), 328–333.

35.

Jensen

J.L.W.V.

, Sur les fonctions convexes et les inégalités entre les valeurs moyennes, Acta Math. 30(0) (1906), 175–193.

36.

Sidky

E.Y.

, Pan

, Image reconstruction in circular cone-beam computed tomography by constrained, total-variation minimization, Phys. Med. Biol. 53(17) (2008), 4777–4807.

37.

Sidky

E.Y.

, Duchin

, Pan

, Ullberg

, A constrained, total-variation minimization algorithm for low-intensity x-ray CT, Med. Phys. 38(SUPPL.1) (2011), 117–125.

38.

Maaß

, Baer

and Kachelrieß

, Image-based dual energy CT using optimized precorrection functions: A practical new approach of material decomposition in image domain, Med. Phys. 36(8) (2009), 3818–3829.

39.

Feldkamp

L.A.

, Davis

L.C.

, Kress

J.W.

, Practical cone-beam algorithm, J. Opt. Soc. Am. A 1(6) (1984), 612.

40.

Seltzer

, Tables of X-ray mass attenuation coefficients and mass energy-absorption coefficients, NIST standard reference database 126, NIST Standard Reference Database 126 (1996).

41.

P. Desprás, “SpectrumGUI,” 2013. [Online]. Available: https://sourceforge.net/projects/spectrumgui/

42.

Boyd

, Vandenberghe

, Convex Optimization. Cambridge University Press (2004).