Image-domain multimaterial decomposition for dual-energy computed tomography with nonconvex sparsity regularization

Qihui Lyu; Daniel O'Connor; Tianye Niu; Ke Sheng

doi:10.1117/1.JMI.6.4.044004

15 October 2019 Image-domain multimaterial decomposition for dual-energy computed tomography with nonconvex sparsity regularization

Qihui Lyu, Daniel O'Connor, Tianye Niu, Ke Sheng

Author Affiliations +

Journal of Medical Imaging, Vol. 6, Issue 4, 044004 (October 2019). https://doi.org/10.1117/1.JMI.6.4.044004

Abstract

Dual-energy computed tomography (CT) has the potential to decompose tissues into different materials. However, the classic direct inversion (DI) method for multimaterial decomposition (MMD) cannot accurately separate more than two basis materials due to the ill-posed problem and amplified image noise. We propose an integrated MMD method that addresses the piecewise smoothness and intrinsic sparsity property of the decomposition image. The proposed MMD was formulated as an optimization problem including a quadratic data fidelity term, an isotropic total variation term that encourages image smoothness, and a nonconvex penalty function that promotes decomposition image sparseness. The mass and volume conservation rule was formulated as the probability simplex constraint. An accelerated primal-dual splitting approach with line search was applied to solve the optimization problem. The proposed method with different penalty functions was compared against DI on a digital phantom, a Catphan^® 600 phantom, a quantitative imaging phantom, and a pelvis patient. The proposed framework distinctly separated the CT image up to 12 basis materials plus air with high decomposition accuracy. The cross talks between two different materials are substantially reduced, as shown by the decreased nondiagonal elements of the normalized cross correlation (NCC) matrix. The mean square error of the measured electron densities was reduced by 72.6%. Across all datasets, the proposed method improved the average volume fraction accuracy from 61.2% to 99.9% and increased the diagonality of the NCC matrix from 0.73 to 0.96. Compared with DI, the proposed MMD framework improved decomposition accuracy and material separation.

1. Introduction

Conventional x-ray computed tomography (CT) projections are acquired with a single-energy spectrum. The reconstructed single-energy CT (SECT) provides the linear attenuation coefficients (LACs) of the scanning object, but the LACs depend on both the effective atomic number and the electron density, making it insufficient to determine the material components. Dual-energy CT (DECT) is acquired with two distinct energy spectra that are attenuated differently by the tissues. Therefore, DECT provides enhanced information to better differentiate and quantify material compositions. DECT shows promise in many clinical applications, including virtual unenhancement imaging,¹^,² liver lesion characterization,³ kidney stone characterization,⁴ oncologic imaging,¹^,⁵ bone removal,¹ etc.

Despite the potential of using DECT for multimaterial decomposition (MMD), the problem is ill-defined for decomposing more than two materials without additional assumption.⁶^,⁷ Existing DECT MMD approaches can be classified into three categories: projection-based method, integrated method, and image-based method. The projection-based method decomposes the two independent sinograms into two basis components by interpolating the lookup table and then performs separate image reconstruction.⁸ The integrated method incorporates the DECT acquisition into the forward projection model and then reconstructs basis material images directly from the DECT projections.⁹ The projection-based method and integrated method are robust to beam hardening artifacts,⁸^,⁹ but there are several limitations. First, the projection-based method is limited to decomposing two basis materials. It is neither straightforward nor desirable to introduce additional constraints in the projection domain. Second, the projection-based method requires strict spatial and temporal consistency between the high-energy (HE) and low-energy (LE) acquisition. The condition is not met by many DECT systems using dual-source or single-source with fast kilovoltage-switching. Third, the integrated method is computationally expensive due to the repeated projection and backprojection during the reconstruction of basis image components. Fourth, the coupling between the reconstruction and image decomposition in the integrated method, especially with multiple basis materials, leads to decomposition results that are sensitive to parameter tuning.

Alternatively, image-based decomposition methods were investigated, which perform basis material decomposition on the HE and LE CT images. This method is applicable to all DECT systems, straightforward, and computationally tractable. Mendonça et al.¹⁰ proposed an image-based MMD method that introduces two assumptions: mass and volume conservation, and that there are no more than three materials in each voxel. For each voxel, the method loops over a material triplet library, identifies the best-fit triplet, and decomposes the voxel into the three materials via direct matrix inversion. The method is termed direct inversion (DI), and it demonstrates the capability of decomposing CT into more than three basis materials. However, the decomposition accuracy is sensitive to the selection of material triplet library, which limits the flexibility of MMD, and the resultant images suffer from substantially amplified noise.

In this study, we propose a DECT MMD method that decomposes the CT volume into multiple basis materials accurately while suppressing the decomposition image noise. Instead of rigidly confining each voxel to contain at most three basis materials, a nonconvex sparsity term is introduced to penalize the number of materials that are simultaneously present in the same voxel.

2. Method

2.1.

Formulation

The image domain DECT MMD problem is formulated as follows:

Eq. (1)

\underset{x}{minimize} \frac{1}{2} {‖ A x - b ‖}_{2}^{2} + λ TV (x) + η \sum_{k, i} {| x_{k i} |}^{α} subject to x \in S_{P},

where the optimization variable

x

is the volume fraction (VF) matrix of dimension

N_{k}

by

N_{i}

. The voxels on the component image are indexed by

i

, and the decomposition materials are indexed by

k

.

x_{k i}

is the VF of the

k

’th material for the

i

’th voxel.

A

is the component-to-attenuation transformation matrix defined as

A = [\begin{matrix} μ_{1 l} & μ_{2 l} & \dots & μ_{N_{k} l} \\ μ_{1 h} & μ_{2 h} & \dots & μ_{N_{k} h} \end{matrix}],

where

μ_{k l}

and

μ_{k h}

are the LACs of the

k

’th material at HE and LE, respectively. Matrix

A

transforms the VF into the LACs at HE and LE for each voxel.

b

is the attenuation image at HE and LE from measurements, defined as

b = [\begin{matrix} b_{1 l} & b_{2 l} & \dots & b_{N_{i} l} \\ b_{1 h} & b_{2 h} & \dots & b_{N_{i} h} \end{matrix}],

where

b_{i l}

and

b_{i h}

denotes the LAC of the

i

’th voxel on the CT image at LE and HE, respectively.

TV (x)

is the isotropic total variation regularization, defined as

