Image domain multi-material decomposition using single energy CT

In the image domain MMD theory, the mass attenuation of a mixture ${\mu _M}\left( E \right)$ in the CT image is assumed to be a linear combination of the mass attenuation of different basis materials and written as (Long and Fessler 2014):

$$\begin{equation}{\mu _M}\left( E \right) = \mathop \sum \limits_{t = 1}^{{T_0}} {\beta _t}{\mu _{tM}}\left( E \right)\!\!,{\,}s.{\,}t.{\,}\mathop \sum \limits_{t = 1}^{{T_0}} {\beta _t} = 1,{\,}{\beta _t} \geqslant 0\end{equation} \tag{ 1 }$$

where ${\mu _{tM}}\left( E \right)$ is the mass attenuation coefficient of the tth basis material at the energy level E. The subscript M indicates the mass attenuation coefficient. ${T_0}$ represents the total number of basis materials. ${\beta _t}$ is the mass fraction of the tth basis material and is written as:

$$\begin{equation}{\beta _t} = \frac{{{m_t}}}{{\mathop \sum \nolimits_{t = 1}^{{T_0}} {m_t}}}\end{equation} \tag{ 2 }$$

where ${m_t}$ is the mass of the tth basis material.

The relationship between the linear attenuation coefficient (LAC)${\,}{\mu _E}{\,}$ and mass attenuation coefficient ${\mu _M}\,$ is expressed as:

$$\begin{equation}{\mu _E} = \,\rho {\mu _M}\left( E \right) = \frac{{\mathop \sum \nolimits_{t = 1}^{{T_0}} {m_t}}}{{\mathop \sum \nolimits_{t = 1}^T {V_t}}}{\mu _M}\left( E \right),\end{equation} \tag{ 3 }$$

where ${V_t}$ is the volume of the tth basis material.

Combining equations (1) and (2) into equation (3), we have:

$$\begin{equation}{\mu _E} = \mathop \sum \limits_{t = 1}^{{T_0}} \left(\frac{{\mathop \sum \nolimits_{t = 1}^{{T_0}} {m_t}}}{{\mathop \sum \nolimits_{t = 1}^{{T_0}} {V_t}}}\frac{{{m_t}}}{{\mathop \sum \nolimits_{t = 1}^{{T_0}} {m_t}}}\frac{{{V_t}{\mu _{tE}}}}{{{m_t}}}\right) = \mathop \sum \limits_{t = 1}^{{T_0}} {\mu _{tE}}\ {x_t}\end{equation} \tag{ 4 }$$

where ${x_t} = \frac{{{V_t}}}{{\mathop \sum \nolimits_{t = 1}^{{T_0}} {V_t}}}$ represents the volume fraction of the tth basis material and ${\mu _{tE}}$ is the LAC of the tth basis material. ${x_t}$ satisfies the sum-to-one and box [0 1] constraints:

$$\begin{equation}\sum\limits_{t = 1}^{{T_0}} {{x_t}} = \sum\limits_{t = 1}^{{T_0}} {\frac{{{V_t}}}{{\sum\nolimits_{t = 1}^{{T_0}} {{V_t}} }}} = 1,{\mkern 1mu} 0 \le {x_t} \le 1,{\mkern 1mu} t = 1, \ldots ,{T_0}\end{equation} \tag{ 5 }$$

The MMD model can be constructed from equations (4) and (5) as:

$$\begin{equation}{\mu _E} = \mathop \sum \limits_{t = 1}^{{T_0}} {x_t}{\mu _{tE}},\,s.\,t.\,\mathop \sum \limits_{t = 1}^{{T_0}} {x_t} = 1,\,{x_t} \ge 0,\forall t\end{equation} \tag{ 6 }$$

In the above equation (6), ${T_0}$ unknows need to be estimated from two equations including a CT measurement ${\mu _E}$ combined from basis material LAC equation and a sum-to-one equation. Additional assumptions are needed to solve the ${T_0}$ unknows. We assume that each pixel in the CT image contains no more than two materials, i.e.

$$\begin{equation}\mathop \sum \limits_{t = 1}^{{T_0}} {I_{\left\{ {{x_{tp}} \ne 0} \right\}}} \le 2,\forall p\end{equation} \tag{ 7 }$$

where ${I_{\left\{ \cdot \right\}}}$ denotes an indicator function which equals 1 if the tth material exists and 0 otherwise. p denotes the pth pixel in the CT image. We define $\Omega$ as a two-material candidates library composed of all the two-material candidates formed from ${T_0}$ basis materials. For each two-material candidate in the library $\Omega$, only two materials need to be estimated from equation (6), indicating that the equation (6) is solvable with the assumption equation (7).

To implement the MMD in the CT image of N_p pixels, equation (6) is rewritten in the matrix–vector form as:

