Statistical image-based material decomposition for triple-energy computed tomography using total variation regularization

Abstract

BACKGROUND:

Triple-energy computed tomography (TECT) can obtain x-ray attenuation measurements at three energy spectra, thereby allowing identification of different material compositions with same or very similar attenuation coefficients. This ability is known as material decomposition, which can decompose TECT images into different basis material image. However, the basis material image would be severely degraded when material decomposition is directly performed on the noisy TECT measurements using a matrix inversion method.

OBJECTIVE:

To achieve high quality basis material image, we present a statistical image-based material decomposition method for TECT, which uses the penalized weighted least-squares (PWLS) criteria with total variation (TV) regularization (PWLS-TV).

METHODS:

The weighted least-squares term involves the noise statistical properties of the material decomposition process, and the TV regularization penalizes differences between local neighboring pixels in a decomposed image, thereby contributing to improving the quality of the basis material image. Subsequently, an alternating optimization method is used to minimize the objective function.

RESULTS:

The performance of PWLS-TV is quantitatively evaluated using digital and mouse thorax phantoms. The experimental results show that PWLS-TV material decomposition method can greatly improve the quality of decomposed basis material image compared to the quality of images obtained using the competing methods in terms of suppressing noise and preserving edge and fine structure details.

CONCLUSIONS:

The PWLS-TV method can simultaneously perform noise reduction and material decomposition in one iterative step, and it results in a considerable improvement of basis material image quality.

Keywords

Triple-energy CT image-based material decomposition penalized weighted least-squares total variation

1 Introduction

X-ray computed tomography (CT) plays an essential role in diagnostic radiology and image-guided procedures by providing dependable insight into anatomical structures of the human body. However, conventional CT cannot identify materials or provide quantitative tissue information of the examined objects [1 –5]. Tissues with different material compositions can show same or very similar intensity in a reconstructed CT image. To overcome this limitation, dual-energy CT (DECT) was developed as an alternative solution [6 –11]. Compared to conventional CT, DECT can obtain x-ray attenuation measurements with two energy spectra, which can be used to decompose the examined object into different basis material images [12 –14]. As a simple technical approach, DECT has achieved considerable success in advanced CT applications, including in automated bone removal in CT angiography, ischemia areas detection in perfusion CT imaging, urinary stone characterization, atherosclerotic plaque removal, and virtual monoenergetic imaging [15]. However, the measurements obtained from DECT are comprised only of two energy spectra information on the examined object, which limits the discrimination capability of the DECT to only two materials without consideration of additional constraints. Therefore, adding a third energy spectrum in the DECT imaging (i.e., three energy spectra) may help further improve the accuracy of material decomposition. Triple-energy computed tomography (TECT) technique has been employed in many applications wherein the discrimination of more than two materials is necessary [16 –19]. Using measurements acquired at three different energy spectra, TECT can generate three basis material images that provides the material composition information of the examined object. However, given the noisy TECT measurements, the basis material image quality would be severely degraded when material decomposition is performed directly using matrix inversion method [16 , 20–25]. To achieve accurate material discrimination, it is imperative to develop superior material decomposition methods for TECT.

Material decomposition in TECT can be divided into direct decomposition, projection-based, and image-based methods. Direct decomposition reconstructs basis material images from measured TECT projections using the given x-ray spectrum [13 , 26–28]. This method has the advantage that the physical model of an x-ray transmission is treated exactly, and therefore, prior information of the basis material image can be incorporated into the decomposition process. However, the computational load of direct decomposition methods would be a challenge for practical application because of the repeated forward and backward projection operations in the decomposition process [21, 27]. A projection-based method first decomposes the TECT measurements into sinograms of the basis material images, and then, it performs image reconstruction for each of the material-specific sinogram. This methods can effectively avoid beam-hardening artifacts in basis material images as sinograms of the basis material images are obtained before image reconstruction [29, 30]. The main challenge to projection-based methods is the need to calibrate the model of x-ray transmission, which is highly nonlinear and computationally expensive. Image-based methods decompose TECT images into basis material images using the linear approximation of the material decomposition process [11 , 31–34]. This method is practical and is performed in the reconstruction of each projection followed by a material decomposition step. Compared to projection-based methods, which can suffer from data inconsistency, image-based methods can be readily applied to clinical settings because they are directly performed on images obtained from commercial CT scanners. The basis material images can be obtained simply by direct matrix inversion decomposition; however, this method is sensitive to the noise in CT image because of the overlap between the different energy spectra [20, 35]. Thus, the image quality of basis material will be severely degraded, and this can compromise the material discrimination ability of TECT.

In this study, we present a new strategy for noise suppression in image-based material decomposition for TECT using a penalized weighted least-squares (PWLS) criteria with total variation (TV) regularization (PWLS-TV). We analyzed the noise statistical properties of the image-based material decomposition process, where the inverse of the estimated noise variance of the TECT images serve as the weight in the weighted least-squares (WLS) term. TV regularization enforces the smoothness of the basis material image by penalizing the differences of local neighboring pixels. We proposed an alternating optimization algorithm for minimize the objective function of the PWLS-TV decomposition, and compared the PWLS-TV method with direct matrix inversion method [36], TV filtration for the direct matrix inversion result [20], and the mean squared difference method [16]. The novelty of the proposed PWLS-TV material decomposition method is threefold. First, we incorporated the noise statistical properties of the TECT images into material decomposition process to improve the quality of basis material image. Second, PWLS-TV method enables noise suppression and material decomposition simultaneously in one iterative step. Third, an efficient alternating optimization algorithm was developed to minimize solve the objective function of the PWLS-TV method.

