Using edge-preserving algorithm with non-local mean for significantly improved image-domain material decomposition in dual-energy CT

Suppose the total number of pixels of one CT image is N, for low- and high-energy CT images and material-specific images, equation (1) can be rewritten in matrix form,

$$\begin{eqnarray}&&\vec{\mu}=A\vec{x}.\end{eqnarray} \tag{ 2 }$$

Here A is a $2N\times 2N$ material decomposition matrix and it can be derived from equation (1) as,

$$\begin{eqnarray}A=\left(\begin{array}{c} {{\mu}_{1H}}I & {{\mu}_{2H}}I \\ {{\mu}_{1L}}I & {{\mu}_{2L}}I \end{array}\right),\end{eqnarray} \tag{ 3 }$$

with I the $N\times N$ identity matrix. $\vec{\mu}$ and $\vec{x}$ are column vectors with dimension of 2N as they consist of high- and low-energy CT image, and decomposed material-specific images with column vector form, respectively. Namely,

$$\begin{eqnarray}&&\vec{\mu}=\left(\begin{array}{c} {{{\vec{\mu}}}_{H}} \\ {{{\vec{\mu}}}_{L}} \end{array}\right),~\vec{x}=\left(\begin{array}{c} {{{\vec{x}}}_{1}} \\ {{{\vec{x}}}_{2}} \end{array}\right),\end{eqnarray} \tag{ 4 }$$

with ${{\vec{\mu}}_{H}}$ and ${{\vec{\mu}}_{L}}$ the measured high- and low-energy CT images, respectively, and ${{\vec{x}}_{1}}$ and ${{\vec{x}}_{2}}$ the two material-specific images. Note that both ${{\vec{\mu}}_{H,L}}$ and ${{\vec{x}}_{1,2}}$ are represented as vectors. To obtain material-specific images from high- and low-energy CT images, one can directly use matrix inversion to solve equation (1), i.e.

$$\begin{eqnarray}&&\vec{x}={{A}^{-1}}\vec{\mu},\end{eqnarray} \tag{ 5 }$$

where

$$\begin{eqnarray}{{A}^{-1}}=C\left(\begin{array}{c} -{{\mu}_{2L}}I & {{\mu}_{2H}}I \\ {{\mu}_{1L}}I & -{{\mu}_{1H}}I \end{array}\right),\end{eqnarray} \tag{ 6 }$$

with matrix determinant $C={{\left({{\mu}_{2H}}{{\mu}_{1L}}-{{\mu}_{1H}}{{\mu}_{2L}}\right)}^{-1}}$. However, material decomposition in this manner can yield significantly increased noise and severely degraded noise-to-signals (NSRs) for the material-specific images. A method taking some measures to reduce the noise and to recover NSRs is highly desirable.

To reduce the amplified image noise in the decomposed material-specific images while keeping the images as accurate as possible (both quantitative measurement and spatial resolution), we developed an image-domain material images denoising algorithm. The algorithm was based on the HYPR-LR (local highlY constrained backPRojection reconstuction) framework (Mistretta et al 2006, Johnson et al 2008, Leng et al 2011), which was first proposed for reconstruction of time-resolved magnetic resonance imaging (MRI) using highly undersampled projection. The HYPR-LR framework has been successfully applied to a broad range of medical imaging applications, such as dynamic MRI (Johnson et al 2008), CT angiography (Supanich et al 2009), dynamic positron emission tomography (Christian et al 2010), spectral CT (Leng et al 2011) and myocardial perfusion imaging (Speidel et al 2013). For noise reduction in spectral CT, the HYPR-LR algorithm treated CT images at different energies as four-dimensional CT images, i.e. the energy dimension was regarded as time dimension. It started with producing a composite image ${{\vec{\mu}}_{c}}$ by averaging the energy images with equal weighting, i.e. yielding an image with lower image noise level by using the images from different energy bins of photon-counting detector CT, or from high- and low-energy images of dual-energy CT. A weight image was then generated by calculating the ratio of two filtered images, which were obtained by filtering both the energy images ${{\vec{\mu}}_{e}}$ and the composite image ${{\vec{\mu}}_{c}}$. The final image ${{\vec{\mu}}_{he}}$ was obtained by multiplying the weight image and the composite image. The mathematical form of HYPR-LR algorithm is as follows,