$$\begin{equation}\vec \mu = A\vec x,\,s.t.\,\left\{ \begin{gathered} \mathop \sum \limits_{t = 1}^{{T_0}} {x_t} = 1,\,{x_t} \ge 0,\forall t \hfill \\ \mathop \sum \limits_{t = 1}^{{T_0}} {I_{\left\{ {{x_{tp}} \ne 0} \right\}}} \le 2,\forall p \hfill \\ \end{gathered} \right.,\end{equation} \tag{ 8 }$$

where $\vec \mu$ is the measured CT images after column vectorization. Similarly, $\vec x = {\left[ {{{\overrightarrow {{x_1}} }^T},{{\overrightarrow {{x_2}} }^T}, \ldots ,{{\overrightarrow {{x_{{T_0}}}} }^T}} \right]^T}$ is the columnized vector with the size of ${T_0} \times {N_p}$ and each entry as the volume fraction image of the tth basis material. $A \in {\mathbb{R}^{{N_p} \times {T_0}{N_p}}}$ is the material composition matrix and is written as:

$$\begin{equation}A = {A_0} \otimes {I_{{N_p}}}\end{equation} \tag{ 9 }$$

where $\otimes$ represents the Kronecker product. ${I_{{N_p}}}$ is a ${N_p} \times {N_p}$ identity matrix. ${A_0} = \left[ {{\mu _1},{\mu _2}, \ldots {\mu _{{T_0}}}} \right]$ is the composition matrix composed of the LACs of the basis materials. The matrix ${A_0}$ is determined by selecting a uniform region of interest (ROI) on the CT images that contain the basis materials and taking the mean value of the ROI as the entry of matrix A.

We formulate the image domain SECT MMD using the least square estimation with regularization as following:

$$\begin{equation}\begin{gathered} \widehat {\vec x} = _{\vec x\,\,subject\,\,to\,\,\left( 5 \right)\left( 7 \right)}^{\,\,\,\,\,\,\,\,\,\,\,\,\,argmin}\frac{1}{2}\left\| {A\vec x - \vec \mu } \right\|_2^2 + R\left( {\vec x} \right) \hfill \\ \hfill \\ \end{gathered} \end{equation} \tag{ 10 }$$

The material-wise edge-preserving regularization, which has distinguished noise suppression ability and edge maintainability, is selected as $R\left( {\vec x} \right)$ (Fessler and Rogers 1996, Long and Fessler 2014):

$$\begin{equation}R\left( {\vec x} \right) = \mathop \sum \limits_{t = 1}^{{T_0}} {\beta _t}{R_t}\left( {{{\vec x}_t}} \right),\end{equation} \tag{ 11 }$$

where the regularizer for the tth material is

$$\begin{equation}{R_t}\left( {{{\vec x}_t}} \right) = \mathop \sum \limits_{p = 1}^{{N_P}} \mathop \sum \limits_{k \in {N_{tp}}} {\psi _t}\left( {{x_{tp}} - {x_{tk}}} \right),\end{equation} \tag{ 12 }$$

where $k$ is the kth neighboring pixel of p. The potential function ${\psi _t}$ is a hyperbola:

$$\begin{equation}{\psi _t}\left( z \right) = \frac{{\delta _t^2}}{3}\left( {\sqrt {1 + 3{{\left( {\frac{z}{{{\delta _t}}}} \right)}^2}} - 1} \right).\end{equation} \tag{ 13 }$$

The regularization parameters ${\beta _t}$ and ${\delta _t}$ are chosen for different materials individually to achieve desired edge preservation and noise-resolution tradeoff for each material image. We refer to this proposed MMD scheme as the edge-preserving (EP) constrained SECT method hereafter.

The objective function equation (10) cannot be minimized directly since the two-material assumption is a non-convex constraint. Similar to our previous work, we apply the optimization transfer principles to minimize to find a series of PWSQS functions at each iteration to decrease the cost function monotonically. The PWSQS function of the cost function at the nth iteration is written as (Xue et al 2017):

$$\begin{equation}{\phi ^{\left( n \right)}}\left( {\vec x} \right) = \mathop \sum \limits_{p = 1}^{{N_p}} \left( {\frac{1}{2}\left\| {{A_0}\overrightarrow {{x_p}} - \overrightarrow {{\mu _p}} } \right\|_2^2 + R_p^{\left( n \right)}\left( {{{\vec x}_p}} \right)} \right)\end{equation} \tag{ 14 }$$

where $R_p^{\left( n \right)}\left( {{{\vec x}_p}} \right)$ is the PWSQS function of the regularization term. Applying the De Pierro's additive convexity trick and Huber's optimal curvature as our previous work, we further derive the PWSQS function ${\phi ^{\left( n \right)}}\left( {\vec x} \right)$ as (De Pierro 1995, Erdogan and Fessler 1999, Fessler 2000, Huber 2011, Xue et al 2017):

$$\begin{equation}{\phi ^{\left( n \right)}}\left( {\vec x} \right) = \mathop \sum \limits_{p = 1}^{{N_p}} \left( {\frac{1}{2}\vec x_p^TH{{\vec x}_p} + {{\vec q}^T}{{\vec x}_p}} \right)\end{equation} \tag{ 15 }$$

where

$$\begin{equation}H = {\,}2A_0^TV_P^{\,- 1}{A_0}\, + H_{Rp}^{\left( n \right)},\end{equation} \tag{ 16 }$$

$$\begin{equation}\vec q = \,2A_0^TV_P^{\,- 1}{A_0}\vec x_p^{\left( n \right)} - 2A_0^TV_P^{\,- 1}{\vec \mu _p} + \,\dot R_p^{\left( n \right)} - {H^T}\vec x_p^{\left( n \right)}\end{equation} \tag{ 17 }$$

where $H_{Rp}^{\left( n \right)}\,$ and $\dot R_p^{\left( n \right)}$ are the Hessian and gradient of the regularization term and is derived as:

$$\begin{equation}\,H_{Rp}^{\left( n \right)} \triangleq diag\left\{ {4{\beta _t}\mathop \sum \limits_{k \in {N_{tp}}} {\omega _{{\psi _t}}}\left( {x_{tp}^{\left( n \right)} - x_{tk}^{\left( n \right)}} \right)} \right\},\end{equation} \tag{ 18 }$$

where${\,}\,{\omega _{{\psi _t}}}\left( z \right) \triangleq {\dot \psi _t}\left( z \right)/z$.

$$\begin{equation}\dot R_p^{\left( n \right)} = {\left[ {{\beta _1}\frac{\partial }{{\partial {x_{1p}}}}{R_1}\left( {\vec x_1^{\left( n \right)}} \right), \ldots ,{\beta _{{T_0}}}\frac{\partial }{{\partial {x_{{T_0}p}}}}{R_{{T_0}}}\left( {\vec x_{{T_0}}^{\left( n \right)}} \right)} \right]^T}\end{equation} \tag{ 19 }$$

where

$$\begin{equation}\frac{\partial }{{\partial {x_{tp}}}}{R_t}\left( {\vec x_t^{\left( n \right)}} \right) = \mathop \sum \limits_{k \in {N_{tp}}} {\dot \psi _t}\left( {x_{tp}^{\left( n \right)} - x_{tk}^{\left( n \right)}} \right),t = 1, \ldots ,{T_0}\end{equation} \tag{ 20 }$$

To satisfy the constraint of each pixel containing at most two material, we loop over all the two-material candidate in the two-material candidate library ${\Omega }$. The candidate with the minimal surrogate objective function (i.e. equation (14)) value is determined to be the optimal two-material candidate for the pixel to be decomposed. For each two-material candidate, $\tau = \left( {i,j} \right) \in {\Omega }$, the surrogate function degenerates to a classical convex quadratic function of a vector with two unknowns, ${\vec x_p}\left( \tau \right) \triangleq {\left[ {{x_{ip,\,}}{x_{jp\,}}} \right]^T}$. Combining the sum-to-one and box [0 1] constraints, the surrogate objective function can be further rewritten as:

$$\begin{equation}\begin{gathered} {\widehat {\vec x}_p}\left( \tau \right) = \,\,_{{{\vec x}_p}\left( \tau \right)\,subject\,to\,\left( 5 \right)\,\& \,\left( 7 \right)}^{\qquad\,argmin}\,\phi _p^{\left( n \right)}\left( {{{\vec x}_p}\left( \tau \right)} \right) \hfill \\ \phi _p^{\left( n \right)}\left( {{{\vec x}_p}\left( \tau \right)} \right) = \frac{1}{2}\vec x_p^T\left( \tau \right)H\left( \tau \right){{\vec x}_p}\left( \tau \right) + {{\vec q}^T}\left( \tau \right){{\vec x}_p}\left( \tau \right) \hfill \\ s.t.\left\{ {\begin{array}{*{20}{c}} {\mathop \sum \limits_{t = 1}^{{T_0}} {x_{tp}} = 1,} \\ {0 \le {x_{tp}} \le 1.} \end{array}} \right. \hfill \\ \end{gathered} \end{equation} \tag{ 21 }$$

where $H\left( \tau \right)$ and $\vec q\left( \tau \right)\,$ are formed from elements in $H\,$ and $\vec q$ with indexes corresponding to $\tau = \left( {i,j} \right)$, respectively. The above convex quadratic programming problem can be solved utilizing the Generalized Sequential Minimization Algorithm (GSMO).

Due to the nonconvex of objective function equation (10), a good initial guess is needed to initialize the minimization. Here we propose a direct inversion method based on the two-material assumption to realize the initialization. Equation (8) for the pth pixel can be re-written as:

$$\begin{equation}{A_0}\overrightarrow {{x_p}} = \overrightarrow {{\mu _p}} ,s.t.\,\mathop \sum \limits_{t = 1}^{{T_0}} {I_{\left\{ {{x_{tp}} \ne 0} \right\}}} \le 2,\,\mathop \sum \limits_{t = 1}^{{T_0}} {x_{pt}} = 1,\,{x_{pt}} \ge 0,\forall t\end{equation} \tag{ 22 }$$

where ${A_0} = \left[ {{\mu _1},{\mu _2}, \ldots {\mu _{{T_0}}}} \right]$ is a $1 \times {T_0}$ matrix. When combining the two-material assumption and the sum-to-one constraint with ${A_0}$, we have the composition matrix $A'$ for a two-material candidate $\tau = \left( {i,j} \right)$ as:

$$\begin{equation}A' = \,\left[ {\begin{array}{*{20}{c}} {{\mu _i}}&{{\mu _j}} \\ 1&1 \end{array}} \right]\end{equation} \tag{ 23 }$$

where ${\mu _i}$ and ${\mu _j}$ are the LACs of the ith and jth basis materials, respectively. The two-material decomposition of the pth pixel can be written as:

$$\begin{equation}\left[ {\begin{array}{*{20}{c}} {{\mu _i}}&{{\mu _j}} \\ 1&1 \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{x_{pi}}} \\ {{x_{pj}}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {{\mu _p}} \\ 1 \end{array}} \right],\,s.t.\,{x_{pi}} \ge 0,{x_{pj}} \ge 0.\end{equation} \tag{ 24 }$$

Equation (24) is initially solved using the direct matrix inversion since the determinant of $A'$ is not equal to zero. The MMD can be achieved by looping over all the two-material candidates with the total number of $C_{{T_0}}^{\,2}$ options in to two-material candidate library $\Omega$ and picking up the optimal solution with minimum Euclidean distance between the LACs of the current pixel and the basis materials. The distance is calculated as:

$$\begin{equation}{\tau ^*} = \arg {\min _{\tau = \left( {i,j} \right)}}\sqrt {{{\left\| {{\mu _p} - {\mu _i}} \right\|}^2} + {{\left\| {{\mu _p} - {\mu _j}} \right\|}^2}} \end{equation} \tag{ 25 }$$

If no direct solution to equation (24) exists, the optimal one is estimated by minimizing the least-square form of equation (24) as:

$$\begin{equation}\left( {{\tau ^*},{x_{pi}}^*{x_{pj}}^*} \right) = arg\mathop {{\text{min}}}\limits_{\begin{array}{*{20}{c}} {\tau = \left( {i,j} \right)} \\[-3pt] {{x_{pi}},{x_{pj}}} \end{array}} {\left\| {\left[ {\begin{array}{*{20}{c}} {{\mu _i}}&{{\mu _j}} \\ 1&1 \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{x_{pi}}} \\ {{x_{pj}}} \end{array}} \right] - \left[ {\begin{array}{*{20}{c}} {{\mu _p}} \\ 1 \end{array}} \right]} \right\|_2},\,s.t.\,\left\{ {\begin{array}{*{20}{c}} {{x_{pi}} \ge 0} \\ {{x_{pj}} \ge 0} \end{array}} \right.\end{equation} \tag{ 26 }$$

Equation (26) can be solved using, e.g. the gradient projection (GP) algorithm (Figueiredo et al 2007). We refer to the proposed SECT initialization scheme as the two-material assumption (TMA) constrained SECT method hereafter.

In summary, the pseudo code of the proposed method is shown in table 1, where the CT images and decomposed material images are denoted by vector signs. The initial value is generated using the proposed TMA method as shown in Lines 1–22. In Lines 11–14 of the multi feasible solutions case, the optimal material candidate is selected by minimizing Euclidean distance between the current pixel and the basis materials. In Lines 15–19 of the no feasible solution case, the optimal solution is estimated by minimizing equation (26). Lines 23–35 show the proposed EP method. The Hessian and gradient of the regularization term are computed in Line 25. The two-material candidate looping operation is performed in Lines 26 to 30. The two-material specific convex quadratic programming problem is solved in Line 29. The optimal result is selected in Line 31.

Table 1. Pseudocode of the proposed MMD algorithm for SECT.

line	program
1	***Initial value generation using the proposed TMA method***
2	//${\mu _p}$ is the LAC of the pth pixel; P is the number of pixels. ${\tau _c}$ is the c-th two-material candidate. C is the number of material candidates. (i,j) is the ith and jth material basis in the material candidate ${\tau _c}$.
3	for p = 1:P
4	count = 0; //the variable count counts the number of feasible solutions
5	for c = 1:C
6	$\left[\!\! {\begin{array}{*{20}{c}} {{x_{pi}}} \\ {{x_{pj}}} \end{array}}\!\! \right] = {\left[\!\! {\begin{array}{*{20}{c}} {{\mu _i}}\quad {{\mu _j}} \\ 1 \quad 1 \end{array}} \right]^{ - 1}}\left[\!\! {\begin{array}{*{20}{c}} {{\mu _p}} \\ 1 \end{array}} \!\!\right]$
7	if ${x_{pi}},{x_{pj}} \ge 0$
8	count = count ± 1;
9	end
10	end
11	if count > 0//Multi solution exist
12	//Calculate the Euclidean distance
13	${\left. {{\tau _c}^{\text{*}}} \right|_{{\tau _c}^{\text{*}} = \left( {{i^*},{j^*}} \right)}} = \arg \mathop {\min }\limits_{{\tau _c} = \left( {i,j} \right)} \sqrt {{{\left\| {{\mu _p} - {\mu _i}} \right\|}^2} + {{\left\| {{\mu _p} - {\mu _j}} \right\|}^2}}$;
14	$x_p^{\text{*}} = {\left( {{x_{p{i^{\text{*}}}}},{x_{p{j^{\text{*}}}}}} \right)^T}$;
15	else //no feasible solution
16	Use GP method to solve:
17	$\left( {{\tau _c}^{\rm{*}},{x_{pi}}^{\rm{*}}{x_{pj}}^{\rm{*}}} \right) = arg\!\!\!\mathop{{\rm{min}}}\limits_{\begin{array}{*{20}{c}} {{\tau _c} = \left( {i,j} \right)}\\ {{x_{pi}},{x_{pj}}} \end{array}} \left[\!\! {\begin{array}{*{20}{c}} {{\mu _i}}\quad{{\mu _j}}\\ 1 \quad 1 \end{array}} \!\!\right]\left[\!\! {\begin{array}{*{20}{c}} {{x_{pi}}}\\ {{x_{pj}}} \end{array}} \!\!\right] - \left[ {\begin{array}{*{20}{c}} {{\mu _p}}\\ 1 \end{array}} \!\!\right]{_2},\{ \begin{array}{*{20}{c}} {{x_{pi}} \ge 0}\\ {{x_{pj}} \ge 0} \end{array}$.
18	$x_p^{\text{*}} = {\left( {{x_{p{i^{\text{*}}}}},{x_{p{j^{\text{*}}}}}} \right)^T}$;
19	end
20	end
21	Obtain $x_p^{\text{*}} \equiv x_p^{\text{*}}$ with padded zeros for $t \notin {\tau _c}$.//t is the tth material
22	Output: ${\vec x^{\left( 0 \right)}} = {\text{col}}{\left( {x_1^{\text{*}},x_2^{\text{*}}, \ldots ,x_p^{\text{*}}} \right)^T}$//initial value
	***************Proposed EP method******************
23	Input: ${\vec x^{\left( 0 \right)}} = {\text{col}}{\left( {x_1^{\text{*}},x_2^{\text{*}}, \ldots ,x_p^{\text{*}}} \right)^T}$//initial value
24	For each iteration n = 1,..., Niter //Niter is the total number of iterations
25	Compute Hessian H and gradient $\vec q$ using equations (16) and (17).
26	For each material candidate ${\tau _c} = \left( {i,j} \right) \in {\Omega }$
27	${\vec x_p}\left( {{\tau _c}} \right) = {\left[ {{x_{ip,\,}}{x_{jp\,}}} \right]^T}$, $\vec x_p^{\left( n \right)}\left( {{\tau _c}} \right) = {\left[ {x_{ip}^{\left( n \right)},\,x_{jp}^{\left( n \right)}} \right]^T}$
28	Form $H\left( {{\tau _c}} \right)$ and $\vec q\left( {{\tau _c}} \right)$ from elements in H and $\vec q$ with indexes corresponding to ${\tau _c} = \left( {i,j} \right)$.
29	Find the optimal ${\widehat {\vec x}_p}\left( {{\tau _c}} \right)$ and the function value $\phi _p^{\left( n \right)}\left( {{{\widehat {\vec x}}_p}\left( {{\tau _c}} \right)} \right)$ of equation (21) using GSMO.
30	end
31	$\widehat {{\tau _c}} = _{\,\tau \in {\Omega }}^{\,argmin}\phi _p^{\left( n \right)}\left( {{{\widehat {\vec x}}_p}\left( {{\tau _c}} \right)} \right)$
32	Obtain ${\widehat {\vec x}_p} \equiv {\widehat {\vec x}_p}\left( {\widehat {{\tau _c}}} \right)$ with padded zeros for $t \notin {\tau _c}$.
33	Update ${\vec x^{\left( {n + 1} \right)}} = \left( {{{\widehat {\vec x}}_1}, \ldots ,{{\widehat {\vec x}}_p}, \ldots ,{{\widehat {\vec x}}_P}} \right)$.
34	end
35	Output: ${\vec x^{\left( {Niter} \right)}}$//The decomposed material images

We evaluate the proposed method using a digital phantom and a Catphan©600 phantom. In the digital phantom study, the three material bases bone, muscle and adipose and the mixture of muscle and adipose are inserted into six rods placed in a digital adipose container. The digital phantom is an ellipse with the long axis length of 220 mm and the short axis length of 160 mm. The LACs of basis materials are obtained from the National Institute of Standards and Technology (NIST) database. To evaluate the proposed method using various imaging doses and energies, the digital phantom data is simulated using different number of incident x-ray photons and energies to mimic the scanning protocols of low- and normal- doses using various energies.

The Catphan©600 phantom data is measured from a tabletop cone-beam CT (CBCT) system with the geometry in accordance with that of a Varian On-Board Imager (OBI, www.varian.com). The system utilizes a CB4030 flat-panel detector (Varian Medical Systems, www.varian.com) which is composed of 1024 × 768 pixels with the pixel size of 0.388 mm × 0.388 mm. To suppress the CBCT scatter photons, the longitudinal beam width of the projection is limited to 15 mm.

The clinical patient data which is scheduled for the abdomial scan in our hospital are acquired on the 256-slice Philips iCT scanner (www.usa.philips.com/healthcare/) using the standard multi-phase protocol including a contrast-free scan and three more contrast-enhanced scans. The contrast-enhanced scans are classified as arterial phase, portal venous phase and delayed phase according to the contrast agent traversal time throughout the patient body. The scanning and reconstruction details of the above three data sets are shown in Table 2.

To evaluate the performance of the proposed method, we compare the proposed algorithm with a classical separate low-pass filtration method (Rutherford et al 1976, Johns and Yaffe 1985) which suppress the decomposition noise using the TMA method as the initialization. For a fair comparison, the image noise standard deviation (STD) using the low-pass filtration method is kept as the same as that using the proposed EP method.

3.3.1. Digital and catphan phantom study

We define the volume fraction accuracy (VFA) of decomposition as:

$$\begin{equation}VFA = \left( {1 - \frac{1}{{{T_0}}}\mathop \sum \limits_{t = 1}^{{T_0}} \frac{{\left| {{x_t}^{truth} - {{\bar x}_t}} \right|}}{{{x_t}^{truth}}}} \right) \times 100{\% }\end{equation} \tag{ 27 }$$

where ${x_t}^{truth}$ and ${\bar x_t}$ are the ground truth of decomposed volume fraction and the mean value within the ROI for the tth material, respectively.

The STD is applied to quantitatively measure each decomposed image noise in an ROI. STD is defined as:

$$\begin{equation}STD = \sqrt {\frac{1}{M}\mathop \sum \limits_{m = 1}^M {{({x_{tm}} - {{\bar x}_t})}^2}} \end{equation} \tag{ 28 }$$

where m is the pixel index within the ROI. ${x_{tm}}$ is the mth pixel value in the tth material image. M is the total number of pixels in the selected ROI. The overall decomposition STD is defined as the averaged value of STDs in all the decomposed images.

To further evaluate the noise frequency characteristics of the low-pass filtration method and the proposed EP method at the comparable noise level, we measure the noise power spectrum (NPS) in a uniform area of the decomposed materials. The 2D NPS is defined as

$$\begin{equation}NPS \approx {\left| {DF{T_2}\left\{ f \right\}} \right|^2}\end{equation} \tag{ 29 }$$

where $f$ denotes the ROI in which gray values are offset to achieve a zero mean, $DF{T_2}\left\{ f \right\}$ is the 2D Discrete Fourier Transform (DFT) on $f$ (Xue et al 2017).

In addition to the NPS, we measure the modulation transfer function (MTF) in the digital and Catphan phantom studies to investigate the decomposed image spatial resolution using the low-pass filtration method and the proposed EP method (Richard et al 2012). The line spread function (LSF) is calculated using the gradient of the object edge profile and is transfered to MTF using the Fourier transform. We average the profile of adjacent boundaries to acquire the resultant MTF with the minimal fluctuation due to image noise. The measured frequencies at MTF magnitude decreased to 0.5 (−3dB) are compared to evaluate the relative spatial resolution (Samei et al 1998).

To evaluate the impact of energy and noise on the composition matrix in the digital phantom study, we calculate the condition number of the composition matrix. The condition number measures the stability of the solution to a linear system of equations when a small perturbation is applied to the system. The higher condition number indicates the more instability of the system after the perturbation (Niu et al 2014).

3.3.2. Clinical data study

The patient data is very complicated and no mature commercial MMD method is applied in the clinic. True value does not exist as a reference to evaluate the proposed method. Instead, we apply the proposed SECT MMD algorithm to acquire the decomposed materials and perform virtual non-contrast imaging (VNC) using the decomposed materials to verify the accuracy of the proposed method indirectly. VNC, which is commonly used in the clinic, is a technique to remove the contrast agent from the original contrast-enhanced CT image to reduce the radiation dose. Currently, VNC is mostly achieved using dual-energy CT technology. In this paper, we apply the proposed method for SECT to acquire the decomposed materials and substitute the contrast agent by the blood with the same volume fraction. The VNC image is synthesized by the linear combination of the substituted materials. The VNC image at a given pixel p and energy E is reconstructed as (Mendonca et al 2014)

$$\begin{equation}\mu _p^{VNC}\left( E \right) = {x_{p{t_1}}}\left( E \right){\mu _{{t_2}}}\left( E \right) + \mathop \sum \limits_{\begin{array}{*{20}{c}} {t = 1\,} \\ {t \ne {t_1}} \end{array}}^{{T_0}} {x_{pt}}\left( E \right){\mu _t}\left( E \right)\end{equation} \tag{ 30 }$$

where the t₁-th material refers to the contrast agent and the t₂-th material refers to the blood (Mendonca et al 2014).

For the evaluation of the clinical data, the CT number percentage error and the root-mean-square-relative error (RMSRE (%)) of the CT number are summarized to quantify the accuracy of VNC generated using the proposed method. The CT number percentage error is calculated as:

$$\begin{equation}{({\text{percentage}}\,\,\,{\text{error}})_t} = \frac{{\left| {{{(\mu )}_t} - {{({\mu ^{truth}})}_t}} \right|}}{{{{({\mu ^{truth}})}_t}}}{\,} \times 100\% ,\end{equation} \tag{ 31 }$$

where ${\left( \mu \right)_t}$ is the average CT number of the tth interested tissue or organ in the VNC image and ${({\mu ^{truth}})_t}$ is the true CT number of the tth interested tissue or organ in the contrast-free image. The RMSRE is calculated as:

$$\begin{equation}{\text{RMSRE}} = \sqrt {\frac{1}{{{T_0}}}\mathop \sum \limits_t^{{T_0}} {{\left( {{{\left( {{\text{percentage error}}} \right)}_t}} \right)}^2}} {\,} \times 100{\text{\% }}\end{equation} \tag{ 32 }$$

where ${T_0}$ denotes the number of interested tissues and organs in the abdomen slice of the human.

The digital phantom on the contrast rods slice consists of four materials including bone, muscle, adipose and air. The high- and low-energy CT images are shown in figure 1. The decomposed basis material images are shown in figure 2. The brightness in the decomposed image indicates the volume fraction of each basis material. The proposed TMA method initially decomposes the CT image into four material images in the high- and low-energy studies. Especially in the ROI3, the proposed TMA method successfully differentiates the mixture area into muscle and adipose, indicating the mixture decomposition capability of the TMA method. On the basis of the TMA method, the proposed EP scheme further suppresses the decomposition noise while maintains the material edge, which is indicated as the smoother area of each material and the comparable material edge.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

For quantitative analysis, we select the ROIs indicated by the dashed circles as shown in figure 1. The average values and noise STDs are summarized in table 3. The VFA of the proposed TMA method is 82.49%. The near to one VFA indicates high decomposition accuracy of the proposed TMA method. Especially, in the mixture region of ROI3, the proposed TMA method decomposes the volume fraction of 0.596 and 0.382. The near to the ground truth (i.e. 0.7 and 0.3) volume fraction of the mixture area indicates the mixture decomposition capability of the proposed TMA method. Compared with the TMA method, the proposed EP method further increases the VFA by a factor of 12.23% and decreases the STD from 0.2492 to 0.0401.

The proposed method suppresses the noise into the comparable level with the low-pass filtration method. We measure the NPS of the decomposed adipose image generated using the TMA method, low-pass filtration method and the EP method to evaluate the noise frequency characteristics. The result is shown in figure 3. Compared with the TMA method, the NPS generated using the low-pass filtration method has gone through a great change while the proposed EP method remains more NPS structure.

Zoom In Zoom Out Reset image size

The MTF is calculated to investigate the spatial resolution maintenance of the proposed EP method. We measure the typical MTFs of muscle and mixture using the low-pass filtration and EP methods. The result is plotted in figure 4. The spatial resolution using the proposed EP method is increased by an overall factor of 0.48 when the MTF magnitude is decreased to 50% compared with that using the low-pass filtration method.

Zoom In Zoom Out Reset image size

To evaluate the impact of energy and noise on the composition matrix, we calculate the condition number of the composition matrix. The condition number result is shown in figure 5 as the red parts. The crosses and marks indicate the condition number in low- and normal- dose studies. In both studies, the condition number of the composition matrix increases along with the energy, indicating the worse decomposition stability at a higher energy. Additionally, the condition number curves in the low and normal dose studies are in high consistence. This observation shows the measurement of composition matrix is independent of the imaging dose in the CT scan.

Zoom In Zoom Out Reset image size

The volume fraction accuracy is calculated to evaluate the impact of energy and noise on the decomposition and is shown as blue part in figure 5. The upper, lower, left and right triangles represent the volume fraction accuracy using TMA and EP method in the low- and normal dose studies, respectively. The blue solid, dotted, short dotted, long and short dotted lines are generated to fit the triangle scatters. Both in the low- and normal dose studies, the linearly fitted curves decrease along with the increasement of the energy, indicating a worse decomposition accuracy at higher energy. This observation is consistent with the condition number increasement along with the energy. The higher energy leads to a larger condition number which means the stability of the decomposition is degraded and thus results in inferior decomposition accuracy.

Additionally, the decomposition accuracy declines severely when the noise is magnified in the low dose scan as shown in figure 5. The reason is that the magnified noise leads to the incorrect material composition selection within the two-material candidates looping operation. The incorrect material composition selection ultimately results in the decomposition bias, namely low decomposition accuracy. As illustrated by the higher dotted blue line than the solid blue line in figure 5, the decomposition accuracy decline is alleviated using the propose EP method compared with that using the TMA method in the low-dose study since the applied edge-preserving achieves noise suppression within the material decomposition.

The basis materials selected in the Catphan phantom on the contrast rods slice are Teflon, Delrin, iodine solution, Polymethylpentene (PMP), inner soft-tissue and air. The reconstructed CT image are shown in figure 6. The decomposed basis material images are shown in figure 7. The brightness in the decomposed image indicates the volume fraction of each basis material. The rods inside the Catphan phantom are composed of the pure material and thus the ideal volume fraction (i.e. ground truth) for each rod is 1. The proposed TMA method initially decomposes the Catphan CT image into six material images in the high- and low-energy studies. Compared with the proposed TMA method, the proposed EP scheme further suppresses the decomposition noise while maintains the material edge.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

For quantitative analysis, we select the ROIs indicated by the dashed circles as shown in figure 6. The average values and noise STDs are summarized in table 4. The VFA of the proposed TMA method is 76.78%. Compared with the TMA method, the proposed EP method further increases the VFA by a factor of 18.9% and decreases the STD from 0.1947 to 0.0406.

The proposed method suppresses the noise into the comparable level with the low-pass filtration method. We further measure the NPS of the decomposed soft-tissue image generated using the TMA method, low-pass filtration method and the EP method to evaluate the noise frequency characteristics. The result is shown in figure 8. Compared with the TMA method, the NPS generated using the low-pass filtration method has gone through a great change while the proposed EP method remains more NPS structure.

Zoom In Zoom Out Reset image size

The MTF is calculated to investigate the spatial resolution maintenance of the proposed EP method. We measure the typical MTFs of Delrin and PMP using the low-pass filtration and EP methods in the high- and low-energy scans. The result is plotted in figure 9. The spatial resolution using the proposed EP method is increased by an overall factor of 1.55 when the MTF magnitude is decreased to 50% compared with that using the low-pass filtration method.

Zoom In Zoom Out Reset image size

The CT images of contrast-free phase, arterial phase, portal venous phase and delayed phase are shown in the first row in the figure 10. We apply bone, muscle, contrast agent, adipose and air as the basis materials to perform the material decomposition using the proposed method for the three contrast-enhanced phases. By applying Eq. (30), we acquire the VNC image of arterial, portal venous and delayed phases as shown in figure 10. The VNC images of three phases generated utilizing the proposed method are faithfully consistent with the contrast-free image (i.e. the ground truth image), especially the kidney structure as pointed by the yellow arrow in figure 10. It should be noted that each two phases have a few seconds time interval. The four phases are not completely consistent due to the organ motion in the scanning interval, for instance the liver region pointed by the red arrow in the contrast-free image in figure 10.

Zoom In Zoom Out Reset image size

For quantitative analysis, we select five regions of interest: bone (ROI1), kidney (ROI2), liver (ROI3), muscle (ROI4) and adipose (ROI5) as shown in the arterial VNC image in figure 10 to calculate the CT number percentage error and RMSRE. The quantitative results are shown in table 5. The RMSREs of arterial, portal venous and delayed phases are 3.6%, 3.0% and 2.2%, respectively. The VNC images generated by the proposed method achieve high accuracy using the contrast-free CT image as the ground truth. The proposed method successfully removes the contrast agent from the contrast-enhanced images and maintains consistent CT number to the contrast-free CT image.

The proposed EP method contains two tunable parameters including ${\beta _t}$ which controls the noise-resolution and ${\delta _t}$ which controls the edge preservation. ${\beta _t}$ and ${\delta _t}$ are material-specific since each material has its unique volume fraction distribution and neighbored pixels. We thus empirically select the optimal combination of ${\beta _t}$ and ${\delta _t}$ to achieve the tradeoff between the noise suppression and spatial resolution maintenance. The regularization coefficients ${\beta _t}$ and the edge-preserving parameters ${\delta _t}$ for each material in the digital and Catphan studies are listed in table 6.

Table 6. The regularization coefficients and edge-preserving parameters.

Data	${\beta _t}$		${\delta _t}$
Digital phantom	0.0010	0.0100	1.0000	1.0000
(for bone, muscle, adipose and air)	0.0100	0.1000	1.0000	0.1000
Catphan phantom	0.0600	0.0500	0.2500	0.2500
(for teflon, delrin, iodine, PMP, soft-tissue and air)	0.0250	0.0100	1.0000	0.6000
	0.0100	0.0700	0.6000	0.6000