2 Methods and materials

2.1 Image-based basis material decomposition for TECT

In image-based basis material decomposition, the linear attenuation coefficient (LAC) of each pixel in a CT image can be formulized as a linear combination of the basis material images [16, 37]. The image-based basis material decomposition for TECT can be written as $(\begin{matrix} f_{E 1} \\ f_{E 2} \\ f_{E 3} \end{matrix}) = (\begin{matrix} μ_{1, E 1} & μ_{2, E 1} & μ_{3, E 1} \\ μ_{1, E 2} & μ_{2, E 2} & μ_{3, E 2} \\ μ_{1, E 3} & μ_{2, E 3} & μ_{3, E 3} \end{matrix}) (\begin{matrix} x_{1} \\ x_{2} \\ x_{3} \end{matrix}),$ (1) where f_E1 = (f_E1,1, f_E1,2, …, f_E1,N), f_E2 = (f_E2,1, f_E2,2, …, f_E2,N), and f_E3 = (f_E3,1, f_E3,2, … , f_E3,N) are the reconstructed TECT images at mean x-ray energy E1, E2, and E3, respectively. x₁ = (x_1,1, x_1,2, …, x_1,N), x₂ = (x_2,1, x_2,2, …, x_2,N), and x₃ = (x_3,1, x_3,2, …, x_3,N) represent the three basis material images (or the decomposed images). N represents the number of pixels in basis material image, and μ_i,j denotes the LAC of basis material i (i = 1, 2, 3) under the energy j (j = E₁, E₂, E₃). We define F = [f_E1, f_E2, f_E3] ^T and X = [x₁, x₂, x₃] ^T, where T denotes the transpose operator. Using the above notations, Equation (1) can be formularized as $F = AX,$ (2) where A represents the material composition matrix defined as $A = (\begin{matrix} μ_{1, E 1} I & μ_{2, E 1} I & μ_{3, E 1} I \\ μ_{1, E 2} I & μ_{2, E 2} I & μ_{3, E 2} I \\ μ_{1, E 3} I & μ_{2, E 3} I & μ_{3, E 3} I \end{matrix}) .$ (3)

In Equation (3), I denotes the identity matrix with size N × N, and the volume fraction X satisfies the following constraint $0 \leq x_{i, j} \leq 1, \forall i, j .$ (4)

Based on Equation (2), the basis material image can be calculated by direct matrix inversion as $X = A^{- 1} F,$ (5) where A^-1 is the material decomposition matrix. Direct matrix inversion material decomposition using Equation (5) is affected by magnified noise in the basis material images because of the large condition number of the material composition matrix A. Therefore, direct matrix inversion material decomposition method is very sensitive to the noise of the TECT images.

2.2 Noise model in image-based material decomposition

We assume that the CT images are independent Gaussian random variables [20] $f_{Ei, k} \tilde{} N ({(a_{Ei})}^{T} b_{k}, var (f_{Ei, k})), i = 1, 2, 3 .$ (6) where f_Ei,k denotes the k-th pixel value at energy Ei, a_Ei = (μ_1,Ei, μ_2,Ei, μ_3,Ei) ^T, b_k = (x_1,k, x_2,k, x_3,k) ^T, and var (f_Ei,k) denotes the variance. The probability density function is $p (f_{Ei, k}; b_{k}) = \frac{1}{\sqrt{2 π var (f_{Ei, k})}} exp (- \frac{{(f_{Ei, k} - {(a_{Ei})}^{T} b_{k})}^{2}}{2 var (f_{Ei, k})}) .$ (7)

The negative log-likelihood function of independent f_Ei,k can be written as $L (X) = - \sum_{i = 1}^{3} \sum_{k = 1}^{N} log (p (f_{Ei, k}; b_{k})) \equiv {(F - AX)}^{T} Σ^{- 1} (F - AX),$ (8) where the symbol ≡ implies that the equation holds within irrelevant constants independent of X. The weight Σ = diag (var (f_E1,1) , …, var (f_E1,N) , var (f_E2,1) , …, var (f_E2,N) , var (f_E3,1) , …, var (f_E3,N)) is a 3N × 3N diagonal matrix, where var (f_E1,k), var (f_E2,k), and var (f_E3,k) are the variance of the k-th pixel in the TECT images at energy E1, E2, and E3, respectively. In this work, we estimate the pixel-wise noise variance by selecting a region of interest (ROI) with size of 3 × 3 that comprises the homogeneous material on the images; then, we calculate the numerical variance in this ROI as the element of the diagonal matrix ∑, which is a method used by other researcher as well [20, 21]. The observations on the noise properties of material decomposition lead to a well-known penalized weighted least-squares cost function $min_{X} {(F - AX)}^{T} Σ^{- 1} (F - AX) + β R (X),$ (9) where R (X) denotes the regularization/penalty term, and β > 0 is a penalty parameter that balances fidelity and regularization terms. The regularization term R (X) is to strengthen the smoothness of the basis material images while preserving the boundary sharpness. In this work, we choose the TV regularization to penalize the differences between one pixel and its local neighboring pixels, which has been widely used in image reconstruction for denoising of CT images while retaining edges [38 –41].

2.3 PWLS-TV for image-based basis material decomposition