$$\begin{eqnarray}&&{{\vec{\mu}}_{he}}=\frac{{{{\vec{\mu}}}_{e}}\otimes K}{{{{\vec{\mu}}}_{c}}\otimes K}\centerdot {{\vec{\mu}}_{c}}\end{eqnarray} \tag{ 7 }$$

with ${{\vec{\mu}}_{e}}$ the image with different energy, and K the low-pass filter kernel (usually an uniform kernel). The symbol $\otimes$ stands for convolution operation. For dual-energy material decomposition, HYPR-LR has been applied to the high- and low-energy CT images to yield noised reduced images, as well as superior material-specific basis images (Leng et al 2011).

In this study, we demonstrated the HYPR-LR algorithm can be applied directly to the dual-energy material-specific images with composite image ${{\vec{\mu}}_{c}}$. For dual-energy high- and low-energy images ${{\vec{\mu}}_{H,L}}$, the processed individual energy images ${{\vec{\mu}}_{he,H}}$ and ${{\vec{\mu}}_{he,L}}$ using the HYPR-LR algorithm are calculated as follows,

$$\begin{eqnarray}&&{{\vec{\mu}}_{he,H}}=\frac{{{{\vec{\mu}}}_{H}}\otimes K}{{{{\vec{\mu}}}_{c}}\otimes K}\centerdot {{\vec{\mu}}_{c}},\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&{{\vec{\mu}}_{he,L}}=\frac{{{{\vec{\mu}}}_{L}}\otimes K}{{{{\vec{\mu}}}_{c}}\otimes K}\centerdot {{\vec{\mu}}_{c}}.\end{eqnarray} \tag{ 9 }$$

Based on equation (5), the basis material image ${{\vec{x}}_{h1}}$ calculated using ${{\vec{\mu}}_{he,H}}$ and ${{\vec{\mu}}_{he,L}}$ is expressed as,

$$\begin{eqnarray}&&{{\vec{x}}_{h1}}=-C\left({{\mu}_{2L}}{{\vec{\mu}}_{he,H}}-{{\mu}_{2H}}{{\vec{\mu}}_{he,L}}\right)\\ &&=-C\left({{\mu}_{2L}}\frac{{{{\vec{\mu}}}_{H}}\otimes K}{{{{\vec{\mu}}}_{c}}\otimes K}\centerdot {{\vec{\mu}}_{c}}-{{\mu}_{2H}}\frac{{{{\vec{\mu}}}_{L}}\otimes K}{{{{\vec{\mu}}}_{c}}\otimes K}\centerdot {{\vec{\mu}}_{c}}\right)\\ &&=-C\left(\frac{\left({{\mu}_{2L}}{{{\vec{\mu}}}_{H}}-{{\mu}_{2H}}{{{\vec{\mu}}}_{L}}\right)\otimes K}{{{{\vec{\mu}}}_{c}}\otimes K}\centerdot {{\vec{\mu}}_{c}}\right)\\ &&=\frac{{{{\vec{x}}}_{1}}\otimes K}{{{{\vec{\mu}}}_{c}}\otimes K}\centerdot {{\vec{\mu}}_{c}}.\end{eqnarray} \tag{ 10 }$$

The above derivation has used the algebraic properties of convolution, i.e. the distributivity and the associativity of convolution with scalar multiplication. Following the fashion in equation (10), we have the same expression for the other basis material image ${{\vec{x}}_{h2}}$. Equation (10) demonstrates the HYPR-LR algorithm can be applied directly to the material-specific image. Thus we have,

$$\begin{eqnarray}&&{{\vec{x}}_{hb}}=\frac{{{{\vec{x}}}_{b}}\otimes K}{{{{\vec{\mu}}}_{c}}\otimes K}\centerdot {{\vec{\mu}}_{c}}.\end{eqnarray} \tag{ 11 }$$

