A Splitting-Based Iterative Algorithm for GPU-Accelerated Statistical Dual-Energy X-Ray CT Reconstruction

TL;DR

This paper introduces a GPU-accelerated, splitting-based iterative algorithm for statistical dual-energy CT reconstruction, significantly improving speed and quality in baggage material classification tasks.

Contribution

It presents a novel splitting-based ADMM algorithm that separates reconstruction and decomposition, enabling faster convergence in dual-energy CT imaging.

Findings

Significant reduction in reconstruction time demonstrated on synthetic and real baggage data.

Improved image quality over traditional methods.

Effective handling of high dynamic range in material properties.

Abstract

When dealing with material classification in baggage at airports, Dual-Energy Computed Tomography (DECT) allows characterization of any given material with coefficients based on two attenuative effects: Compton scattering and photoelectric absorption. However, straightforward projection-domain decomposition methods for this characterization often yield poor reconstructions due to the high dynamic range of material properties encountered in an actual luggage scan. Hence, for better reconstruction quality under a timing constraint, we propose a splitting-based, GPU-accelerated, statistical DECT reconstruction algorithm. Compared to prior art, our main contribution lies in the significant acceleration made possible by separating reconstruction and decomposition within an ADMM framework. Experimental results, on both synthetic and real-world baggage phantoms, demonstrate a significant reduction in time required for convergence.

Figures

Tables

Table 1. Table 1 : Minimum number of operations per ADMM iteration. In LM( m 𝑚 m ) × \times CG( n 𝑛 n ) , n 𝑛 n CG iterations are taken to compute the update per LM iteration. Note that LM may require more operations than listed to tune the damping parameter.

	$𝑹$	${𝑹}^T$	$f(\cdot )$
LM($m$)$\times$CG($n$)	$2m\left(n+1\right)$	$2m\left(n+1\right)$	$4m$
UDM($m$)-PCG($n$)	$2n$	$2\left(n+1\right)$	$4m$

Equations

$$\mu (E,x)=x_c{f}_{\text{KN}}(E)+x_p{f}_p(E),$$

$${\mathbf{m}}_{h/l}=-\ln \int S_{h/l}(E)e^{-\mathbf{R}\mathbf{\mu }(E)}dE+\ln \int S_{h/l}(E)dE,$$

$$f_{h/l}(\tilde{\mathbf{R}}\mathbf{x})=C_{h/l}-\ln \int S_{h/l}{e}^{-f_{\text{KN}}\mathbf{R}{\mathbf{x}}_c-f_p\mathbf{R}{\mathbf{x}}_p}dE$$

$${\hat{\mathbf{x}}}_{\text{MAP}}=\underset{\mathbf{x}}{\text{mbox}argmin}\frac{1}{2}{\Vert f(\tilde{\mathbf{R}}\mathbf{x})-\mathbf{m}\Vert }_{\Sigma }^2+\lambda |\mathbf{x}{|}_{\text{TV}}+g(\mathbf{x}),$$

$${\hat{\mathbf{x}}}_{\text{MAP}}=$$

$$\text{s.t.}\left[\begin{array}{c}{\mathbf{y}}_{c/p}\\ {\mathbf{z}}_{c/p}\end{array}\right]=\left[\begin{array}{c}\mathbf{D}\\ \mathbf{I}\end{array}\right]{\mathbf{x}}_{c/p}=\mathbf{C}{\mathbf{x}}_{c/p},$$

$$L(\mathbf{x},\mathbf{y},\mathbf{z})=$$

$$+\sum _{\beta \in \{c,p\}}\frac{{\rho }_{\beta }}{2}{\Vert \left[\begin{array}{c}{\mathbf{y}}_{\beta }\\ {\mathbf{z}}_{\beta }\end{array}\right]-\mathbf{C}{\mathbf{x}}_{\beta }+{\mathbf{u}}_{\beta }\Vert }^2,$$

$${\hat{\mathbf{x}}}_c=$$