Based on the aforementioned noise model, we propose a PWLS-TV basis material decomposition method. The cost function of the PWLS-TV basis material decomposition is written as $min_{X} {(F - AX)}^{T} W^{- 1} (F - AX) + β_{2} TV (X),$ (10) where $W = \frac{1}{β_{1}} A A^{T} + Σ$ , TV (X) is the total variation of X, β₁ and β₂ are positive regularization parameters. The main difference between Equations (9) and (10) is that a generalized WLS term is used in the new cost function. Moreover, it is obvious that when β₁ is sufficiently large, the generalized WLS term in Equation (10) is nearly the same as that in Equation (9). Furthermore, we note that $min_{U} {(F - AU)}^{T} Σ^{- 1} (F - AU) + β_{1} {∥ U - X ∥}_{2}^{2} = {(F - AX)}^{T} W^{- 1} (F - AX),$ (11) where U is an intermediate variable. Therefore, the cost function of the PWLS-TV basis material decomposition in Equation (10) can be rewritten as $\begin{matrix} min_{X} {min_{U} {(F - AU)}^{T} Σ^{- 1} (F - AU) + β_{1} {∥ U - X ∥}_{2}^{2}} + β_{2} TV (X) \\ = min_{X, U} {(F - AU)}^{T} Σ^{- 1} (F - AU) + β_{1} {∥ U - X ∥}_{2}^{2} + β_{2} TV (X)) . \end{matrix}$ (12)

2.4 Optimization algorithm

To solve the cost function in Equation (12), we alternatively minimize U and X as $U = arg min_{U} {(F - AU)}^{T} Σ^{- 1} (F - AU) + β_{1} {∥ U - X ∥}_{2}^{2},$ (13) $X = arg min_{X} β_{1} {∥ U - X ∥}_{2}^{2} + β_{2} TV (X) .$ (14)

The purpose of this two-step minimization is straightforward. We first obtain an initial solution from Equation (13), which may be very noisy. Then, TV regularization is used to suppress the noise and artifacts, thereby effectively improving the quality of the basis material image.

2.4.1 Minimization with respect to U

Equation (13) is a quadratic optimization problem, and its solution can be obtained from the optimal condition $(A^{T} Σ^{- 1} A + β_{1} I) U = A^{T} Σ^{- 1} F + β_{1} X .$ (15)

Equation (15) admits a closed-form solution $U = {(A^{T} Σ^{- 1} A + β_{1} I)}^{- 1} (A^{T} Σ^{- 1} F + β_{1} X) .$ (16)

Although Equation (16) is an exact analytic solution, the enormous size of the matrix A^T Σ ^-1A + β₁I makes it impossible to store and invert A^T Σ ^-1A + β₁I exactly. Therefore, we propose employing a conjugate gradient (CG) method [42, 43] to solve the solution of Equation (15).

2.4.2 Minimization with respect to X

The minimum of the problem in Equation (14) can be solved by the Chambolle gradient projection (CGP) algorithm [44]. In CGP algorithm, we first solve the following constrained optimization problem $min_{| Y_{i, j} | \leq 1} {∥ U - β div Y ∥}_{2}^{2},$ (17) where Y_i,j is the dual variable at the (i, j)th pixel position, Y is the concatenation of all Y_i,j, β = 2β₂/β₁ , and div denotes the discrete divergence operator. When the minimum Y^*. of the constrained minimization problem in Equation (17) is determined, the solution of the problem in Equation (14) can be obtained as $X = U - β div Y^{*} .$ (18)

The iterative scheme for computing Y^* is given as $Y_{i, j}^{l + 1} = \frac{Y_{i, j}^{l} + λ β {(\nabla (β div Y^{l} - U))}_{i, j}}{1 + λ β | {(\nabla (β div Y^{l} - U))}_{i, j} |} .$ (19) where $Y_{i, j}^{l}$ is the lth iteration for the minimum, and λ > 0 is the step size used in the CGP algorithm; see reference [44] for more details.

2.4.3 Implementation

The implementation of the PWLS-TV material decomposition method is described in the following steps:

Initialization: X⁰, U⁰;

Initialization: β₁, β₂, λ, and k : =0

While stop criterion is not met;

Solve the problem in Equation (14): $U^{k + 1} = arg min_{U} {(F - AU)}^{T} Σ^{- 1} (F - AU) + β_{1} {∥ U - X^{k} ∥}_{2}^{2}$ with the initial value U^k using the CG algorithm;

U^k+1 = max (min(U^k+1, 1) , 0);

Solve the problem in Equation (15): $X^{k + 1} = arg min_{X} β_{1} {∥ U^{k + 1} - X ∥}_{2}^{2} + β_{2} TV (X)$ with the initial value X^k using the CP algorithm;

X^k+1 = max (min(X^k+1, 1) , 0);

End if criterion is satisfy.

In Step 1, the initial images X⁰ and U⁰ are set to zero. In Step 2, three parameters (i.e., β₁, β₂, and λ) should be initialized before implementation. The outer loop (Steps 4–7) is performed with two separated steps, i.e., CG step for a trial solution and CGP step for noise and artifacts suppression. Steps 5 and 7 are used to enforce the constraint in Equation (4). In Step 8, an adequate stopping criterion is selected to terminate the iteration.

2.5 Parameter selection

2.5.1 Selection of penalty parameters β₁ and β₂

There are two penalty parameters (β₁ and β₂) in the PWLS-TV material decomposition method. These parameters determine the smoothness of the decomposed result, and the values of β₁ and β₂ depend on the noise level of the measured TECT images. Clearly, a considerably small penalty parameter would result in a very noisy result, whereas a considerably large penalty parameter will unavoidably lead to an over smoothed result. Searching for the optimal penalty parameters remains a difficult task in CT imaging [45 –49].