Here ${{\vec{x}}_{b}}$ is the material-specific basis image, i.e. ${{\vec{x}}_{1}}$ or ${{\vec{x}}_{2}}$, and ${{\vec{x}}_{hb}}$ is the processed basis image. In order to further reduce decomposed images noise without loss of high-frequency edge information, an edge-preserving non-local mean (NLM) (Buades et al 2005) was introduced into the HYPR framework. Edge-preserving techniques were usually applied to yield noise reduced, spatial resolution well preserved images (Chen et al 2008, Wang et al 2009, Zhu et al 2009, Chen et al 2012, Chun et al 2014). In this study, NLM was employed to generate the weight image for the HYPR framework and the resulting algorithm was referred as HYPR-NLM, hence equation (11) was rewritten in pixel-wise fashion as,

$$\begin{eqnarray}&&{{x}_{hb}}(i)=\frac{\underset{j\in {{\Omega }_{i}}}{\mathop \sum}\,\omega (i,j){{x}_{b}}(\,j)}{\underset{j\in {{\Omega }_{i}}}{\mathop \sum}\,\omega (i,j){{\mu}_{c}}(\,j)}{{\mu}_{c}}(i),\end{eqnarray} \tag{ 12 }$$

where the weight $\omega (i,j)$ depicts the similarity between the pixels i and j, and it satisfies the constraint conditions $0\leqslant \omega (i,j)\leqslant 1$ and $\sum_{j}\omega (i,j)=1$ (as demonstrated later). The pixel dependent summation domain ${{\Omega }_{i}}$ denotes a search-window centered at the pixel i and it is usually a square neighborhood with fixed size. Thus each pixel of the filtered image is a weighted summation of a square neighborhood. For the right hand side of equation (12), $\omega (i,j)$ is applied to the material images in the numerator, and to the composite image in the denominator.

The similarity between pixels i and j can be measured using a weighted Euclidean distance of two square neighborhoods centered at pixels i and j, thus the weight $\omega (i,j)$ is defined as,

$$\begin{eqnarray}&&\omega (i,j)=\frac{1}{Z(i)}\exp \left(-\frac{\parallel \vec{\mu}\left({{\Theta}_{i}}\right)-\vec{\mu}\left({{\Theta}_{j}}\right)\parallel _{2,a}^{2}}{{{d}^{2}}}\right),\end{eqnarray} \tag{ 13 }$$

where ${{\Theta}_{i}}$ and ${{\Theta}_{j}}$ are two square neighborhoods centered at pixels i and j, respectively. The parameter d controls the decay of the exponential function and it acts as a degree of filtering. The parameter a is the standard deviation (SD) of a Gaussian kernel which has the same size as the square neighborhood. During numerical implementation, the Gaussian kernel gives weights to the Euclidean distance of the two square neighborhoods. Z(i) is the normalization constant calculated as,

$$\begin{eqnarray}&&Z(i)=\underset{j}{\mathop \sum}\,\exp \left(-\frac{\parallel \vec{\mu}\left({{\Theta}_{i}}\right)-\vec{\mu}\left({{\Theta}_{j}}\right)\parallel _{2,a}^{2}}{{{d}^{2}}}\right).\end{eqnarray} \tag{ 14 }$$

Based on the above definition, $\omega (i,j)$ satisfies the constraint conditions.

Since dual-energy material decomposition can yield noise amplified images, for comparison, we also used an iterative image-domain material decomposition method which significantly reduced the noise of the material-specific images with superior performance on image spatial resolution and low-contrast detectability (Niu et al 2014). The method balanced the data fidelity of the value of the material image and a quadratic penalty using an optimization framework, and the unconstrained optimization problem was solved by the nonlinear conjugate gradient method (see the appendix for details). In this work, we referred this method as Iter-DECT.

To evaluate the proposed HYPR-NLM dual-energy material decomposition algorithm, we first used simulated phantom data. For all of the simulation studies, the high- and low-energy CT images were first employed to generate the material-specific basis images, which were then processed using HYPR-NLM method to yield noise reduced images. These images were compared to the images generated using direct matrix inversion, Iter-DECT and HYPR-LR methods and the results were quantitatively analyzed.