TV (x) = \sum_{i k} \sqrt{{(\nabla_{r} x)}_{i k}^{2} + {(\nabla_{c} x)}_{i k}^{2}},

where

\nabla_{r}

and

\nabla_{c}

are the derivative matrix along the row and column directions, respectively.

S_{P}

is the probability simplex set that enforces the mass and volume conservation rule on each voxel, such that each column of the VF matrix

x

satisfies the sum-to-one and non-negative constraints:

Eq. (2)

x \in S_{P} \leftrightarrow {\begin{cases} x_{k i} \geq 0, & \forall i, k \\ \sum_{k} x_{k i} = 1, & \forall i \end{cases} .

The first term is the quadratic data fidelity term that minimizes the difference between the measured DECT image and the estimated DECT image calculated from the VF matrix $x$ . The second term is the isotropic TV regularization term, applied on each component image to encourage image smoothness while preserving image edges. The last term is the penalty function in the form of ${| \cdot |}^{α}$ ( $0 \leq α \leq 1$ ), which promotes sparsity on the number of materials that simultaneous present in the same voxel. In this study, we specifically focus on four cases with $α = 0, \frac{1}{2}, \frac{2}{3}, 1$ . In the case of $α = 1$ , the sparsity term reduces to a constant under the mass and volume conservation constraint $(x \in S_{P})$ . In other cases, the penalty functions are nonconvex with the nonconvexity increases as $α$ goes to 0.

2.2.

Algorithm

This study utilizes an accelerated primal-dual splitting approach with line search for both convex and nonconvex problems,¹¹^,¹² which solves the optimization problem of the form:

Eq. (3)

minimize F (K x) + G (x),

where

G

is convex,

F

possibly nonconvex, and

K

a linear operator.

The optimization problem in Eq. (1) is formulated into the canonical form shown in Eq. (3) by defining

K = [\begin{matrix} A \\ D \\ I \end{matrix}], G (x) = I_{S_{p}} (x) = {\begin{cases} \infty, & x \notin S_{p} \\ 0, & x \in S_{p} \end{cases},

F ([\begin{matrix} {\hat{z}}_{1} \\ {\hat{z}}_{2} \\ {\hat{z}}_{3} \end{matrix}]) = F_{1} ({\hat{z}}_{1}) + F_{2} ({\hat{z}}_{2}) + F_{3} ({\hat{z}}_{3}),

F_{1} ({\hat{z}}_{1}) = \frac{1}{2} {‖ {\hat{z}}_{1} - b ‖}_{2}^{2},

F_{2} ({\hat{z}}_{2}) = λ {‖ {\hat{z}}_{2} ‖}_{1},

{(F_{3} ({\hat{z}}_{3}))}_{i} = η {| {\hat{z}}_{3 i} |}^{α} .

F (K x)

is equivalent to the objective function in Eq. (1).

G (x)

is an indicator function that equals to 0 if

x \in S_{p}

and infinity if

x \notin S_{p}

, which enforces the sum-to-one and non-negative constraints.

The accelerated primal-dual algorithm with line search¹² is presented in Algorithm 1, where the key steps are the evaluations of the proximal operator of functions $G$ and $F$ . The proximal operator of a function $H$ with step size $t$ is defined as follows:¹³

Eq. (4)

{Prox}_{t H} (z) = \underset{x}{\arg \min} (H (x) + \frac{1}{2 t} {‖ x - z ‖}_{2}^{2}) .

Following the definition, the proximal operator of

G

reduces to the projection onto the probability simplex

S_{p}

, which can be efficiently solved by sorting and thresholding the input vector, as presented in Algorithm 2.¹⁴

Algorithm 1

Accelerated primal-dual algorithm with line search.

Initialization:

x_{0} ≔ 0 \in X

,

z_{0} ≔ 0 \in Z

,

t_{0} > 0

,

θ_{0} ≔ 1

,

β_{0} ≔ 1

,

γ > 0

,

r ≔ 0.8

For

k = 1,2, \dots

, do

x_{k} ≔ {Prox}_{t_{k - 1}, G} (x_{k - 1} - t_{k - 1} K^{T} z_{k - 1})

β_{k} ≔ β_{k - 1} (1 + {γ t}_{k - 1})

t ≔ t_{k - 1} \sqrt{\frac{β_{k - 1}}{β_{k}} (1 + θ_{k - 1})}

Repeat

θ_{k} ≔ \frac{t}{t_{k - 1}}

{\bar{x}}_{k} ≔ x_{k} + θ_{k} (x_{k} - x_{k - 1})

s_{k} ≔ β_{k} t

{\bar{z}}_{k} ≔ z_{k - 1} + s_{k} K {\bar{x}}_{k}

z_{k} ≔ {\bar{z}}_{k} - s_{k} {Prox}_{s_{k}^{- 1}, F} (\frac{{\bar{z}}_{k}}{s_{k}})

t_{k} ≔ t

Break if

\sqrt{β_{k}} t {‖ K^{T} (z_{k} - z_{k - 1}) ‖}_{2} \leq {‖ z_{k} - z_{k - 1} ‖}_{2}

t ≔ t_{k} * r

End

Algorithm 2

Proximal operator evaluation of G (ProxtG(x)).

Input:

x \in R^{N_{k}} \times R^{N_{i}}

For

i = 1,2, \dots, N_{i}

do

Sort

{x_{1 i}, x_{2 i}, \dots, x_{k i}, \dots, x_{N_{k} i}}

into

x^{'}

such that

{x^{'}}_{1 i} \geq {x^{'}}_{2 i} \geq \dots \geq {x^{'}}_{N_{k} i}

Find

J_{i} ≔ \max {1 \leq k \leq K | {x^{'}}_{k i} + \frac{1}{k} (1 - \sum_{j = 1}^{k} {x^{'}}_{j i})}

Define

{\hat{x}}_{i} ≔ \frac{1}{J} (1 - \sum_{j = 1}^{J} {x^{'}}_{j i})

For

k = 1, 2, \dots, N_{k}

do

{({Prox}_{t G} (x))}_{i k} ≔ \max {x_{i k} + {\hat{x}}_{i}, 0}

Output:

{Prox}_{t G} (x)

With the separable sum rule, the evaluation of the proximal operator of $F$ reduces to evaluating the proximal operator with respect to each variable:

{Prox}_{w F} ([\begin{matrix} {\hat{z}}_{1} \\ {\hat{z}}_{2} \\ {\hat{z}}_{3} \end{matrix}]) = ([\begin{matrix} {Prox}_{w F_{1}} ({\hat{z}}_{1}) \\ {Prox}_{w F_{2}} ({\hat{z}}_{2}) \\ {Prox}_{w F_{3}} ({\hat{z}}_{3}) \end{matrix}]) .

Following the definition of proximal operator in Eq. (4), the proximal operators of

F_{1}

and

F_{2}

are

{Prox}_{w F_{1}} ({\hat{z}}_{1}) = \frac{{\hat{z}}_{1} + w b}{w + 1},

{({Prox}_{w F_{2}} ({\hat{z}}_{2}))}_{i k} = {\begin{cases} (1 - \frac{w λ}{{‖ {({\hat{z}}_{2})}_{i k} ‖}_{2}}) {({\hat{z}}_{2})}_{i k}, & {‖ {({\hat{z}}_{2})}_{i k} ‖}_{2} > w λ \\ 0, & {‖ {({\hat{z}}_{2})}_{i k} ‖}_{2} \leq w λ \end{cases} .

Due to the separability, the proximal operator of

F_{3}

reduces to pointwise proximal operator evaluation of

η {| \cdot |}^{α}

, as follows:

Eq. (5)

{({Prox}_{w F_{3}} ({\hat{z}}_{3}))}_{i} = \underset{x}{argmin} (η {| x |}^{α} + \frac{1}{2 w} {(x - {\hat{z}}_{3 i})}^{2}) .

Exact analytic solutions to Eq. (5) exist for the scenarios considered in this study

(α = 0, \frac{1}{2}, \frac{2}{3}, 1)

, which can be found in Sec. 6.1. For other values of

α

, the proximal operators could be evaluated numerically with an iterative approach, such as Newton’s method.¹⁵

With the evaluations of the proximal operator of functions $G$ and $F$ , the DECT MMD optimization problem in Eq. (1) could be solved following Algorithm 1.

2.3.

Evaluation

The proposed framework with different sparsity parameters in the penalty functions was evaluated on a digital phantom, a Catphan^® 600 phantom, a quantitative imaging phantom, and a pelvis patient, and is compared with the DI method.¹⁰ The quantitative imaging phantom data were acquired on a Siemens SOMATOM Force DECT. The Catphan and pelvis patient data were acquired on a Siemens SOMATOM Definition Flash. For the pelvis patient data, the mAs and CT dose index for the [LE, HE] acquisitions were [170 mAs, 131 mAs] and [4.92 mGy, 3.66 mGy], respectively. The LE and HE CT images, as shown in Figs. 1, 3, 5, and 7, were reconstructed using the standard filtered backprojection (FBP). The HE LAC and LE LAC of the basis material were computed as the average values of the HE LACs and LE LACs of the region of interest (ROIs), as shown on the LE CT images for all cases. These LAC values are presented in Table 4 in Sec. 6.

Fig. 1

(a) The LE: 75 kVp and (b) the HE: 140 kVp CT image of the digital phantom. The components of the ROIs are bone (ROI1), iodine (ROI2), water (ROI3), and air (ROI4). The displaying window is $[0.01, 0.065] {mm}^{- 1}$ .

For quantitative evaluation of the material decomposition accuracy, the mean and standard deviation (STD), electron density, as well as the VF, were computed within each ROI on the decomposition component image. The electron density $ρ_{i}$ at voxel $i$ is calculated by

ρ_{i} = \sum_{k = 1}^{N_{k}} x_{k i} ρ (k),

where

ρ (k)

is the electron density of the

k

’th material. The VF accuracy in a uniform ROI is defined as follows:

VF = (1 - \frac{{‖ \bar{x} - x_{g} ‖}_{2}}{{‖ x_{g} ‖}_{2}}) \times 100 %,

where

\bar{x}

is the mean material component vector over all voxels within the uniform ROI and

x_{g}

is the material component vector of the ground truth decomposition.

To quantify the amount of overlap between different material decompositions across the whole image, the normalized cross correlation (NCC) coefficients at zero lag are evaluated for every pair of materials. The NCC coefficient at zero lag of material $k_{1}$ and $k_{2}$ is defined as follows:

Eq. (6)

R_{k_{1} k_{2}} = \frac{\sum_{i = 1}^{N_{i}} x_{k_{1} i} x_{k_{2} i}}{\sqrt{\sum_{i = 1}^{N_{i}} x_{k_{1} i}^{2}} \sqrt{\sum_{i = 1}^{N_{i}} x_{k_{2} i}^{2}}},

which equals to 0 if the two materials are completely separated and equals to 1 if they have identical distribution on the whole image. The NCC matrix is defined as a matrix with entries of

R_{k_{1} k_{2}}

. The diagonal elements of the NCC matrix are equal to 1 by definition, and the other elements are between 0 and 1, indicating the extent of material separation. If all materials are completely separated for all voxels, then all the off-diagonal elements should be 0. The diagonality

D

of the NCC matrix is computed using the Pearson correlation coefficient:

D = \frac{(\sum_{k_{1}, k_{2}} R_{k_{1} k_{2}}) (\sum_{k_{1}, k_{2}} k_{1} k_{2} R_{k_{1} k_{2}}) - (\sum_{k_{1}, k_{2}} k_{1} R_{k_{1} k_{2}}) (\sum_{k_{1}, k_{2}} k_{2} R_{k_{1} k_{2}})}{\sqrt{(\sum_{k_{1}, k_{2}} R_{k_{1} k_{2}}) (\sum_{k_{1}, k_{2}} k_{1}^{2} R_{k_{1} k_{2}}) - {(\sum_{k_{1}, k_{2}} k_{1} R_{k_{1} k_{2}})}^{2}} \sqrt{(\sum_{k_{1}, k_{2}} R_{k_{1} k_{2}}) (\sum_{k_{1}, k_{2}} k_{2}^{2} R_{k_{1} k_{2}}) - {(\sum_{k_{1}, k_{2}} k_{2} R_{k_{1} k_{2}})}^{2}}},

where

D

equals to 1 for the diagonal matrix,

- 1

for the antidiagonal matrix, and 0 for the uniform matrix. In the case of every voxel on the CT image being composed of only a single material, the diagonality of the NCC matrix is 1. On the other hand, in the case where every voxel contains an equal amount of all basis materials, the diagonality of the NCC matrix is 0. The diagonality

D

of the NCC matrix summarizes the decomposition separability over the whole image and across all materials.

3. Results

3.1.

Digital Phantom

Figure 1 shows the LE and HE CT image of the digital phantom. This simple digital phantom is made up of four basis materials, including bone, iodine, water, and air, corresponding to the four ROIs indicated by the rectangular area. Figure 2 shows the decomposition results on the digital phantom using the proposed framework with different $α$ values and the classic DI method. In the cases when $α = 0, \frac{1}{2}, \frac{2}{3}$ , the proposed framework distinctly separated all basis materials and with low noise, achieving a clear border between different materials. On the contrary, neither DI nor the proposed framework with $α = 1$ achieved the desired sparsity or assigned the correct material to each ROI. The material separation capability of the proposed framework with $α = 0, \frac{1}{2}, \frac{2}{3}$ is further confirmed with the NCC map in Fig. 2, showing little or no cross talks between two materials. With DI, the iodine component is mixed up with the bone and water. In the case of $α = 1$ , none of the component is separated.

Fig. 2

Decomposition component images of (1) bone, (2) iodine, (3) water, and (4) air, decomposed using the proposed framework when (a) $α = 0$ , (b) $α = \frac{1}{2}$ , (c) $α = \frac{2}{3}$ , (d) $α = 1$ , and (e) the DI method. The last column is the $4 \times 4$ NCC map of the decomposition in the same row, with each square showing the corresponding entries in the NCC matrix, where the basis materials are bone, iodine, water, and air, from top to bottom and from left to right.

3.2.

Catphan^® 600 Phantom

Figure 3 shows the LE and HE CT images of the Catphan^® 600 phantom with contrast rods, which are made of six basis materials, including Teflon, Delrin, iodine solution of $10 mg / ml$ , polymethyl pentene (PMP), inner soft tissue, and air, corresponding to the six ROIs on the LE CT image [Fig. 3(a)]. The materials in the labeled rods on the HE CT [Fig. 3(b)] are Teflon, Delrin, iodine solution of $10 mg / ml$ , polystyrene, low-density polyethylene, PMP, and iodine solution of $5 mg / ml$ , respectively. The VF accuracy was evaluated on the six ROIs, and the electron density was evaluated on the seven contrast rods. Figure 4 shows the decomposition image and the NCC map. The proposed framework with $α = 0, \frac{1}{2}, \frac{2}{3}$ successfully separated the phantom into the six basis components with minimal cross talk between different materials. With the classical DI method, the off-diagonal elements of the NCC matrix are up to 0.36, showing that the two corresponding basis materials, Delrin and PMP, are not well-separated. The proposed method with $α = 1$ is unable to achieve material separation at all. Table 1 shows the evaluated electron densities for the seven contrast rods. The mean square error of the electron density was reduced by 72.6% for the proposed framework with $α = 0,1 / 2, 2 / 3$ compared with DI.

Fig. 3

(a) The LE: 75 kVp and (b) the HE: 125 kVp CT image of the Catphan^® 600 phantom. The material components of the ROIs on the low energy CT image are Teflon (ROI1), Delrin (ROI2), iodine of $10 mg / ml$ (ROI3), PMP (ROI4), inner soft tissue (ROI5), and air (ROI6). The labeled contrast rods on the HE CT image are composed of (1) Teflon, (2) Delrin, (3) iodine solution of $10 mg / ml$ , (4) polystyrene, (5) low-density polyethylene, (6) PMP, and (7) iodine solution of $5 mg / ml$ . The displaying window is $[0.01, 0.04] {mm}^{- 1}$ .

Fig. 4

Decomposition component images of (1) Teflon, (2) Delrin, (3) iodine solution of $10 mg / ml$ , (4) PMP, (5) inner soft tissue (LDPE), (6) air, decomposed using the proposed framework when (a) $α = 0$ , (b) $α = \frac{1}{2}$ , (c) $α = \frac{2}{3}$ , and (d) $α = 1$ ; (e) the DI method. The last column is the $6 \times 6$ NCC map of the decomposition in the same row, with each square showing the corresponding entries in the NCC matrix, where the basis materials are Teflon, Delrin, iodine solution of $10 mg / ml$ , PMP, LDPE, and air, from top to bottom and from left to right.

Table 1

The electron densities measured on the Catphan contrast rods labeled in Fig. 3(b). The RMSE is evaluated for each method as the mean square error of the seven rods.

Method		Rod 1 Teflon	Rod 2 Delrin	Rod 3 Iodine (10 mg/ml)	Rod 4 polystyrene	Rod 5 LDPE	Rod 6 PMP	Rod 7 Iodine (5 mg/ml)	Average NMSE
Ground truth	$ρ_{e}$	6.240	4.525	3.368	3.400	3.155	2.851	3.356	\
$α = 0$	$ρ_{e}$	6.240	4.525	3.368	3.104	3.008	2.851	3.354	1.92%
$α = 0$	NMSE	0.00%	0.00%	0.00%	8.70%	4.67%	0.00%	0.06%	1.92%
$α = 1 / 2$	$ρ_{e}$	6.218	4.514	3.357	3.119	2.997	2.854	3.352	2.06%
$α = 1 / 2$	NMSE	0.35%	0.24%	0.33%	8.26%	5.01%	0.12%	0.13%	2.06%
$α = 2 / 3$	$ρ_{e}$	6.211	4.507	3.346	3.113	2.984	2.856	3.346	2.27%
$α = 2 / 3$	NMSE	0.46%	0.41%	0.66%	8.43%	5.43%	0.16%	0.30%	2.27%
$α = 1$	$ρ_{e}$	6.169	4.203	4.011	2.961	2.770	2.509	3.578	10.16%
$α = 1$	NMSE	1.14%	7.11%	19.11%	12.91%	12.22%	12.01%	6.63%	10.16%
DI	$ρ_{e}$	4.847	4.174	3.699	3.251	3.134	2.887	3.523	7.32%
DI	NMSE	22.33%	7.76%	9.83%	4.38%	0.67%	1.28%	4.96%	7.32%

3.3.

Quantitative Imaging Phantom

Figure 5 show the DECT image for the quantitative imaging phantom, which consists of 12 basis materials, including $2 mg / ml$ iodine solution (ROI1), $5 mg / ml$ iodine solution (ROI2), $10 mg / ml$ iodine solution (ROI3), $15 mg / ml$ iodine solution (PMP) (ROI4), $50 mg / ml$ calcium solution (ROI5), $100 mg / ml$ calcium solution (ROI6), $300 mg / ml$ calcium solution (ROI7), HE blood 70 (ROI8), HE blood 100 (ROI9), adipose (ROI10), water (ROI11), and brain (ROI12). In this phantom study, to avoid being trapped in undesired local minima with the increased number of basis materials, the algorithm is initialized by setting $VF = 1$ for the basis material that is closest to each pixel. Figure 6 shows the decomposition images with the corresponding NCC map. The proposed framework with $α = 0, \frac{1}{2}, \frac{2}{3}$ decomposed the DECT images into 12 different basis components plus air and improved the VF accuracy from 51% using DI method to 100%. The nondiagonal entries of the NCC map are close to 0 for the proposed framework with $α = 0, \frac{1}{2}, \frac{2}{3}$ , showing clean separation of the basis materials. The off-diagonal entries are close to 1 for the $α = 1$ scenario. The classical DI method is unable to separate similar materials. For example, the HE blood 100 (ROI9), adipose (ROI10), and water (ROI11) have similar LAC values, resulting in substantial nonzero NCC elements $R_{k_{1} k_{2}}$ using DI.

Fig. 5

(a) The LE: 100 kVp and (b) the HE: 150 kVp CT image of the quantitative imaging phantom. The material components of the ROIs are $2 mg / ml$ iodine solution (ROI1), $5 mg / ml$ iodine solution (ROI2), $10 mg / ml$ iodine solution (ROI3), $15 mg / ml$ iodine solution (PMP) (ROI4), $50 mg / ml$ calcium solution (ROI5), $100 mg / ml$ calcium solution (ROI6), $300 mg / ml$ calcium solution (ROI7), HE blood 70 (ROI8), HE blood 100 (ROI9), adipose (ROI10), water (ROI11), brain (ROI12), and air (ROI13).

Fig. 6

Decomposition component images of (1) $2 mg / ml$ iodine solution, (2) $5 mg / ml$ iodine solution (ROI2), (3) $10 mg / ml$ iodine solution (ROI3), (4) $15 mg / ml$ iodine solution (PMP) (ROI4), (5) $50 mg / ml$ calcium solution (ROI5), (6) $100 mg / ml$ calcium solution (ROI6), (7) $300 mg / ml$ calcium solution (ROI7), (8) HE blood 70 (ROI8), (9) HE blood 100 (ROI9), (10) adipose (ROI10), (11) water (ROI11), (12) brain (ROI12), and (13) air (ROI13), decomposed using the proposed framework when (a) $α = 0$ , (b) $α = \frac{1}{2}$ , (c) $α = \frac{2}{3}$ , and (d) $α = 1$ , and (e) the DI method. The last column is the NCC map of the corresponding decomposition, with each square showing one NCC matrix element.

3.4.

Pelvis Patient

Figures 7 and 8 show the DECT image and the decomposition images with the corresponding NCC map, respectively, for the pelvis patient. The proposed framework with $α = 0, \frac{1}{2}, \frac{2}{3}$ decomposed the DECT images into bone (ROI1), iodine (ROI2), muscle (ROI3), fat (ROI4), and air, and achieved an NCC matrix with nondiagonal coefficients close to 0.

Fig. 7

(a) The LE: 100 kVp and (b) the HE: 140 kVp CT image of the pelvis patient. The material components of the ROIs are bone (ROI1), iodine (ROI2), muscle (ROI3), fat (ROI4), and air (ROI5). The displaying window is $[0.01, 0.04] {mm}^{- 1}$ .

Fig. 8

Decomposition component images of (1) bone (2) iodine (3) muscle (4) fat and (5) air, decomposed using the proposed framework when (a) $α = 0$ , (b) $α = \frac{1}{2}$ , (c) $α = \frac{2}{3}$ , and (d) $α = 1$ ; (e) the DI method. The last column is the $5 \times 5$ NCC map of the decomposition in the same row, with each square showing the corresponding entries in the NCC matrix, where the basis materials are bone, iodine, muscle, fat, and air, from top to bottom and from left to right.

Table 2 summarizes the VF accuracy and diagonality of the NCC matrix for all datasets. The proposed method with $α = 0, \frac{1}{2}, \frac{2}{3}$ achieves a significantly higher VF accuracy and diagonality than the comparison methods and approaching nearly perfect material decomposition. Across all datasets, the proposed method improved the average VF accuracy from 61.2% to 99.9% and increased the diagonality of the NCC matrix from 0.73 to 0.96.

Table 2

The VF accuracy and diagonality of the NCC matrix for all datasets.

		Digital phantom	Catphan	Quantitative imaging phantom	Pelvis patient
VF accuracy	$α = 0$	100.0%	99.9%	100%	NA
	$α = 1 / 2$	100.0%	99.6%	100%
	$α = 2 / 3$	100.0%	99.3%	100%
	$α = 1$	62.1%	45.3%	15%
	DI	88.4%	44.1%	51%
Diagonality	$α = 0$	1.00	0.95	1.00	0.95
	$α = 1 / 2$	0.99	0.94	0.99	0.89
	$α = 2 / 3$	0.99	0.93	0.99	0.87
	$α = 1$	0.33	0.14	0.03	0.53
	DI	0.87	0.72	0.54	0.8

Table 3 presents the runtime and the hyperparameters used in all cases. Despite that our algorithm used MATLAB built-in GPU computing tools for acceleration, the proposed framework takes 12 min on average, which is slower than the DI method, requiring only 3 min on CPU. The hyperparameters were tuned case-by-case to achieve visually desired sparseness and smoothness. The sparseness and the smoothness can be promoted by increasing $η$ and $λ$ , respectively. $γ$ and $t$ are related to the step sizes in the algorithms, which were tuned in a trial-and-error way for faster and more stable convergence.

Table 3

The runtime and the hyperparameters used in all cases.

		Digital phantom	Catphan	Quantitative imaging phantom	Pelvis patient
$η$	$α = 0$	0.03	0.006	0.0009	0.003
	$α = 1 / 2$	0.05	0.002	0.0009	0.001
	$α = 2 / 3$	0.1	0.003	0.0009	0.002
	$α = 1$	0.5	0	0.1	0
$λ$	$α = 0$	0.05	0.002	0.02	0.0007
	$α = 1 / 2$	0.25	0.003	0.005	0.0012
	$α = 2 / 3$	0.05	0.0025	0.003	0.0005
	$α = 1$	0.01	0.0005	0	0.0005
$γ$	$α = 0$	0.01	0.0001	1.00E-05	0.0001
	$α = 1 / 2$	0.01	0.0001	1.00E-05	0.0001
	$α = 2 / 3$	0.001	0.001	0.001	1.00E-05
	$α = 1$	0.1	0.001	1	0.001
$t$	$α = 0$	100	1	1	1
	$α = 1 / 2$	100	1	1	1
	$α = 2 / 3$	100	100	0.5	100
	$α = 1$	100	100	0.01	100
Runtime (s)	$α = 0$	63.4	791.5	861.7	738.8
	$α = 1 / 2$	71.2	881.4	489.1	1640.1
	$α = 2 / 3$	81.5	1059.4	3542.8	962.1
	$α = 1$	81.3	95.1	323.9	294.7
	DI	36.2	133.6	490.1	30.9

Figure 9 shows the convergence plots for different $α$ values on the Catphan and the pelvis patient case. The objective values for different cases were not comparable since different hyperparameter values were used in different cases. However, it is worth noting that the plots show different converging patterns. For the convex case $α = 1$ , the objective goes down nicely with a convex-shaped convergence curve, showing stable and robust convergence. For the nonconvex case $α = 0$ , the curve goes down with sudden changes. For the $α$ values in between, the $α = \frac{2}{3}$ curve patterns are more similar to that of the $α = 1$ curves, whereas the $α = \frac{1}{2}$ curves are more irregular and bumpy.

Fig. 9

Convergence plots for different $α$ values on the Catphan and the pelvis patient case.

4. Discussions

The standard DI method for DECT-based MMD imposes sparsity constraint by enforcing a hard ceiling of three on the number of materials in each voxel and then solves the basis material components by direct matrix inversion. However, the three-material constraint is arbitrary and unrealistically rigid. Moreover, the direct matrix inversion inevitably amplifies the image noise, as shown in previous publications¹⁰ and the current study. Our proposed DECT MMD framework utilizes a TV regularization that regulates the decomposition image noise and uses a sparsity regularization to penalize the number of materials that are simultaneously present at the same voxel. The soft sparsity regularization allows the number of basis materials to vary across different voxels.

The sparsity term is in the form of ${| x |}^{α}$ ( $0 \leq α \leq 1$ ), where $α = 0, \frac{1}{2}, \frac{2}{3}, 1$ were specifically studied. The sparsity term reduces to the L1 norm of the $x$ with $α = 1$ and further reduces to a constant under the mass and volume constraint. Therefore, the sparsity term with $α = 1$ does not promote material sparsity despite its desirable mathematical properties of being convex. When $α = 0, \frac{1}{2}, \frac{2}{3}$ , the corresponding sparsity term has a closed-form proximal operator, which could be difficult to evaluate for other values of $α$ . The differences between $α = 0$ , $α = \frac{1}{2}$ , and $α = \frac{2}{3}$ are subtle with respect to the decomposition results. The material penalty term with $α = 0$ is also referred to as the counting norm and a mathematically rigorous description of the material sparsity. This is reflected in the best quantitative performance achieved with this norm. However, due to its extreme nonconvexity, it is more challenging to tune the parameters for convergence to an acceptable local minimum. The sparsity terms with $α = \frac{1}{2}$ and $α = \frac{2}{3}$ are better behaving nonconvex functions, showing more robust performance to optimization parameter selection.

The nonconvex optimization problem in this study was solved using an accelerated primal-dual algorithm with line search proposed in Ref. 12. The primal-dual algorithm reduces to the primal-dual hybrid gradient (PDHG) method¹⁶^,¹⁷ when applied to convex optimization problems. PDHG is one of the first-ordered algorithms that are well-suited for large-scaled optimization problems. Other first-ordered algorithms such as the fast iterative shrinkage-thresholding algorithm,¹⁸ which achieves a convergence rate of $O (\frac{1}{k^{2}})$ , substantially faster than the $O (\frac{1}{k})$ convergence rate of the PDHG, will be investigated in future studies on the nonconvex optimization problems.

5. Conclusions

The proposed method accurately decomposed the DECT phantom and patient images up to 12 basis materials, markedly reduced cross talk among materials, suppressed decomposition image noise, and retained image spatial resolution. The proposed method using nonconvex material sparsity penalty outperforms convex penalty and the standard DI method.

6. Appendix

6.1.

Evaluate the Proximal Operator of $η {| \cdot |}^{α}$ for Different $α$ Values

For $α = 0$ , the proximal operator in Eq. (5) reduces to a hard thresholding operation:

{({Prox}_{w F_{3}} ({\hat{z}}_{3}))}_{i} = {\begin{matrix} {\hat{z}}_{3 i}, & if | {\hat{z}}_{3 i} | \geq \sqrt{2 η w} \\ 0, & if | {\hat{z}}_{3 i} | < \sqrt{2 η w} \end{matrix} .

For $α = \frac{1}{2}$ , the closed-form solution was proposed by McKelvey:¹⁹

{({Prox}_{w F_{3}} ({\hat{z}}_{3}))}_{i} = {\begin{cases} \frac{4}{3} \sin^{2} (\frac{1}{3} \arccos (\frac{3 \sqrt{3}}{4} u) + \frac{π}{2}), & if u \leq \frac{2 \sqrt{6}}{9} \\ 0, & if u > \frac{2 \sqrt{6}}{9} \end{cases},

where

u = η w {| {\hat{z}}_{3 i} |}^{- \frac{3}{2}}

.

For $α = \frac{2}{3}$ , the solution of the proximal operator is a root of the quartic polynomial:

Eq. (7)

x^{4} - 3 {\hat{z}}_{3 i} x^{3} + 3 {\hat{z}}_{3 i}^{2} x^{2} - {\hat{z}}_{3 i}^{3} x + \frac{8 η^{3} w^{3}}{27} = 0,

for which the analytic solution exists and could be evaluated following the work in Ref. 20. The details on the proximal operator evaluation for

α = \frac{2}{3}

are shown in Algorithm 3.

Algorithm 3

Proximal operator evaluation (α=2/3), the solution to Eqs. (5) and (7).

Compute intermediary terms

a_{1}

,

a_{2}

,

a_{3}

,

a_{4}

:

a_{1} ≔ \frac{8 w^{3} η^{3}}{27 {\hat{z}}_{3 i}^{4}}, a_{2} ≔ \sqrt[3]{\frac{a_{1}}{16} + \sqrt{- \frac{a_{1}^{3}}{27} + \frac{a_{1}^{2}}{256}}}, a_{3} ≔ \frac{5}{8} + 2 a_{2} + \frac{2 a_{1}}{3 a_{2}}, a_{4} ≔ \sqrt{a_{3} - \frac{3}{8}}

Compute 4 roots

r_{1}

,

r_{2}

,

r_{3}

,

r_{4}

:

r_{1} ≔ \frac{3}{4} + (a_{4} + \sqrt{\frac{9}{8} - a_{3} - \frac{1}{4 a_{4}}}) / 2, r_{2} ≔ \frac{3}{4} + (a_{4} - \sqrt{\frac{9}{8} - a_{3} - \frac{1}{4 a_{4}}}) / 2 r_{3} ≔ \frac{3}{4} + (- a_{4} + \sqrt{\frac{9}{8} - a_{3} + \frac{1}{4 a_{4}}}) / 2, r_{4} ≔ \frac{3}{4} + (- a_{4} - \sqrt{\frac{9}{8} - a_{3} + \frac{1}{4 a_{4}}}) / 2

Pick the root

r

as the maximum of

{0, r_{1}, r_{2}, r_{3}, r_{4}}

that satisfies

1.

r

is a real number

2.

0.5 < r < 1

Output:

{({Prox}_{{w F}_{3}} ({\hat{z}}_{3}))}_{i} = r {\hat{z}}_{3 i}

For $α = 1$ , the proximal operator results in soft thresholding:

{({Prox}_{w F_{3}} ({\hat{z}}_{3}))}_{i} = {\begin{cases} {\hat{z}}_{3 i} - η, & if {\hat{z}}_{3 i} \geq η, \\ 0, & if | {\hat{z}}_{3 i} | < η, \\ {\hat{z}}_{3 i} + η, & if {\hat{z}}_{3 i} \leq - η \end{cases} .

6.2.

LAC Values

Table 4 presents the HE LAC and LE LAC values of the basis materials, which were computed as the average values of the HE LACs and LE LACs within the ROIs shown on the LE CT images in Figs. 1, 3, 5, and 7.

Table 4

The HE LAC and LE LAC values of the basis materials in all cases.

ROI #	Digital phantom		Catphan		Quantitative imaging phantom		Pelvis patient
ROI #	LE (mm−1)	HE (mm−1)	LE (mm−1)	HE (mm−1)	LE (mm−1)	HE (mm−1)	LE (mm−1)	HE (mm−1)
1	0.059	0.033	0.043	0.038	0.024	0.019	0.059	0.035
2	0.045	0.022	0.029	0.027	0.026	0.020	0.031	0.021
3	0.023	0.017	0.033	0.027	0.029	0.021	0.025	0.019
4	0.001	0.000	0.018	0.017	0.031	0.022	0.020	0.016
5			0.024	0.022	0.027	0.022	0.000	0.000
6			0.003	0.003	0.030	0.024
7					0.043	0.030
8					0.025	0.020
9					0.025	0.021
10					0.021	0.018
11					0.023	0.019
12					0.023	0.019
13					0.000	0.000

Disclosures

No conflicts of interest, financial or otherwise, are declared by the authors.

Acknowledgments

This research is supported by DOE under Grant Nos. DE-SC0017057 and DE-SC0017687, and NIH under Grant Nos. R44CA183390, R01CA188300, R01CA230278, R21CA228160, and R43CA183390. We would like to thank Shuai Leng at Mayo Clinic for providing the quantitative imaging phantom DECT data. An earlier version of this work was published in the conference proceedings of SPIE Medical Imaging 2019.²¹

References

1.

M. D. Agrawal et al., “Oncologic applications of dual-energy CT in the abdomen,” RadioGraphics, 34 (3), 589 –612 (2014). https://doi.org/10.1148/rg.343135041 Google Scholar

2.

H. A. Lee et al., “Comparison of virtual unenhanced images derived from dual-energy CT with true unenhanced images in evaluation of gallstone disease,” Am. J. Roentgenol., 206 (1), 74 –80 (2016). https://doi.org/10.2214/AJR.15.14570 AJROAM 0092-5381 Google Scholar

3.

Q. Wang et al., “Quantitative analysis of the dual-energy CT virtual spectral curve for focal liver lesions characterization,” Eur. J. Radiol., 83 (10), 1759 –1764 (2014). https://doi.org/10.1016/j.ejrad.2014.07.009 EJRADR 0720-048X Google Scholar

4.

D. T. Boll et al., “Renal stone assessment with dual-energy multidetector CT and advanced postprocessing techniques: improved characterization of renal stone composition—pilot study,” Radiology, 250 (3), 813 –820 (2009). https://doi.org/10.1148/radiol.2503080545 RADLAX 0033-8419 Google Scholar

5.

C. N. De Cecco et al., “Dual-energy CT: oncologic applications,” Am. J. Roentgenol., 199 (Suppl. 5), S98 –S105 (2012). https://doi.org/10.2214/AJR.12.9207 AJROAM 0092-5381 Google Scholar

6.

R. E. Alvarez and A. MacOvski, “Energy-selective reconstructions in x-ray computerized tomography,” Phys. Med. Biol., 21 733 –744 (1976). https://doi.org/10.1088/0031-9155/21/5/002 PHMBA7 0031-9155 Google Scholar

7.

A. Macovski et al., “Energy dependent reconstruction in x-ray computerized tomography,” Comput. Biol. Med., 6 325 –334 (1976). https://doi.org/10.1016/0010-4825(76)90069-X CBMDAW 0010-4825 Google Scholar

8.

C. H. McCollough et al., “Dual- and multi-energy CT: principles, technical approaches, and clinical applications,” Radiology, 276 (3), 637 –653 (2015). https://doi.org/10.1148/radiol.2015142631 RADLAX 0033-8419 Google Scholar

9.

Y. Long and J. A. Fessler, “Multi-material decomposition using statistical image reconstruction for spectral CT,” IEEE Trans. Med. Imaging, 33 (8), 1614 –1626 (2014). https://doi.org/10.1109/TMI.2014.2320284 ITMID4 0278-0062 Google Scholar

10.

P. R. S. Mendonça, P. Lamb and D. V. Sahani, “A flexible method for multi-material decomposition of dual-energy CT images,” IEEE Trans. Med. Imaging, 33 (1), 99 –116 (2014). https://doi.org/10.1109/TMI.2013.2281719 ITMID4 0278-0062 Google Scholar

11.

T. Moellenhoff et al., “Low rank priors for color image regularization,” Lect. Notes Comput. Sci., 8932 126 –140 (2014). https://doi.org/10.1007/978-3-319-14612-6 LNCSD9 0302-9743 Google Scholar

12.

Y. Malitsky and T. Pock, “A first-order primal-dual algorithm with linesearch,” SIAM J. Optim., 28 (1), 411 –432 (2018). https://doi.org/10.1137/16M1092015 SJOPE8 1095-7189 Google Scholar

13.

N. Parikh and S. Boyd, “Proximal algorithms,” Found. Trends Optim., 1 (3), 123 –231 (2013). https://doi.org/10.1561/2400000003 Google Scholar

14.

W. Wang and M.-Á. Carreira-Perpiñán, “Projection onto the probability simplex: an efficient algorithm with a simple proof, and an application,” (2013). Google Scholar

15.

S. Boyd and L. Vandenberghe, “Convex optimization,” Cambridge University Press, Cambridge (2004). Google Scholar

16.

T. Pock et al., “An algorithm for minimizing the Mumford-Shah functional,” in Proc. IEEE Int. Conf. Comput. Vision, (2009). https://doi.org/10.1109/ICCV.2009.5459348 Google Scholar

17.

A. Chambolle and T. Pock, “A first-order primal-dual algorithm for convex problems with applications to imaging,” J. Math. Imaging Vision, 40 (1), 120 –145 (2011). https://doi.org/10.1007/s10851-010-0251-1 JIMVEC 0924-9907 Google Scholar

18.

A. Beck and M. Teboulle, “A fast iterative shrinkage-thresholding algorithm for linear inverse problems,” SIAM J. Imaging Sci., 2 (1), 183 –202 (2009). https://doi.org/10.1137/080716542 Google Scholar

19.

J. P. McKelvey, “Simple transcendental expressions for the roots of cubic equations,” Am. J. Phys., 52 (3), 269 –270 (1984). https://doi.org/10.1119/1.13706 AJPIAS 0002-9505 Google Scholar

20.

D. Krishnan, R. Fergus, “Fast image deconvolution using hyper-Laplacian priors,” in Proc. 22nd Int. Conf. Neural Inf. Process. Syst., 1033 –1041 (2009). Google Scholar

21.

Q. Lyu et al., “Image-domain multi-material decomposition for dual-energy CT with non-convex sparsity regularization,” Proc. SPIE, 10949 1094903 (2019). https://doi.org/10.1117/12.2508037 PSISDG 0277-786X Google Scholar

Biography

Qihui Lyu received her BS degree in physics from Nanjing University in Nanjing, China. Currently, she is a PhD candidate in the Physics and Biology in Medicine Program at University of California, Los Angeles. She is interested in advanced optimization algorithms for inverse treatment planning, medical image reconstruction, and dual-energy CT.

Daniel O’Connor received his PhD in mathematics from University of California, Los Angeles. He is an assistant professor at University of San Francisco. His postdoctoral research in the UCLA Department of Radiation Oncology focused on applications of optimization in radiation therapy, medical imaging, and machine learning. His research interests are in large scale optimization modeling and algorithms.

Tianye Niu received his PhD in physics and electronics from the University of Science and Technology of China. He is a professor at the Institute of Translational Medicine at Zhejiang University. His research interests include compressed sensing, low-dose reconstruction, dual-energy CT imaging, spectral CT imaging, IGRT, proton therapy, radiomics, and cone-beam CT instrument.

Ke Sheng received his PhD in medical physics from the University of Wisconsin–Madison. He is a professor and an associate vice chair of radiation oncology at University of California, Los Angeles. His main research areas include inverse optimization, 4π noncoplanar radiotherapy, image guided radiotherapy, biological outcome modeling, and small animal irradiation.

Citation Download Citation

Qihui Lyu, Daniel O'Connor, Tianye Niu, and Ke Sheng "Image-domain multimaterial decomposition for dual-energy computed tomography with nonconvex sparsity regularization," Journal of Medical Imaging 6(4), 044004 (15 October 2019). https://doi.org/10.1117/1.JMI.6.4.044004

Received: 8 April 2019; Accepted: 23 September 2019; Published: 15 October 2019

Access the abstract

JOURNAL ARTICLE
12 PAGES

DOWNLOAD PAPER SAVE TO MY LIBRARY

GET CITATION

CITATIONS

Cited by 6 scholarly publications.

Explore citations on Lens.org

RIGHTS & PERMISSIONS

Get copyright permission Get copyright permission on Copyright Marketplace

KEYWORDS

Computed tomography

Iodine

Bone

Calcium

Blood

Tissues

Nickel

Show All Keywords

CHORUS Article. This article was made freely available starting 14 October 2020

1.

Introduction

2.

Method

2.1.

Formulation

Eq. (1)

Eq. (2)

2.2.

Algorithm

Eq. (3)

Eq. (4)

Algorithm 1

Algorithm 2

Eq. (5)

2.3.

Evaluation

Fig. 1

Eq. (6)

3.

Results

3.1.

Digital Phantom

Fig. 2

3.2.

Catphan® 600 Phantom

Fig. 3

Fig. 4

Table 1

3.3.

Quantitative Imaging Phantom

Fig. 5

Fig. 6

3.4.

Pelvis Patient

Fig. 7

Fig. 8

Table 2

Table 3

Fig. 9

4.

Discussions

5.

Conclusions

6.

Appendix

6.1.

Evaluate the Proximal Operator of η|·|α for Different α Values

Eq. (7)

Algorithm 3

6.2.

LAC Values

Table 4

Disclosures

Acknowledgments

References

Biography

Show All Keywords

Keywords/Phrases

Search In:

Publication Years

Catphan^® 600 Phantom

Evaluate the Proximal Operator of $η {| \cdot |}^{α}$ for Different $α$ Values