In the experimental studies, we found that decomposed basis material images are less sensitive to the value of penalty parameter β₁ within the range of 1.0 × 10^-4 ≤ β₁ ≤ 1.0 × 10^-3. The value of β₂ considerably affects the smoothness of the basis material images; further, the value of β₂ is highly related to the noise level of the TECT images, and thus, it needs to be determined case by case.

2.5.2 Selection of iteration numbers for CG and CGP steps

The iteration numbers for each sub-iteration step is an essential factor for achieving high-quality basis material images. Although lager iteration numbers can obtain a convergent solution, the computational process may take a long time. Considering the characteristic of the alternating optimization algorithm, an adequate iteration number of CG steps is usually enough to generate a satisfactory result for the following CGP step. In this study, for all experimental data, the iteration numbers of CG and CGP steps are set to 2 and 10, respectively.

2.5.3 Stopping criterion

In this study, the relative error (RE) is used as the stopping criterion for the PWLS-TV material decomposition method. The RE is written as $RE = \sqrt{\frac{\sum_{n} {(U^{k + 1} (n) - U^{k} (n))}^{2}}{\sum_{n} {(U^{k} (n))}^{2}}}$ (20) where U^k+1 (n) and U^k (n) denote the value of the basis material image U^k+1 and U^k at pixel n, respectively. The iteration process is stopped when RE ≤ 1 . 0 × 10^-3 or the number of iteration reaches 200.

2.6 Experimental data

2.6.1 Digital phantom