In all of the simulations, a 2D fan-beam CT geometry was performed. The distance from the source to the center of rotation was 785 mm and the distance from the source to the detector was 1200 mm. A circular scan was simulated and a total of 720 projections per rotation were acquired in an angular range of ${{360}^{{}^\circ}}$. The detector pixel size was 0.388 mm and the detector had 1024 pixels. The dual-energy spectra were 100 kVp and 140 kVp, which were generated using the SpekCalc software (Poludniowski et al 2009) with 12 mm Al and 0.4 mm Sn  +  12 mm Al filtration, respectively. The phantom was a water cylinder with inserts that contains a series of six solutions of varying iodine concentrations (range, 0–20 mg ml⁻¹). The diameters of the water cylinder and the inserts were 198 mm and 22.5 mm, respectively.

To be realistic, Poisson noise was considered during the numerical simulations. The fan beam CT projection data was created by polychromatic forward projecting the numerical phantom. Specifically, after introducing Poisson noise (Nuyts et al 2013), the projection data can be represented as:

$$\begin{eqnarray}&&I=\text{Poisson}\left\{N\int_{0}^{{{E}_{\text{max}}}}\text{d}E\,\Omega (E)\,\eta (E)\,\text{exp}\left[-\int_{0}^{l}\mu (E,s)\text{d}s\right]\right\},\end{eqnarray} \tag{ 15 }$$

with N the total number of photons and was set to $3\times {{10}^{5}}$ and $1.5\times {{10}^{5}}$ per ray for the low- and high-energy CT scans, respectively. $\eta (E)$ was the energy dependent response of the detector, which was considered to be proportional to photon energy E for energy-integrating detectors. ${{E}_{\text{max}}}$ was the maximum photon energy of the polychromatic spectrum $\Omega (E)$. $\mu (E,s)$ was the energy-dependent linear attenuation coefficient and was obtained from the National Institute of Standards and Technology (NIST) database. l was the propagation path length for each ray and can be calculated using either analytical methods or numerical methods. The x-ray spectra and the phantom are shown in figure 1. The linear attenuation coefficients of the iodine concentrate inserts calculated using the mixture rule, are show in figure 2. Note that the y axis is plotted on a logarithmic scale. First order beam hardening correction (water correction) was performed by simply mapping the polychromatic raysum to the corresponding monochromatic value with the incorporation of the 100 kVp and 140 kVp energy spectra.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The calculated iodine concentrations in the phantom using different algorithms were compared with known true iodine concentrations. For comparison, direct matrix inversion without noise was also performed and the results were regarded as absolute truth. Difference image between the absolute truth and the images processed using different kinds of algorithms were also generated to emphasize resolution preservation and noise reduction of these algorithms.

2.4.1. Parameters selection.

2.4.2. Iteration formulation.

2.4.3. Effect of dose level.

To evaluate the proposed method, two retrospective clinical patient studies were also performed with a dual-source DECT scanner (Somatom Definition, Siemens Healthcare). The system acquired high- and low-energy data with two x-ray tubes with corresponding detector rings mounted onto a rotating gantry with angular offset of ${{90}^{{}^\circ}}$. The two tubes (tube 1 and 2) operated independently with regard to tube voltage, tube current as well as tube filtration. Dual-energy data were acquired using the following parameters: tube 1, 140 kVp, 287 mA and tube 2, 100 kVp, 412 mA for patient 1; tube 1, 140 kVp, 345 mA and tube 2, 100 kVp, 500 mA for patient 2. The 140 kV spectrum was filtered with a tin filter. Images were reconstructed using the commercial software with the convolution kernel B25f for patient 1 and B26f for patient 2. All of the image slices were reconstructed in 0.75 mm slice thickness. The matrix size for each slice was $512\times 512$, and the pixel size was $0.38\times 0.38$ mm². Blending images were also generated with a linear combination (with weighting factors of 0.3 and 0.7) of the 100 kV and 140 kV images. To obtain contrast-enhanced CT images, the contrast agent of 300 mg I ml⁻¹ was injected in an antecubital vein (300 mg I ml⁻¹, Ultravist 300, Bayer HealthCare).