$$+\frac{{\rho }_c}{2}{\Vert \left[\begin{array}{c}{\mathbf{y}}_c\\ {\mathbf{z}}_c\end{array}\right]-\mathbf{C}{\mathbf{x}}_c+{\mathbf{u}}_c\Vert }^2,$$

$${\hat{\mathbf{x}}}_p=$$

$$+\frac{{\rho }_p}{2}{\Vert \left[\begin{array}{c}{\mathbf{y}}_p\\ {\mathbf{z}}_p\end{array}\right]-\mathbf{C}{\mathbf{x}}_p+{\mathbf{u}}_p\Vert }^2.$$

$${\hat{\mathbf{y}}}_{c/p}$$

$${\hat{\mathbf{z}}}_{c/p}$$

$${\hat{\mathbf{u}}}_{c/p}={\mathbf{u}}_{c/p}+\left[\begin{array}{c}{\hat{\mathbf{y}}}_{c/p}\\ {\hat{\mathbf{z}}}_{c/p}\end{array}\right]-\mathbf{C}{\hat{\mathbf{x}}}_{c/p}.$$

$${\hat{\mathbf{x}}}_{c/p}=\underset{{\mathbf{x}}_{c/p}}{\text{mbox}argmin}\frac{{\rho }_{c/p}}{2}{\Vert \left[\begin{array}{c}{\mathbf{a}}_{c/p}\\ {\mathbf{y}}_{c/p}\\ {\mathbf{z}}_{c/p}\end{array}\right]-\left[\begin{array}{c}\mathbf{R}\\ \mathbf{C}\end{array}\right]{\mathbf{x}}_{c/p}+\left[\begin{array}{c}{\mathbf{u}}_{c/p}^a\\ {\mathbf{u}}_{c/p}^y\\ {\mathbf{u}}_{c/p}^z\end{array}\right]\Vert }^2.$$

$$\left[\begin{array}{c}{\hat{\mathbf{a}}}_c\\ {\hat{\mathbf{a}}}_p\end{array}\right]=$$

$$+\sum _{\beta \in \{c,p\}}\frac{{\rho }_{\beta }}{2}{\Vert {\mathbf{a}}_{\beta }-\mathbf{R}{\hat{\mathbf{x}}}_{\beta }+{\mathbf{u}}_{\beta }^a\Vert }^2$$

$$\xi (\mathbf{x})=20{\log }_{10}\left(\frac{{\Vert \mathbf{x}-{\mathbf{x}}^{\ast }\Vert }_2}{{\Vert {\mathbf{x}}^{\ast }\Vert }_2}\right).$$

Full text

**Index Terms— ** DECT, ADMM, GPU

1 Introduction

X-ray Computed Tomography (CT) is a widely deployed technique for threat detection in baggage at airport security checkpoints. However, using a single X-ray spectrum limits its reconstruction to only that of LAC (linear attenuation coefficient) or HU (Hounsfield Unit) images which are, at best, an approximation to the underlying energy-dependent characteristics of the materials. While LAC might suffice for discriminating between different types of tissues in medical imaging, in the adversarial case of airport baggage scanning, materials used commonly for homemade explosives can possess LAC values that are nearly the same as for benign materials. As a result, the technique of DECT, where projections from two X-ray spectra are collected simultaneously at each angle, has been proposed for material discrimination because it allows for the material property to be recovered in an additional dimension by using an energy-dependent attenuation model.

The most commonly used DECT model represents X-ray attenuation as the combined effect of Compton scattering and photoelectric absorption, which can be written as:

[TABLE]

where $f_{\text{KN}}(E)$ and $f_{p}(E)$ denote the energy-dependent multipliers111See [1] for a detailed formulation of $f_{\text{KN}}(E)$ and $f_{p}(E)$. to the Compton coefficient $x_{c}$ and the photoelectric (PE) coefficient $x_{p}$, respectively. Therefore, the goal of DECT is to recover both the Compton and the PE coefficients of the object using the projections $\boldsymbol{m}_{h},\boldsymbol{m}_{l}\in\mathbb{R}^{M}$ measured at two different X-ray spectra. This amounts to solving simultaneously for $\boldsymbol{x}_{c},\boldsymbol{x}_{p}\in\mathbb{R}^{N}$ from the following two equations that can be shown compactly by:

[TABLE]

where $\boldsymbol{R}$ denotes the forward projection matrix, and $S_{h}(E)$ and $S_{l}(E)$ denote, respectively, the photon distribution across energy levels in the high- or low-energy spectra.

Solving the two equations shown above is made challenging by the fact the two phenomena that contribute to X-ray attenuation — Compton and PE, occur at grossly different scales. At most applicable energy levels, Compton scattering is the dominant contributor to attenuation and this disparity between the two phenomena becomes worse with increasing energy since $f_{p}(E)$ has a cubic decay with $E$. As a result, stable recovery of PE coefficients requires more sophisticated inversion algorithms, such as those described in [2, 3, 4, 5].

Unfortunately, the algorithms cited above tend to be iterative and, with run-of-the-mill computing hardware, take a long time to return the results. Therefore, a straightforward implementation of these algorithms is not appropriate for the end-goals that motivate our research — high-throughput baggage screening at airports. The focus of the concepts presented in this paper is on improving the computational efficiency of an ADMM-based statistical inversion algorithm.

2 Related Work

DECT involves two key steps: dual-energy decomposition and tomographic reconstruction. The two tasks can be done either sequentially, as in projection-wise decomposition, or in a unified step using iterative statistical approaches.

Projection-wise Decomposition: One of the earliest approaches for DECT decomposition, the Constrained Decomposition Method (CDM) [1] involves directly decomposing the dual energy projections to Compton and PE line integrals followed by the FBP reconstructions of the two. In [6], CDM was also extended to operate for Multi-Energy CT. A major disadvantage of this method is that it guards itself poorly against artifacts, especially in PE coefficients but it is still a preferred approach as it enables parallel implementation.

Iterative Statistical Approaches: The statistical methods for DECT solve for the MAP estimates, finding the Compton and PE coefficients that best correspond to the measurements and any prior knowledge. The literature on Multi-Energy CT has focused largely on designing the models and the priors that best leverage the structural similarity across bases. In [4], Compton/PE images are reconstructed on a set of exponential basis functions with an edge-correlation penalty term. This idea of encouraging geometric similarity between the more stable Compton image and the PE image is further explored in [5]. For this purpose, the authors have proposed a new Non-Local Mean (NLM) regularizer on the PE image and laid out an ADMM formation that scaled up to a problem size of practical interest. By treating the energy level as a third dimension, the contributions in [2, 3] adopt a tensor-based model for reconstruction using sparse priors.

Despite being able to produce high-quality DECT, statistical reconstructions generally fail drastically with regard to the timing constraints of practical applications. Compared to LAC reconstructions, DECT has to deal with the added computational burden associated with the decomposition step. As a result, in existing approaches, solving decomposition and reconstruction in one step becomes inefficient due to the combined complexity and high-dimensionality of the problem.

We therefore propose a statistical DECT approach that combines the best of the projection-wise decomposition methods and the iterative statistical methods. More specifically, we employ a splitting-based MAP estimation method embedded in an ADMM framework. As we will show in Section 4, our new splitting scheme not only provides a better convergence rate, but also allows for powerful hardware-enabled acceleration.

3 Proposed Method

3.1 Problem Formulation

To describe our proposed method, we first define the forward model for X-ray attenuation, i.e., the nonlinear transform $f(\cdot)$ from Compton/PE coefficients to logarithmic projections:

[TABLE]

where $C_{h/l}$ are constants, $\boldsymbol{x}=\begin{bmatrix}\boldsymbol{x}^{T}_{c};\boldsymbol{x}^{T}_{p}\end{bmatrix}^{T}$ and $\boldsymbol{\widetilde{R}}=\begin{bmatrix}\boldsymbol{R};\boldsymbol{R}\end{bmatrix}$. Given a pair of line-integral measurements $\boldsymbol{m}=\begin{bmatrix}\boldsymbol{m}^{T}_{h};\boldsymbol{m}^{T}_{l}\end{bmatrix}^{T}$ corrupted by Poisson photon noise and possibly other artifacts, we construct a MAP estimate with the Total Variation (TV) term and the non-negativity prior:

[TABLE]

where $\Sigma$ is a diagonal matrix with photon counts as trace elements as proposed in [7], $g(x_{i})=\left\{\begin{matrix}0,x_{i}\in\mathbb{R}^{+}\\ \infty,x_{i}\notin\mathbb{R}^{+}\end{matrix}\right.$ and $\lambda$ is the regularization parameter. While the unconstrained optimization in (4) is highly inefficient to solve directly, the Alternating Direction Method of Multipliers (ADMM) method provides a flexible splitting-based framework for both optimality and convergence [5].

3.2 Formulation for ADMM

ADMM begins with the conversion of (4) to its constrained equivalent. By introducing two new auxiliary variables $\boldsymbol{y}$ and $\boldsymbol{z}$ for the TV and the non-negativity terms, respectively, and by posing $\boldsymbol{x}$ as the primal variable, the MAP minimization can be transformed into:

[TABLE]

where $\boldsymbol{D}$ denotes the finite difference operator. Subsequently, the corresponding Augmented Lagrangian (AL) can be written as:

[TABLE]

where $\boldsymbol{u}$ denotes the dual variable, or the scaled Lagrangian multiplier, and $\rho$ denotes the penalty parameter. ADMM splits the minimization of AL into subproblems that are solved separately in an iterative framework. One such intuitive splitting scheme proposed in [5] is for the iterative updates to be carried out according to (7) through (11) shown below. First, the primal variables are updated using:

[TABLE]

[TABLE]

Note that the update to $\boldsymbol{x}_{p}$ is made subsequent to $\boldsymbol{\widehat{x}}_{c}$ since we can expect $\boldsymbol{\widehat{x}}_{c}$ to stabilize the recovery of the more noise-prone $\boldsymbol{x}_{p}$. Secondly, we update the auxiliary variables $\boldsymbol{y}$ and $\boldsymbol{z}$ corresponding to the TV and non-negativity terms, respectively, by using

[TABLE]

where the shrinkage function is defined in [8, Eq. 9.35]. Then the dual variable $\boldsymbol{u}$ is updated with

[TABLE]

Lastly, we update $\rho$ adaptively using a method described in [9] that is based on the primal and dual residual. While (9)-(11) can be realized by straightforward element-wise operations, solving (7)-(8) requires a nonlinear least squares algorithm such as Levenberg-Marquardt (LM). Despite its robustness, LM, is a sequential algorithm, offering little room for parallelization. Intuitively, the overall computational inefficiency can be attributed to (7) or (8) which at once addresses two problems of very different nature: nonlinear dual-energy decomposition and regular linear tomographic reconstruction.

3.3 Proposed Splitting Scheme

For the new sped-up implementation presented in this paper, we further decompose (7)-(8) into two simpler subproblems: tomographic reconstruction followed by dual-energy decomposition. This is achieved by viewing measurement fitting as an additional constraint: $\boldsymbol{a}=\boldsymbol{\widetilde{R}}\boldsymbol{x}$, where $\boldsymbol{a}$ is a new auxiliary variable for that purpose. Incorporating this constraint with the others, for the first subproblem we pose tomographic reconstruction as an unweighted unconstrained least squares problem:

[TABLE]

and we can now separately deal with the decomposition, i.e. solving for $\boldsymbol{a}$, in the following subproblem. We have chosen to use a CG solver for the minimization shown above because our system is huge and sparse. Additionally, since the projection-wise weights $\Sigma$ are now absorbed by the decomposition step to be described later, the linear system corresponding to (12) is shift-invariant. Therefore, for improved convergence rate, each CG iteration lends itself well to preconditioning using a high-passing filter [10]. In our implementation, we experimented with the ramp filter as the preconditioner.

The second subproblem, that of dual-energy decomposition, can be solved by finding the pair of Compton/PE coefficients that minimizes the cost function at each ray projection independently. Therefore, the decomposition is achieved by using the Unconstrained Decomposition Method (UDM):

[TABLE]

Compared to (7) and (8), the update formulas in (12) and (13) have the following advantages: First, treating reconstruction and decomposition separately allows us to use two entirely different algorithms, with each tailored for the corresponding subproblem, which in our case are Preconditioned CG (PCG) and UDM. As an additional benefit, while LM is still used in UDM, backtracking of the damping parameter, which is necessary for dealing with nonlinear optimization, no longer involves expensive tomographic operations. Lastly, computational efficiency is improved in terms of both the number of operations required, as listed in Table 1, and parallelization. We can now not only independently execute (12) for both bases, but also leverage massively parallelized hardware such as GPUs to solve (13).

4 Experimental Results

We now present qualitative and quantitative results on three phantoms: a simulated phantom (sim18), a real-world water bottle phantom (Water), and an actual baggage phantom (Clutter). We have implemented the following algorithms and evaluated them on these three phantoms: (i) CDM-FBP: CDM [1] with 16 CPU threads; (ii) LM($m$)$\times$CG($n$): ADMM as in [5]; and (iii) UDM($m$)-(P)CG($n$): proposed ADMM method with $m$ UDM iterations and $n$ (P)CG iterations.

Our Python implementations use the Astra toolbox [11] for GPU-based calculations of forward and backward projections. Additionally, as a key factor for acceleration, we use Gpufit [12] to parallelize the decomposition step. All experiments are carried out on a single computing cluster node with 16 cores, 20GB RAM and an Nvidia 1080Ti GPU. Regarding initialization, in our experiment we have chosen the CDM output as the initial Compton estimate, since it is generally stable. For initial PE coefficients, we use a scaled version of the Compton estimate, similar to what is done in [5]. For both Compton and PE, we used $\lambda=10^{-5}$ as the TV regularization parameter and $\rho=10^{-3}$ as the initial penalty parameter. Furthermore, for scale-invariant evaluation, quantitative reconstruction quality is assessed with the normalized $\ell_{2}$-distance between $\boldsymbol{x}$ and ground truth $\boldsymbol{x}^{*}$:

[TABLE]

The first phantom, sim18, contains seven circular regions of different materials. The reconstructed images, of size $512\times 512$, are constructed from 720 angles, each containg 725 parallel ray projections. The first row in Figure 1 shows the reconstructions by CDM-FBP and the proposed UDM-PCG algorithm. In Figure 2a, we display the average computational time per iteration inside parentheses in the inset box and plot $\xi(\boldsymbol{x})$ versus iteration to compare the convergence rates. In general, the ADMM algorithm with the proposed splitting scheme results in significantly shorter total execution time than LM$\times$CG to reach the same error.

For results on the two real-world phantoms, the parallel beam projection data for the two was collected on the Imatron C300 CT scanner. The X-ray source emits $1.8\times 10^{5}$ and $1.7\times 10^{5}$ photons per ray with two energy spectra at 95keV and 130keV, respectively. The high- and the low-energy sinograms are subsampled by two, resulting in 360 angles for each and with 512 bins for each angle. The reconstruction results for the Water phantom are shown in the second row of Figure 1. For quantitative evaluation, in Figure 2b we plot $\xi(\boldsymbol{x})$ versus iteration within the ROI – the central circular region that is occupied by distilled water. Our proposed UDM-PCG algorithm is significantly more efficient not only in terms of the number of the iterations needed, but also in terms of the time per iteration as compared to LM$\times$CG.

The real-world baggage phantom Clutter is particularly challenging because it contains metallic objects, and illustrates well the need for a statistical reconstruction. As shown in the third row of Figure 1, while the CDM-FBP PE reconstruction is completely overshadowed by streaking artifacts, UDM1-PCG1 is able to recover object shapes to a reasonable degree.

5 Conclusions and Future Work

In this paper, we have proposed a new splitting scheme for implementing ADMM to reconstruct Compton and PE coefficient images using dual-energy projection data. By separating the the reconstruction and decomposition steps, the proposed GPU-accelerated ADMM algorithm achieves a significant speedup in time when compared to the prior state-of-the-art. Future work will be aimed at generalization to 3D reconstruction and improving initialization heuristics.