The proposed PWLS-TV material decomposition method was first evaluated using a digital phantom as shown in Fig. 1. This phantom is composed of four materials: soft tissue (#1), Gadolinium (#2), bone (#3), and water (#4). The TECT images were obtained at 70 kVp, 100 kVp, and 140 kVp. The x-ray spectra of the corresponding energy were generated by the TASMIC method [50]. The normalized x-ray spectra are shown in Fig. 2. A total number of 1.0 × 10⁵ incident photons are assumed to be emitted from the tube of a fan-beam CT scanner. The x-ray source to the detector distance is 949 mm, and the rotation center to the detector distance is 408 mm. The detector has 888 channels with a size of 1.0 mm. A total number of 984 projection views are evenly acquired from a full scan. The TECT images are reconstructed by a conventional filtered back-projection (FBP) algorithm with a ramp filter. The dimension of the reconstructed image is 256×256, and the physical size is 1.0 mm for each pixel.

Fig.1

Digital phantom (left) with six objects and the geometric setting (right) of the phantom.

Fig.2

X-ray spectra used in the digital phantom study.

2.6.2 Mouse thorax phantom

The proposed PWLS-TV material decomposition method was evaluated using a realistic mouse thorax phantom [51], and 2.0% iodinated contrast agent (see Fig. 3). The x-ray spectra of the three different energies were generated by the TASMIC method [50]. The corresponding normalized x-ray spectra are shown in Fig. 4. A total of 5.0 × 10⁴ incident photons are assumed to be emitted from the tube of a fan-beam CT scanner. The imaging geometry of the CT scanner is the same as the digital phantom. The dimension of the reconstructed TECT image is 256×256 with the pixel size of 1.0×1.0 mm.

Fig.3

A slice of the thorax phantom (left) and iodinated contrast agent (right) distribution of this phantom.

Fig.4

X-ray spectra used in the mouse thorax phantom study.

2.7 Performance evaluation

2.7.1 Material decomposition accuracy

In this study, relative root mean square error (RMSE) is used to measure the material decomposition accuracy. RMSE can be expressed as $RMSE = \sqrt{\frac{\sum_{n} {(x (n) - x_{truth} (n))}^{2}}{\sum_{n} {(x_{truth} (n))}^{2}}}$ (21) where x (n) and x_truth (n) denote the intensities of the decomposed results and ground truth, respectively, at pixel n.

2.7.2 Image similarity measurement

Structural similarity (SSIM) [52] and feature similarity (FSIM) [53] are two widely used image quality assessment metrics, and they are used to measure the similarity between the decomposed results and the reference. The values of SSIM and FSIM range from 0 to 1; a larger value indicates a better similarity between the decomposed basis material image and the reference. In this study, we used the default parameters described in the corresponding references to calculate SSIM and FSIM, respectively.

2.7.3 Comparison methods

To further validate and evaluate the performance of the PWLS-TV material decomposition method, we compared it with the direct matrix inversion (DMI) method described in Equation (5), the TV filtration method for DMI decomposition results (TV-DMI), and the mean squared difference (MSD) method [16]. The TV-DMI method for material decomposition can be formulated as $min_{X} {(X - A^{- 1} F)}^{T} (X - A^{- 1} F) + β TV (X)$ (22) where β > 0 is the regularization parameter. The cost function of MSD material decomposition method is given as $min_{X} {∥ F - AX ∥}_{2}^{2}$ (23) where ∥ · ∥ ₂ denotes Euclidean norm.

3 Results

3.1 Digital phantom study

3.1.1 Visual evaluation

Figure 5 shows the reconstructed TECT images of the circle phantom with 70, 100, and 140 kVp x-ray tube voltages, respectively. The four circular objects show similar or same x-ray attenuations, which are considerably difficult to distinguish for different objects. In this section, we choose three basis materials namely, soft tissue, bone, and gadolinium contrast agent. Figure 6 shows the decomposed results, with columns 1–4 showing three basis material images and the corresponding color image, respectively. The results decomposed from the phantom by MSD method were viewed as the references. It can be observed that these two methods result in high noise level in decomposed results, especially in the soft tissue and color images. The noise has been effectively reduced compared with the basis material images decomposed using DMI and MSD methods; however, some undesired noise and artifacts (as indicated by the arrows in Fig. 6) still exist in TV-DMI results. The fifth row depicts decomposed results from TECT images by the proposed PWLS-TV method, where the severe noise is reduced while the edges information is well preserved. The presented PWLS-TV material decomposition method achieves the best performance in terms of noise suppression and edge detail preservation compared to the other methods.

Fig.5

TECT images of the digital phantom: (a) 70 kVp; (b) 100 kVp; and (c) 140 kVp. From left to right, the display winds are [0, 0.05], [0, 0.045], and [0, 0.04] mm^–1, respectively.

Fig.6

Basis material images from the digital phantom. From the left to right, the first three columns are soft tissue, bone, and gadolinium images, respectively. The fourth column is the color image where red, green, and blue components correspond to the three basis materials: soft tissue, bone, and Gadolinium, respectively. The first row is the results decomposed from phantom by MSD method. From the second to bottom rows, the results were decomposed by DMI, MSD, TV-DMI, and PWLS-TV methods, respectively.

3.1.2 Noise suppression

To quantitatively evaluate the performance of these different material decomposition methods, the RMSE values of the decomposed results are depicted in Fig. 7. Figure 7 indicates that the PWLS-TV decomposition has the smallest RMSE value, and DMI decomposition has the largest RMSE value. The average RMSE values of the three basis material images decomposed by the DMI, MSD, TV-DMI, and PWLS-TV methods are 20.85%, 15.87%, 9.58%, and 8.06%, respectively. The results illustrate that the PWLS-TV material decomposition method achieves considerable gains compared with the DMI, MSD, and TV-DMI methods in the aspect of noise suppression.

Fig.7

RMSE values of the decomposed results from digital phantom.

3.1.3 SSIM study

The SSIM values of the decomposed results are shown in Fig. 8 in which PWLS-TV method has the largest SSIM value, and DMI method has the smallest SSIM value. The average SSIM values of the three basis material images decomposed by the DMI, MSD, TV-DMI, and PWLS-TV methods are 0.25, 0.52, 0.45, and 0.77, respectively. The results demonstrate that PWLS-TV method achieves the best image quality over the DMI, MSD, and TV-DMI methods in terms of structure similarity.

Fig.8

SSIM values of the basis material images from digital phantom.

3.1.4 FSIM study

The FSIM values of the decomposed results are shown in Fig. 9. It can be seen that PWLS-TV decomposition has the largest FSIM value, followed by the TV-DMI and MSD methods, and DMI decomposition has the smallest FSIM value. The average FSIM values of the three basis material images decomposed by the DMI, MSD, TV-DMI, and PWLS-TV methods are 0.50, 0.62, 0.81, and 0.97, respectively. We confirm that PWLS-TV material decomposition method has the best performance in terms of feature similarity.

Fig.9

FSIM values of the decomposed results from digital phantom.

3.1.5 Convergence analysis

Figure 10 shows the convergence curves of the PWLS-TV decomposition with respect to RE and RMSE. Figure 10(a) and (b) depict the RE and RMSE curves with respect to iteration number, respectively. This figure indicates that the errors decrease gradually and remain stable after enough iterations.

Fig.10

Convergence curves of the PWLS-TV decomposition: (a) RE values versus iteration number; (b) RMSE values versus iteration number.

3.2 Mouse thorax phantom study

3.2.1 Visual evaluation

Figure 11 shows the TECT images of the mouse thorax phantom at 30, 50, and 70 kVp x-ray tube voltages, respectively. In this section, we choose soft tissue, bone, and iodinated contrast agent as the basis materials. Figure 12 shows the decomposed results, with columns 1–4 showing three basis material images and the corresponding color image, respectively. The first row depicts the decomposed results from the phantom by the MSD method, which is used as the reference. The second row shows the decomposed results from TECT images by DMI method, which results in high noise level in the decomposed results, especially in the soft tissue image, where fine structures in lung (as indicated by the red arrows in Fig. 12) is hardly distinguished from noise. The third row represents the decomposed results from TECT images by the MSD method, which suffers from severe degradation of image quality. Moreover, many pixels in the bone images (as indicated by the white arrows in Fig. 12) decomposed by DMI and MSD methods are misclassified as containing bone. The fourth row shows the decomposed results from TECT images by the TV-DMI method. The noise in these basis material images is effectively reduced; however, the edge of the resulting bone image are oversmoothed (as indicated by the blue arrow in Fig. 12). The fifth row shows decomposed results from TECT images by the proposed PWLS-TV method; not only is the severe noise reduced, but the edge and fine structure are well preserved. The abovementioned results reveal that PWLS-TV decomposition can obtain much better basis material images as compared to the result of the competing decomposition methods. The observation is consistent with our previous results in the digital phantom study.

Fig.11

TECT images of the mouse thorax phantom: (a) 30 kVp; (b) 50 kVp; (c) 70 kVp. From left to right, the display winds are [0, 0.16], [0, 0.07], and [0, 0.055] mm^–¹, respectively.

Fig.12

Basis material images from the mouse thorax phantom. From the left to right, the first three columns are the soft tissue, bone, and iodine images, respectively. The fourth column is the color image where red, green, and blue components correspond to the three basis materials: soft tissue, bone, and iodine, respectively. The first row is the results decomposed from phantom by MSD method. From the second to bottom rows, the results were decomposed by DMI, MSD, TV-DMI, and PWLS-TV methods, respectively.

3.2.2 Noise suppression

Figure 13 shows the RMSE values of the basis material images in Fig. 12. It can be seen that PWLS-TV decomposition has the smallest RMSE value. The average RMSE values of the three basis material images decomposed by the DMI, MSD, TV-DMI, and PWLS-TV methods are 29.03%, 25.73%, 26.33%, and 8.49%, respectively. That is, the material decomposition accuracy of TECT is greatly improved by the PWLS-TV method compared with DMI, MSD, and TV-DMI methods.

Fig.13

RMSE values of the decomposed results from the Mouse thorax phantom.

3.2.3 SSIM study

To further illustrate the differences between decomposed results, a zoomed ROI (indicated by the red square in Fig. 3) is shown in Fig. 14. The corresponding SSIM values of the zoomed ROI are shown in Fig. 15. For each basis material, PWLS-TV decomposition has the largest SSIM value, and the SSIM value decreases gradually from TV-DMI decomposition, MSD decomposition, to DMI decomposition. The average SSIM values of the three basis material images from DMI, MSD, TV-DMI, and PWLS-TV methods are 0.26, 0.33, 0.75, and 0.89, respectively. The results indicate that PWLS-TV decomposition has the best performance in terms of structure similarity.

Fig.14

Zoomed-in views of the ROI (indicated by the red square in Fig. 3) in decomposed results. The column and row have the same meaning as those in Fig. 12.

Fig.15

SSIM values of the images in Fig. 14.

3.2.4 FSIM study

The FSIM values of the zoomed ROI are shown in Fig. 16. It can be find that the PWLS-TV decomposition has the largest FSIM value, followed by the TV-DMI and MSD methods, and DMI decomposition has the smallest FSIM value in all basis materials. The average FSIM values of the three basis material images from DMI, MSD, TV-DMI, and PWLS-TV methods are 0.58, 0.65, 0.85, and 0.96, respectively. The results indicate that PWLS-TV decomposition achieves the best image quality over the DMI, MSD, and TV-DMI methods in terms of feature similarity.

Fig.16

FSIM values of the images in Fig. 14.

4 Discussion and conclusion

In this study, we presented an image-based material decomposition method for TECT using PWLS criteria with TV regularization. In contrast to existing material decomposition strategies, the PWLS-TV method enables simultaneous material decomposition and noise suppression in an iterative process. Moreover, the noise statistical properties of the image-based material decomposition process were fully exploited to enhance the quality of basis material images. The results from experimental studies demonstrate that the presented PWLS-TV method shows superior performance in noise reduction without sacrificing edge and fine structure details in the decomposed results.

In the digital phantom study, the PWLS-TV decomposition greatly improved material decomposition accuracy of TECT compared with DMI, MSD, and TV-DMI methods. In terms of the visual evaluation, PWLS-TV decomposition greatly reduced the noise in basis material images and generated satisfactory color images. The corresponding quantitative evaluation (RMSE, SSIM, and FISM studies) further showed the superior performance of PWLS-TV decomposition. In the mouse thorax phantom study, PWLS-TV method provided the most accurate basis material and color images compared with DMI, MSD, and TV-DMI methods. The PWLS-TV decomposition effectively reduced magnified noise and generated satisfactory results inside of the lung and boundary of the spine (see magnified view in Fig. 14) with visual inspection. For quantitative evaluation, the results of PWLS-TV decomposition were consistent with that of digital phantom study.

The convergence analysis in digital phantom study illustrated that the PWLS-TV decomposition will converge to a satisfactory solution with reasonable parameter settings. In the PWLS-TV decomposition, there are five main tuning parameters: the iteration number of the CG algorithm, the iteration number of the CGP algorithm, the step size λ of the CGP algorithm, β₁, and β₂. The iteration number of each sub-iteration plays an essential role in obtaining satisfactory results within an efficient time. A large iteration step can obtain a convergent solution; however, it may take a considerable amount of time to execute. In this study, for all experimental data, the iteration steps of the CG algorithm and the CGP algorithm were set to 2 and 10, respectively. It can be seen that satisfactory basis material images were obtained in the results section, which shows that the selected iteration step for each sub-iteration is suitable for PWLS-TV decomposition. The step size λ in the CGP algorithm can be determined according to the approach suggested in a previous study [44]. The penalty parameters β₁ and β₂ are highly related to the smoothness of the basis material images; a large one will generate over smoothing results, while a small one will obtain a noisy result. In this study, we first chose a large penalty value and then gradually decrease it to obtain a satisfactory result. This approach can be viewed as a trial-and-error process. Accurately determining the appropriate value for each parameter remains a challenge in CT imaging [45 , 54]. The strategy to adaptively determine the penalty parameters can be conducted to further improve the performance of the PWLS-TV decomposition.

In this study, the statistical weight was determined by the noise variance of the TECT images. In the implementation, we assumed a uniform distribution of the noise in TECT images, and calculated the noise variance of a manually selected ROI. Such an assumption has also been used in some previous image-domain based decomposition studies [20, 21]. Determining the exact distribution of the noise in TECT images remains a challenge, as it depends on a number of factors such as reconstruction algorithms. The performance of PWLS-TV decomposition method could be further improved using a pixel-dependent noise variation in the implementation of PWLS-TV decomposition. In the past few years, theoretical studies on computing pixel-dependent noise variance from CT images have been conducted. For example, Wunderlich and Noo developed a novel approach to analytically calculate the noise variance between pixel values in reconstructed CT images [55]. Based on the noise distribution of the projection data, Li et al. presented a computationally efficient method to reconstruct the noise map in the CT images [56]. Incorporating the pixel-dependent noise variance estimation into the presented PWLS-TV decomposition method will be an interesting topic for future study.

In summary, we developed a PWLS-TV material decomposition method to improve the quality of basis material images. The PWLS-TV method can simultaneously perform noise reduction and material decomposition in one iterative step. The PWLS-TV method effectively explores the noise statistical properties of the material decomposition process, and it results in a considerable improvement of basis material image quality.

Footnotes

Acknowledgments

This work was supported in part by the National Natural Science Foundation of China (11701097, 11661007), Training Plan for Excellent Young Talents in Jiangxi Province (20192BCB23019), Jiangxi Double Thousand Plan, Science and Technology Innovation Talent Program of Ganzhou, Project of Jiangxi Key Laboratory of Numerical Simulation Technology, US National Institutes of Health (R01 EB020366), Cancer Prevention and Research Institute of Texas (RP160661).

References

Niu

, et al., Total image constrained diffusion tensor for spectral computed tomography reconstruction, Applied Mathematical Modelling 68 (2019), 487–508.

Shikhaliev

P.M.

, Fritz

S.G.

, Photon counting spectral CT versus conventional CT: comparative evaluation for breast imaging application, Phys Med Biol 56(7) (2011), 1905–1930.

Willemink

M.J.

, et al., Photon-counting CT: technical principles and clinical prospects, Radiology 289(2) (2018), 293–312.

, et al., Spatial-spectral cube matching frame for spectral CT reconstruction, Inverse Problems 34(10) (2018), 104003.

, et al., Low-dose spectral CT reconstruction using image gradient ℓ0–norm and tensor dictionary, Applied Mathematical Modelling 63 (2018), 538–557.

Graser

, et al., Dual energy CT: preliminary observations and potential clinical applications in the abdomen, Eur Radiol 19(1) (2009), 13–23.

Johnson

T.R.

, et al., Material differentiation by dual energy CT: initial experience, Eur Radiol 17(6) (2007), 1510–1517.

Alvarez

R.E.

, Macovski

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

Kalender

W.A.

, et al., Evaluation of a prototype dual-energy computed tomographic apparatus, I. Phantom studies, Med Phys 13(3) (1986), 334–339.

10.

, Zhao

, Zhang

, Accurate iterative FBP reconstruction method for material decomposition of dual energy CT, IEEE Trans Med Imaging 38(3) (2018), 802–812.

11.

Zhang

, et al., Iterative reconstruction for dual energy CT with an average image-induced nonlocal means regularization, Phys Med Biol 62(13) (2017), 5556–5574.

12.

Patino

, et al., Material separation using dual-energy CT: current and emerging applications, RadioGraphics 36(4) (2016), 1087–1105.

13.

Long

, Fessler

J.A.

, Multi-material decomposition using statistical image reconstruction for spectral CT, IEEE Trans Med Imaging 33(8) (2014), 1614–1626.

14.

Mendonca

P.R.

, Lamb

, Sahani

D.V.

, A flexible method for multi-material decomposition of dual-energy CT images, IEEE Trans Med Imaging 33(1) (2014), 99–116.

15.

McCollough

C.H.

, et al., Dual- and multi-energy CT: principles, technical approaches and clinical applications, Radiology 276(3) (2015), 637–653.

16.

Granton

P.V.

, et al., Implementation of dual- and triple-energy cone-beam micro-CT for postreconstruction material decomposition, Med Phys 35(11) (2008), 5030–5042.

17.

Sato

, et al., Triple-energy high-count-rate X-ray computed tomography scanner using a cadmium telluride detector, Health and Technology 8(3) (2018), 197–203.

18.

Zhao

, et al., Material decomposition using triple-energy CT for accurate proton therapy dose calculation, International Journal of Radiation Oncology Biology Physics 99(2) (2017), S176.

19.

Zhao

, et al., A unified material decomposition framework for quantitative dual- and triple-energy CT imaging, Med Phys 45(7) (2018), 2964–2977.

20.

Niu

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

21.

Xue

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

22.

Harms

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

23.

Tao

, et al., Material decomposition with prior knowledge aware iterative denoising (MD-PKAID), Phys Med Biol 63(19) (2018), 195003.

24.

Zeng

, et al., Penalized weighted least-squares approach for multienergy computed tomography image reconstruction via structure tensor total variation regularization, Comput Med Imaging Graph 53 (2016), 19–29.

25.

Zeng

, et al., Spectral CT Image Restoration via an Average Image-Induced Nonlocal Means Filter, IEEE Trans Biomed Eng 63(5) (2016), 1044–1057.

26.

Foygel Barber

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

27.

Liu

, et al., TICMR: total image constrained material reconstruction via nonlocal total variation regularization for spectral CT, IEEE Trans Med Imaging 35(12) (2016), 2578–2586.

28.

, et al., Improved material decomposition with a two-step regularization for spectral CT, IEEE Access 7 (2019), 158770–158781.

29.

Sidky

E.Y.

, Zou

, Pan

, Impact of polychromatic x-ray sources on helical, cone-beam computed tomography and dual-energy methods, Phys Med Biol 49(11) (2004), 2293–303.

30.

Ducros

, et al., Regularization of nonlinear decomposition of spectral x-ray projection images, Med Phys 44(9) (2017), e174–e187.

31.

Zhao

, et al., Using edge-preserving algorithm with non-local mean for significantly improved image-domain material decomposition in dual-energy CT, Phys Med Biol 61(3) (2016), 1332–1351.

32.

Niu

, et al., Nonlocal low-rank and sparse matrix decomposition for spectral CT reconstruction, Inverse Problems 34(2) (2018), 024003.

33.

Niu

, et al., Iterative reconstruction for photon-counting CT using prior image constrained total generalized variation, Computers in Biology and Medicine 103 (2018), 167–182.

34.

Yao

, et al., Multi-energy computed tomography reconstruction using a nonlocal spectral similarity model, Phys Med Biol 64(3) (2019), 035018.

35.

Dong

, Niu

, Zhu

, Combined iterative reconstruction and image-domain decomposition for dual energy CT using total-variation regularization, Med Phys 41(5) (2014), 051909.

36.

Heismann

B.J.

, Spectral Computed Tomography. 2012.

37.

Janos

, Ragnar

, Tommy

, Measurement of bone mineral using multiple-energy x-ray absorptiometry, Phys Med Biol 43(2) (1998), 379.

38.

Sidky

E.Y.

, Kao

C.M.

, Pan

X.H.

, Accurate image reconstruction from few-views and limited-angle data in divergent-beam CT, J Xray Sci Technol 14(2) (2006), 119–139.

39.

Sidky

E.Y.

, Pan

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

40.

Tang

, Nett

B.E.

, Chen

G.H.

, Performance comparison between total variation (TV)-based compressed sensing and statistical iterative reconstruction algorithms, Phys Med Biol 54(19) (2009), 5781–5804.

41.

Chen

G.H.

, Tang

, Leng

, Prior image constrained compressed sensing (PICCS): a method to accurately reconstruct dynamic CT images from highly undersampled projection data sets, Med Phys 35(2) (2008), 660–663.

42.

, Guan

, Wei

, Globally convergent Polak–Ribière–Polyak conjugate gradient methods under a modified Wolfe line search, Applied Mathematics and Computation 215(8) (2009), 3082–3090.

43.

, Zhao

, Wei

, A descent nonlinear conjugate gradient method for large-scale unconstrained optimization, Applied Mathematics and Computation 187(2) (2007), 636–643.

44.

Chambolle

, An algorithm for total variation minimization and applications, Journal of Mathematical Imaging and Vision 20(1–2) (2004), 89–97.

45.

, et al., Optimizing a Parameterized Plug-and-Play ADMM for iterative low-dose CT reconstruction, IEEE Trans Med Imaging 38(2) (2019), 371–382.

46.

Shen

, et al., Intelligent parameter tuning in optimization-based iterative CT reconstruction via deep reinforcement learning, IEEE Trans Med Imaging 37(6) (2018), 1430–1439.

47.

Wang

, Guan

, Solberg

, Inverse determination of the penalty parameter in penalized weighted least-squares algorithm for noise reduction of low-dose CBCT, Med Phys 38(7) (2011), 4066–4072.

48.

, et al., Projection decomposition algorithm for dual-energy computed tomography via deep neural network, J Xray Sci Technol 26(3) (2018), 361–377.

49.

Zeng

, Guo

, Liu

, Limited-angle cone-beam computed tomography image reconstruction by total variation minimization and piecewise-constant modification, Journal of Inverse and Ill-Posed Problems 21(6) (2013), 735–754.

50.

Hernandez

A.M.

, Boone

J.M.

, Tungsten anode spectral model using interpolating cubic splines: unfiltered x-ray spectra from 20 kV to 640 kV, Med Phys 41(4) (2014), 042101.

51.

Segars

W.P.

, et al., Development of a 4-D digital mouse phantom for molecular imaging research, Molecular Imaging & Biology 6(3) (2004), 149–159.

52.

Wang

, et al., Image quality assessment: from error visibility to structural similarity, IEEE Trans Image Process 13(4) (2004), 600–612.

53.

Zhang

, et al., FSIM: a feature similarity index for image quality assessment, IEEE Trans Image Process 20(8) (2011), 2378–2386.

54.

Niu

, et al., Iterative reconstruction for sparse-view x-ray CT using alpha-divergence constrained total generalized variation minimization, J Xray Sci Technol 25(4) (2017), 673–688.

55.

Wunderlich

, Noo

, Image covariance and lesion detectability in direct fan-beam x-ray computed tomography, Phys Med Biol 53(10) (2008), 2471–2493.

56.

, et al., Adaptive nonlocal means filtering based on local noise level for CT denoising, Med Phys 41(1) (2014), 011908.

Statistical image-based material decomposition for triple-energy computed tomography using total variation regularization

Abstract

BACKGROUND:

OBJECTIVE:

METHODS:

RESULTS:

CONCLUSIONS:

Keywords

1 Introduction

2 Methods and materials

2.1 Image-based basis material decomposition for TECT

2.5 Parameter selection

2.5.1 Selection of penalty parameters β1 and β2

2.5.2 Selection of iteration numbers for CG and CGP steps

2.5.3 Stopping criterion

2.6.1 Digital phantom

2.7.1 Material decomposition accuracy

2.7.3 Comparison methods

3.1 Digital phantom study

3.1.1 Visual evaluation

3.2.1 Visual evaluation

Footnotes

Acknowledgments

References

2.5.1 Selection of penalty parameters β₁ and β₂