The noise in the raw CT images was quantified using region-of-interest (ROI) analysis. Two ROIs were selected to calculate noise-to-signal ratio (NSR) of myocardium and ventricle, respectively. Myocardial perfusion imaging was conducted and iodine concentration map was calculated using direct matrix inversion, HYPR-LR, Iter-DECT and HYPR-NLM. We conducted all of the evaluations on a personal laptop which featured four Intel Core i7-4700HQ CPUs, and all of the material decomposition algorithms were implemented using MATLAB (The MathWorks, Inc., Natick, MA, United States).

Image reconstruction results of the numerical iodine concentrate phantom are shown in figure 3. Water and the 20 mg ml⁻¹ iodine concentrate are chosen as the basis materials for dual-energy material decomposition. Other iodine concentrates are expected to show up in the decomposed image corresponding to the contribution of their mass attenuation coefficients to the mass attenuation coefficient of the chosen basis material. In the first row, the 100 kV and 140 kV images and the difference between these two images are depicted. The 100 kV image shows superior iodine contrast, as expected. Note that the images of both energies have been corrected using the first order beam hardening correction algorithm and there are residual high order beam hardening artifacts (streaks) between high attenuation iodine concentrates, especially for the 100 kV image where the spectrum is much softer.

Zoom In Zoom Out Reset image size

The second and the third rows show dual-energy material decomposition results using direct matrix inversion with and without water correction. As can be seen, dual-energy material decomposition yields residual cupping artifacts in the material images without water correction, and the water image and the iodine image are inhomogeneous. The difference images further clearly show the cupping artifacts and inhomogeneity. This is because the image-domain matrix inversion just linearly combine the high- and low-energy images, and does not take the nonlinear attenuation process of x-ray projection into account.

Figure 4 shows the noise of the ROIs (labelled in figure 1(a)) of the iodine images processed using HYPR-NLM method with different parameters. The noise reduction of HYPR-NLM with different search window sizes are depicted in figure 4(a). For comparison, the results of direct matrix inversion and Iter-DECT are also present. For non-iterative HYPR-NLM formulation (i.e. the iteration number N  =  1), noise magnitude is significantly reduced using HYPR-NLM method and there is no significant change between HYPR-NLM with different search windows sizes, namely, with search window size from $11\times 11$ to $19\times 19$. Thus in the following studies, we use search window size of $11\times 11$ for computational consideration. The results of noise reduction using iterative HYPR-NLM with different iterations are shown in figure 4(b). As can been seen, with two iterations (N  =  2), HYPR-NLM provides results that have comparable noise magnitude as the results of Iter-DECT method.

Zoom In Zoom Out Reset image size

Figure 5(a) shows the results of the iodine images processed using HYPR-LR and HYPR-NLM with different window sizes. Different from HYPR-LR method, which results in spatial blurring as the uniform kernel size increases from $5\times 5$ to $11\times 11$, the edges do not have significant change for the HYPR-NLM with search window sizes from $11\times 11$ to $19\times 19$ and iteration numbers from one to two. As indicated in the difference images, there are clear edge effects for the HYPR-LR method, while there are marginal edge effects for the HYPR-NLM method, suggesting the edge of the iodine images are well preserved using the HYPR-NLM method with different search window sizes. Line profiles of the iodine images processed using HYPR-LR with different window sizes, and HYPR-NLM with different window sizes as well as iteration numbers are illustrated in figures 5(b) and (c). These profiles clearly demonstrate HYPR-NLM can provide superior images with respect to edge-preservation, and the method is robust against parameters selection.

Zoom In Zoom Out Reset image size

Figure 6 shows the results of dual-energy material decomposition images using HYPR-NLM method at different number of photons and quantitative analysis of the ROIs. Again, the difference image is the subtraction of the image obtained using noise-free direct matrix inversion and the HYPR-NLM image. As can be seen, HYPR-NLM works well for different scenarios. ROIs analyses indicate HYPR-NLM significantly reduces noise while preserving quantitative iodine concentration. The standard deviation of the ROI reduces as the number of photons increases. Note that H and L stand for the number of high- and low-energy photons, respectively.

Zoom In Zoom Out Reset image size

Material-specific images of the numerical iodine concentrate phantom obtained using direct matrix inversion, HYPR-LR, Iter-DECT and HYPR-NLM algorithms are shown in figure 7. As can be seen, both HYPR-LR and Iter-DECT yield residual edge effects in the difference images, suggesting there are certain degree of spatial resolution degradation. However, there are marginal edge effects for the results of HYPR-NLM algorithm, showing there is minimal loss of spatial resolution. Iodine concentrations of the ROIs (labeled in figure 1(a)) measured using different algorithms are shown in table 1. Compared to the true values of the iodine concentrate inserts, the quantitative measurements of iodine concentration using the four algorithms are well preserved.

Zoom In Zoom Out Reset image size

Figure 8 shows the low-energy (100 kV), high-energy (140 kV) CT images and their linear blending images with weighting factors of 0.3 and 0.7 for the myocardial perfusion imaging of patient 1. Two regions-of-interests (ROI) A and B are labelled respectively on myocardium and ventricle to calculate their NSRs. The NSRs of the myocardium for the 100 kV and 140 kV CT images are 34.6% and 27.9%, respectively. While the NSRs of the ventricle for the 100 kV and 140 kV CT images are 5.3% and $5\%$, respectively. For both ROIs, the NSRs of the 100 kV images are higher than the 140 kV images. This may be the reason why the 140 kV image has a higher weighting factor in the blending image. Water image, and iodine contrast image decomposed using the high- and low-energy images with different algorithms are present in figure 9. As can be seen, dual-energy material decomposition using direct matrix inversion yields increased noise. HYPR-LR reduces the noise level to a certain extent. Iter-DECT significantly reduces the amplified noise. For both of the water and iodine images, HYPR-NLM shows superior image quality. Quantitative measurements of iodine concentration and noise reduction of the labelled ROIs are depicted in table 2.

Table 2. The mean values and standard deviations of the ROIs (labeled as #1 and #2 in figures 9 and 11) of the myocardial data.

		Direct inversion	HYPR-LR	Iter-DECT	HYPR-NLM
Patient 1	Water image	$1.16\pm 0.36$	$1.16\pm 0.16$	$1.16\pm 0.12$	$1.16\pm 0.08$
	Iodine image	$0.99\pm 0.31$	$0.99\pm 0.14$	$0.99\pm 0.12$	$0.99\pm 0.08$
Patient 2	Water image	$1.04\pm 0.34$	$1.04\pm 0.18$	$1.04\pm 0.08$	$1.04\pm 0.07$
	Iodine image	$0.97\pm 0.28$	$0.98\pm 0.14$	$0.98\pm 0.06$	$0.98\pm 0.05$

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The images of patient 2, a case with motion artifact presented in figures 10 and 11, provides a challenging situation for dual-energy material decomposition. This is because myocardial imaging is often compromised by motion blur, leading to DECT artefacts. The blurring edge labelled in the 100 kV image of figure 10 does not have a consistent profile compared to the 140 kV image, which shows a sharper edge here. This happens for the labelled edge for the 140 kV image as well. Again, the blending image shows the linear combination of the 100 kV and 140 kV CT images.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Water images and iodine contrast images, as shown in figure 11, clearly depict the error induced by motion. Still, there are marginal edge effects for HYPR-NLM processed images, suggesting improvements compared to HYPR-LR and Iter-DECT. Dual-energy decomposition using HYPR-NLM leads to more homogeneous results and more significant noise reduction.

Table 2 shows the mean values and standard deviations of the ROIs measured using direct matrix inversion, HYPR-LR, Iter-DECT and HYPR-NLM algorithms. Note that the NSR for direct matrix inversion is approximately $30\%$ for both basis material images and is reduced by factors for approximately 2, 3 and 5 for HYPR-LR, Iter-DECT and HYPR-NLM, respectively. The computational time for each slice per iteration for HYPR-NLM is about 280 s, compared to 300 s for Iter-DECT. There are no significant differences between the mean values measured using the four algorithms, suggesting the quantitative measurements of iodine concentration are acceptable for all four approaches. For noise reduction, direct inversion and HYPR-NLM yield the highest and the lowest noise level, respectively. HYPR-NLM outperforms HYPR-LR in both resolution preservation and noise reduction.