Iterative reconstruction for low dose dual energy CT using information-divergence constrained spectral redundancy information

Abstract

Dual energy computed tomography (DECT) can improve the capability of differentiating different materials compared with conventional CT. However, due to non-negligible radiation exposure to patients, dose reduction has recently become a critical concern in CT imaging field. In this work, to reduce noise at the same time maintain DECT images quality, we present an iterative reconstruction algorithm for low-dose DECT images where in the objective function of the algorithm consists of a data-fidelity term and a regularization term. The former term is based on alpha-divergence to describe the statistical distribution of the DE sinogram data. And the latter term is based on the redundant information to reflect the prior information of the desired DECT images. For simplicity, the presented algorithm is termed as “AlphaD-aviNLM”. To minimize the associative objective function, a modified proximal forward-backward splitting algorithm is proposed. Digital phantom, physical phantom, and patient data were utilized to validate and evaluate the presented AlphaD-aviNLM algorithm. The experimental results characterize the performance of the presented AlphaD-aviNLM algorithm. Speficically, in the digital phantom study, the presented AlphaD-aviNLM algorithm performs better than the PWLS-TV, PWLS-aviNLM, and AlphaD-TV with more than 49%, 34%, and 40% gains for the RMSE metric, 1.3%, 0.4%, and 0.7% gains for the FSIM metric and 13%, 8%, and 11% gains for the PSNR metric. In the physical phantom study, the presented AlphaD-aviNLM algorithm performs better than the PWLS-TV, PWLS-aviNLM, and AlphaD-TV with more than 0.55%, 0.07%, and 0.16% gains for the FSIM metric.

Keywords

Dual energy computed tomography low-dose renconstruction alpha-divergence redundant information

1 Introduction

Dual energy computed tomography (DECT) has attracted increasing attention because of its ability to provide two sets of energy-specific data at the same time [1]. And DECT can derive additional attenuation properties of the object, including the attenuation values at different energies as well as discrimination and quantification of materials within the object, which may provide diagnostic information beyond that is obtainable with conventional CT [2, 3]. However, there is a trade-off between DECT images quality and the amount of radiation exposure in DECT imaging to patients. Without adequate treatments, the conventional filtered back-projection (FBP) algorithm would produce images with increasing noise from the DECT projection data under low-milliampere-second (mAs) level. Besides, the direct inversion operation that derives the material-selective images from DECT images is highly sensitive to noise and could result in excessive noise-induced artifacts in the material-selective images. Therefore, to address these issues, in this work, we focus on robust low-dose DECT images reconstruction.

Recent efforts have devoted to radiation exposure reduction in DECT imaging while maintaining the spatial resolution and quantitative accuracy. Numerous advanced DECT images processing algorithms have been proposed to obtain high-quality DECT images at low radiation dosage, such as projection- and image-domain filtering algorithms [4 –8], and statistical iterative reconstruction (SIR) algorithms [9 –13]. Among them, the SIR algorithms have shown superior performance in reducing noise and improving resolution. In general, the SIR algorithms often incorporate an accurate system model, statistical noise model, and prior model and are formulated by an objective function consisting of a data-fidelity term and a regularization term [12 , 14– 17]. The data-fidelity term is a fitting model of the observed DECT projection data and quantifies how well the reconstructed DECT images fit to the DECT measurements [18]. The accurate DECT projection measurement statistic distribution is the prerequisite for the data-fidelity term building. For polyenergetic X-ray source, the CT measurements follow compound Poisson distribution plus Gaussian electronic noise [19, 20], so that the model under polyenergetic X-ray source can be well-approximated by a simple Poisson plus Gaussian model. Under this assumption, the data fidelity can be built and it can measure the difference between the measured DECT projections and corresponding estimates. It is reported that many criteria can be used for describing the data fidelity, such as weighted least squares (WLS) criterion, maximum likelihood (ML) criterion and information divergence criterion [18, 21]. In this work, we propose the use of information divergence theory to describe the statistical distribution of the measured DECT projection data.

On the other hand, the regularization term plays a critical role for successful DECT images reconstruction. One primary strategy for regularizing reconstruction from DECT projection data is based on total variation (TV) models, which is high effectiveness to preserve edges and recover the smooth regions by exploiting local structure modes [11 , 23]. For example, T. Szczykutowicz et al. achieved substantial noise suppression on DECT images by using a prior image constrained compressed sensing (PICCS) algorithms [11]. The PICCS enforces gradient sparsity relative to a prior reconstruction, recognizing that most image features are strongly similar (redundant) among DECT images. However, the TV models usually lose the image details and tend to over-smooth effects due to the piecewise constant assumption [14 , 23– 25]. Another class family of the regularization term is based on patch-based models wherein small patches are extracted from an image and processed individually, and then the individually processed patches are reassembled to obtain the final output image [9 , 26]. For instance, Ma et al. proposed a prior-image induced NLM-based (ndiNLM) regularization for perfusion CT [27]. The patch similarity weights were calculated between patches of the low-dose reconstructed images and those of the prior images. Compared with pixel-based models, the patch-based models effectively captures image features by employing geometrical similarities existing in the desired images [28, 29]. Moreover, the patch-based models reported very good reconstruction results, outperforming conventional regularization terms such as TV models [27].

This work develops an iterative algorithm for radiation reduction in the low-mAs DECT imaging by using information-divergence constrained spectral redundancy information. In particular, the DECT images are yielded by minimizing one objective function which consists of an alpha-divergence [21] based data-fidelity term to describe the statistical distribution of the DE sinogram data, and a DECT images redundant information [10, 27] based regularization term to reflect the prior information of the desired DECT images. To distinguish it from the ndiNLM model and emphasize the use of strong similarity among DECT images, the presented algorithm is termed as “AlphaD-aviNLM” for simplicity. Subsequently, to minimize the associative objective function, a modified proximal forward-backward splitting algorithm is proposed.

2 Methods and materials

2.1 Statistical model in DECT sinogram domain

It is noted that DECT system scans the object using two energy spectra with different kVp settings and it can obtain high- and low- energy sinogram data simultaneously. In this study, the DECT measurements at distinct projections are assumed to be conditionally independent. Mathematically, the X-ray CT measurement at either high- or low- energy spectra can be approximately expressed a following discrete linear system [14, 27]: $\bar{y} = Hx,$ (1) where $\bar{y} = ({\bar{y}}_{1}, {\bar{y}}_{2}, . . ., {\bar{y}}_{M})^{'}$ denotes the vector of expected line integrals data, and x = (x₁, x₁, . . . , x_N) ′ represents the vector of attenuation coefficients, where M and N denote the number of line integral measurements and image pixels, respectively. H denotes the system matrix with the size M × N, and its element H_i,j is typically calculated as the intersection length of projection ray i with j.

In our previous studies of low-dose CT data statistics, the line integral measurements follow a non-stationary Gaussian distribution and can be treated as normally distributed with a nonlinear signal-dependent variance [20]. Assuming the line integral measurements among all the bins are statistically independent, the likelihood function of the joint probability distribution can be written as follows [15]: $φ (y) = \prod_{i} \exp [- \frac{(y_{i} - {\bar{y}}_{i})^{2}}{2 σ_{i}^{2}}],$ (2) where y = (y₁, y₂, . . . , y_M) ′ denotes the vector of measured line integral data. $σ_{i}^{2}$ denotes the electronic noise background and can be determined by the mean-variance relationship proposed by Ma et al. [20]. The statistical model in Eq. (2) would become more complicated when the polychromatic nature of the X-ray source is considered. Alternatively, the information divergence criterion is an attractive measure to describe the statistical properties of the sinogram data and can be introduced into objective function building.

2.2 Information divergence

As one type of information divergence, Alpha-divergence is particularly suitable to fit the low-dose CT sinogram data because of its adaptive nature and it has been proven that the Alpha-divergence can be adapted to measure the discrepancy between sinogram data and corresponding estimates for CT image reconstruction [21]. In this study, to find an approximate y = Hx, the Alpha-divergence theory is utilized to measure discrepancy between y and Hx, which is written as follows: $D_{α} (y, Hx) = \frac{1}{α (1 - α)} \sum_{i = 1}^{M} [α y_{i} + (1 - α) (Hx)_{i} - y_{i}^{α} {(Hx)}_{i}^{1 - α}],$ (3) where Alpha-Divergence is parameterized by α. In this work, α is set to be 1.2 empirically by trying several values.

Although an iterative procedure can effectively suppress image noise to some extent, it could fail when a smoother estimate is required. Ma et al. proved that Eq. (3) could favor a particular solution based on a prior knowledge [21, 30]. Consequently, the cost function with a regularization term R (x) can be formulated as follows: $min_{i \geq 0} D_{α} (y, Hx) + β R (x),$ (4) where β denotes a trade-off parameter to balance the data-fidelity term and regularization term. Since we assume that the measurements at high- and low- energy are made independently. In this study, we propose the following Alpha-divergence constrained spectral redundancy information minimization for DECT image reconstruction: $min_{X \geq 0} D_{α} (Y, HX) + \overset{⇀}{β} R (X),$ (5) where Y = [y_h, y_l] ^T and y_h, y_l denote the measured high and low -energy sinogram data, respectively. X = [x_h, x_l] ^T and x_h, x_l denote obtained high and low -energy image, respectively. $\overset{⇀}{β} = {[β_{h}, β_{l}]}^{T}$ and β_h and β_l denote the trade-off parameters.

2.3 aviNLM regularization

It is noted that the choice of the penalty term in the objective function greatly affects the quality of the resultant estimates. Ma et al. proposed a previous normal-dose scan induced NLM-based regularization to improve the follow-up low-dose CT image quality [27]. Subsequently, Zeng et al. developed an averaged-image induced nonlocal means (aviNLM) regularization for DECT image reconstruction [8, 10], where the regularization utilized the redundant spectral information among DECT images because of strong similarity (redundancy) among DECT images [8, 11]. It is noted that the same characteristics existing in these NLM-based models is that they incorporate the redundant information in high-quality prior images, which can obtain remarkable improvements. Specifically, the aviNLM models can be described as follows: $aviNLM (x_{h / l} (i)) = \sum_{j \in N_{i}} w_{h / l} (i, j) x_{av} (j),$ (6) where N_i represents the search window in the image x_av which is produced by averaging the DECT images acquired at two distinct X-ray energy spectra. The weight w (i, j) is a decreasing function of the similarity between two local neighborhoods centered at the voxel i in the object image x_h/l and the voxel j in the prior image x_av. Mathematically, w (i, j) can be defined as follows: $w_{h / l} (i, j) = \frac{C (i, j) \cdot \exp {- {∥ x_{h / l} (V_{i}) - C \cdot x_{av} (V_{j}) ∥}_{2, a}^{2} / p^{2}}}{\sum_{j \in N_{i}} \exp {- {∥ x_{h / l} (V_{i}) - C \cdot x_{av} (V_{j}) ∥}_{2, a}^{2} / p^{2}}},$ (7) where V_i and V_j denote two similarity patch-windows at the voxel i in the object image x_h/l and the voxel j in the prior image x_av. The notation ∥ · ∥ _2,α denotes a Gaussian-weighted Euclidean distance metric between two similarity patch-windows, where a denotes the standard deviation of Gaussian function. The parameter p denotes the decay of weighting coefficient as the function of two patches distance.

In DECT imaging, the different kVp leads to the attenuation coefficient changes in the energy-specific image. Therefore, a local compensation factor C is needed in Eq. (7) to account for local intensity change [8]. In this study, the compensation factor C is defined as follows: $C (x_{h / l} (V_{i}), x_{av} (V_{j})) = \frac{E (x_{h / l} (V_{i}))}{E (x_{av} (V_{j}))},$ (8) where E (·) denotes the expected value or mean of the intensity in the patch-window V.

In this work, the aviNLM regularization is introduced and it is defined as follows: $R (X) = φ_{p} (X - aviNLM (X)),$ (9) where φ_p denotes positive potential function. Specifically, in this study, a non-quadratic function is utilized to preserve the edge details, which can be defined as follows: $φ_{p} (Δ) = ∥ Δ ∥_{p}^{p} 1 < p < 2,$ (10) where ∥ · ∥ _p denotes the p-norm of the discrete magnitude of the image. p can be regarded as a shape factor controlling the cost of abrupt edges and in this study, p is set to 1.5.

2.4 Optimization approach

To find the solution of the cost function in Eq. (5), a modified proximal forward-backward splitting algorithm is utilized in this study and the convergence result of the proximal forward-backward splitting algorithm has been given in Combettes and Wajss work [31]. Mathematically, minimizing Eq. (5) can be split as two sub-problems: $\begin{matrix} (P 1) : [\begin{matrix} u_{h}^{k + 1} \\ u_{l}^{k + 1} \end{matrix}] = [\begin{matrix} x_{h}^{k} - \frac{1}{μ_{h}} \nabla D_{α} (y_{h}, {Hx}_{h}^{k}) \\ x_{l}^{k} - \frac{1}{μ_{l}} \nabla D_{α} (y_{l}, {Hx}_{l}^{k}) \end{matrix}] \\ (P 2) : [\begin{matrix} x_{h}^{k + 1} & x_{l}^{k + 1} \end{matrix}] = arg min_{x_{h}, x_{l}} [\frac{μ_{h}}{2 β_{h}} {∥ x_{h} - u_{h}^{k + 1} ∥}_{2}^{2} + R (x_{h}) \\ + \frac{μ_{l}}{2 β_{l}} {∥ x_{l} - u_{l}^{k + 1} ∥}_{2}^{2} + R (x_{l})] \end{matrix},$ (11)

where k denotes the iteration index. μ_h and μ_l are auxiliary vectors and they are positive and constant. By the intuitive perspective, the solution of (P1) is used as the initial value in minimizing the objective function of (P2) and vice versa. As a result, Eq. (5) can obtain acceptable solution with reasonable noise-resolution trade-off. In the implementation, a surrogate function (SF) algorithm was taken to minimize the objective function of (P1), and the implementation details of the SF algorithm were performed as follows: $u_{h / l, j}^{k + 1} = x_{h / l, j}^{k} {\sum_{i = 1}^{M} h_{ij} {(\frac{y_{h / l, i}}{{({Hx}_{h / l}^{k})}_{i}})}^{α}}^{1 / α}, j = 1, 2, . . ., N .$ (12) If α = 0, Eq. (12) can be updated as follows: $u_{h / l, j}^{k + 1} = x_{h / l, j}^{k} \prod_{i = 1}^{M} {(\frac{y_{h / l, i}}{({Hx}_{h / l}^{k})_{i}})}^{h_{ij}}, j = 1, 2, . . ., N .$ (13) The subproblem (P2) is essentially an aviNLM denoising problem [8]. After obtaining optimal reconstructed DECT images, an image-based basic material decomposition algorithm is applied to yield basic material images [1].

2.4.1 Selection of the parameter α

The parameter α is important parameter in the objective function in Eq. (5). Therefore, how to select the optimal value of α is a vital task. It is noted that if α = 2, the Alpha-divergence can be turned into the WLS model. And if α = 1, the Alpha-divergence can be turned into the KL divergence based on the Minka’s fixed-point scheme [32]. In this study, we only focus on the α value between 1 and 2. However, it is difficult to define the “optimal” α value at this moment, and in this study, we have set the α value in an empirical fashion.

2.4.2 Selection of the parameters in the aviNLM model

Determining the optimal parameters in the aviNLM model is not a trivial task. In this work, all the parameters in the aviNLM model were empirically determined by trying several combinations of parameters through extensive experiments.

2.4.3 Selection of the parameters $\overset{⇀}{β}$

The parameter $\overset{⇀}{β}$ is used to control the tradeoff between the datafidelity and the regularization term in the objective function. As is known, it is a difficult task to optimize it in the DECT images reconstruction. The parameter $\overset{⇀}{β}$ is chosen by a trial-end-error fashion for best image quality. The range of the parameter β_h is between 0.05 and 0.15, and the range of the parameter β_l is between 0.03 and 0.12 in this study.

2.5 Experimental data acquisition

In this study, a digital XCAT phantom [33], a physical phantom and real patient data were utilized to validate and evaluate the performance of the AlphaD-aviNLM algorithm in low-dose DECT images reconstruction.

2.5.1 XCAT phantom data

A digital XCAT phantom [33] was used in this study as shown in Fig. 1. The phantom consists of three different materials mimicking bone, average soft tissue, and air. Their densities are listed in Table 1. In our study, the low-dose dual energy (i.e., 140 and 80 kVp) sinogram data of the XCAT phantom were acquired using the simulator model described in [34] with a scan geometry. The scan geometry was chosen according to the Siemens Somatom Sensation 16 CT scanner (Siemens Healthcare, Forchheim, Germany), which is illustrated in Table 2.

Fig.1

Illustration of a digital XCAT phantom (a) and physical phantom (b).

Table 1

Material types in the digital XCAT phantom shown in Fig. 1

Number	Material	Density (g/cm³)
1	Air	0.00
2	Tissue	0.99
3	Bone	1.90

Table 2

The system geometry of the Siemens Somatom Sensation 16 CT scanner

Parameter
Detector bin spacing	1.407 mm
Number of detector bins	672
Projection views per rotation	1160
Source-to-detector distance	1040 mm
Object-to-detector distance	470 mm

2.5.2 Physical phantom data

In this study, the physical phantom (Fig. 1(b)) projection data were collected on a commercial CT scanner (Brilliance Big Bore, Philips) with two different tube potentials (90 and 120 kVp) at two noise levels (282 and 56 mAs), respectively. The phantom consists of a circular water background with the diameter of 15 cm, eight large circular inserts (C1: Acrylic, 1.147 g/cm³, C2: Delrin^TM, 1.368 g/cm³, C3: Teflon, 1.868 g/cm³, C4, and B8: Air, 0.00 g/cm³, C5: PMP, 0.858 g/cm³, C6: LDPE, 0.945 g/cm³, C7: Polystyrene, 0.998 g/cm³), and four small circular inserts (S1, S3, and S4: Air, 0.00 g/cm³; S2: Teflon pin, 1.868 g/cm³). In this study, the DECT projection data at lower noise level (i.e. 282 mAs) are selected as the normal dose one which can be served as the ground truth.

2.5.3 Patient data

The patient data were obtained with patient consents for two chest CT studies (patient 1 with suspected coronary atherosclerotic plaque, patient 2 with pulmonary nodules) for medical reasons and the clinical data serve as a pilot clinical study. The two patients were scanned by the GE Discovery CT750 HD scanner. The virtual monochromatic spectral (VMS) images were generated using the GE commercial software at 10 keV monochromatic energy level increments from 40-140 keV. In patient 1 study, the VMS images at 40 keV and 120 keV were specifically selected as the DECT images, and in patient 2 study, the VMS images at 70 keV and 140 keV were selected. The low-dose DECT sinogram data were acquired using the simulation method in [33].

2.6 Comparison Methods

To validate and evaluate the performance of the presented AlphaD-aviNLM algorithm, the TV-based penalized weighted least-squares (PWLS-TV), PWLS-aviNLM [10], and AlphaD-TV algorithms were adopted for comparison, and the cost function of the PWLS-TV, PWLS-aviNLM, and Alpha-TV algorithm can be written as follows: $\begin{matrix} min_{X \geq 0} {(Y - HX)}^{T} W (Y - HX) + {\overset{⇀}{η}}_{1} TV (X), \end{matrix}$ (14) $\begin{matrix} min_{X \geq 0} {(Y - HX)}^{T} W (Y - HX) + {\overset{⇀}{η}}_{2} aviNLM (X), \end{matrix}$ (15) $min_{X \geq 0} D_{α} (Y, HX) + {\overset{⇀}{η}}_{3} TV (X),$ (16) where $W = {[\sum_{h}^{- 1}, \sum_{l}^{- 1}]}^{T}$ and ∑_h/l is energy-specific diagonal matrix with the ith element of $σ_{h / l, i}^{2}$ which is the variance of sinogram data ${\overset{⇀}{η}}_{t} = {[η_{th}, η_{tl}]}^{T} (t = 1, 2, 3)$ denotes the trade-off parameters.

3 Results

3.1 XCAT Phantom Study

Figure. 2 shows the ground truth and the digital XCAT phantom results reconstructed by different algorithms from low-dose DECT measurements, including the FBP, PWLS-TV, PWLS-aviNLM, AlphaD-TV and presented AlphaD-aviNLM algorithms. It can be seen that the statistical iterative reconstruction (SIR) algorithms have shown very promising results in suppressing the noise the artifacts in low-dose DECT compared with conventional FBP algorithm wherein severe noise-induced artifacts exist in the FBP DECT images. The TV-based SIR algorithms, including the PWLS-TV and AlphaD-TV algorithms, can reduce such noise-induced artifacts to some extent. However, it is noticeable that some of the anatomical structures are blurred due to some patchy artifacts indicated by the red arrows. We also can see that the aviNLM-based SIR algorithms outperform the TV-based SIR algorithms in terms of noise-induced artifacts reduction and image quality improvement. Moreover, the presented AlphaD-aviNLM reconstructed DECT images are sharper than the TV-based ones and slightly sharper than PWLS-aviNLM ones. To further illustrate the gains of the presented AlphaD-aviNLM algorithm, two regions of interest (ROIs) indicated by the red squares were selected, which have abundant detailed features. Figure. 3 shows the magnified ROIs for a more detailed analysis. From the results, it can be concluded that the presented Alphad-aviNLM algorithm performs much better than the TV-based SIR algorithms and slightly better than the PWLS-aviNLM algorithms in terms of effective noise suppression and low-contrast objects preservation.

Fig.2

The digital XCAT DECT images reconstructed by the FBP (second row), PWLS-TV (third row), PWLS-aviNLM (fourth row), AlphaD-TV (fifth row), and presented AlphaD-aviNLM (sixth row) algorithms from low-dose measurements. All the images are displayed in the same window.

Fig.3

The two magnified ROIs as indicated by the red squares in Fig 2: the top two rows are from ROI 1, and the bottom two rows are from ROI 2. From A to F: Phantom, the FBP, PWLS-TV, PWLS-aviNLM, AlphaD-TV, and presented AlphaD-aviNLM algorithms. All the images are displayed in the same window.

The image resolution can be characterized by the horizontal profile labeled by a cyan line in Fig 2. From the Fig 4, it can be observed that the presented AlphaD-aviNLM algorithm follows the closest to the phantom profile among all the algorithms. To give a more detailed comparison, an enlarged part of profile indicated by the cyan squares can demonstrate more noticeable edge preservation can be achieved by the presented AlphaD-aviNLM algorithm than other four algorithms.

Fig.4

Profile comparison along the horizontal cyan line labeled in Fig 2: (a): 140 kVp, and (b): 80 kVp

To provide quantitative and statistical comparison results, we compared the root mean squared error (RMSE) and the feature similarity (FSIM) assessments among the four SIR algorithms in terms of global image and local ROIs indicated by the red squares in Figure. 2. Figure. 5 illustrates the quantitative RMSE and FSIM assessments of the reconstructions in terms of global image and local ROIs. The results can demonstrate that the level of agreement with the phantom data is, in descending order: the presented AlphaD-aviNLM, PWLS-aviNLM, AlphaD-TV and PWLS-TV for all cases. Moreover, the presented AlphaD-aviNLM algorithm exhibits an average of more than 34%, 40%, and 49% gains for the RMSE metric and 0.4%, 0.7%, and 1.3% gains for the FSIM metric over the PWLS-aviNLM, AlphaD-TV and PWLS-TV, respectively. The RMSE and FSIM measurements are consistent with the observations in the visualization-based evaluation study.

Fig.5

The RMSE and FSIM measurements of four different SIR algorithms on global image and local ROIs indicated by the red squares in Fig. 2.

To evaluate the noise reduction performance of the different SIR algorithms, two quantitative metrics were employed and the corresponding global peak signal-to-noise ratio (PSNR) and local standard deviation (STD) scores are listed in Tables 3 and 4. From the results, the presented AlphaD-aviNLM algorithm can obtain the largest PSNR and smallest STD values in all cases. Speficically, the presented AlphaD-aviNLM algorithm performs better than the PWLS-TV, PWLS-aviNLM, and AlphaD-TV with more than 13%, 8%, and 11% gains for the PSNR metric. Therefore, the results demonstrate the presented AlphaD-aviNLM algorithm has the best performance in noise reduction, confirming the visual observations.

Table 3

The Global PSNR measurements of the different SIR algorithms for DECT images

	FBP	PWLS-TV	PWLS-aviNLM	AlphaD-TV	AlphaD-aviNLM
140 kVp	29.21	32.06	33.30	32.18	37.95
80 kVp	27.90	33.53	35.60	34.21	37.46

Table 4

The local STD measurements of the different SIR algorithms for DECT images

		FBP	PWLS-TV	PWLS-aviNLM	AlphaD-TV	AlphaD-aviNLM
ROI A	140 kVp	5.35e-3	3.13e-3	2.15e-3	2.60e-3	1.41e-3
	80 kVp	5.76e-3	3.20e-3	2.90e-3	2.89e-3	2.35e-3
ROI B	140 kVp	6.44e-3	3.35e-2	2.41e-2	3.45e-2	1.61e-2
	80 kVp	6.01e-3	3.48e-2	3.13e-2	3.62e-2	2.27e-2

To further demonstrate the performance of the presented AlphaD-aviNLM algorithm, Fig. 6 shows the decomposed images of basic materials from different algorithms. In the implementation, the tissue and bone were chosen as the basis materials. In can be seen that the direct matrix inversion of the low-dose FBP images results in severe noise-induced artifacts in the decomposed images, but the AlphaD-aviNLM algorithm performs much better than the other three algorithms in terms of noise suppression and edge preservation, which is indicated by the red arrows. In addition, the mean and STDs of the ROIs indicated by the blue boxes in Fig. 6 are calculated and the related results are listed in Table 5. The results demonstrated that the presented AlphaD-aviNLM algorithm can obtain the lowest noise level in the basis material images among all the algorithms. In other words, the presented AlphaD-aviNLM algorithm can outperform the other algorithms in differentiating materials with the low-dose data acquisition.

Fig.6

The bone and tissue images reconstructed by different algorithms. All the basic material images are displayed in the same window.

Table 5

The means and STDs measurements of the ROIs indicated by the blue boxes within tissue images

Methods	FBP	PWLS-TV	PWLS-aviNLM	AlphaD-TV	AlphaD-aviNLM
Density	0.968±0.265	0.972±0.132	0.979±0.102	0.982±0.124	0.986±0.096

3.2 Physical phantom study

Figure. 7 shows the physical phantom results reconstructed by different algorithms from normal-and low- dose DECT measurements. The first row shows the physical phantom image reconstructed by the FBP algorithm from the normal dose scan, which serves as the ground truth image for evaluation. It can be seen that the SIR reconstructed results are of good quality. Specifically, the TV-based results are slightly smoother, especially in the vicinity of edges, than the aviNLM-based results. In addition, the presented AlphaD-aviNLM results achieve slight gain over the PWLS-aviNLM results in visual inspection. To further demonstrate the performance of the presented AlphaD-aviNLM algorithm, Table 6 lists the FSIM measurements of different algorithms for DECT images. It can be observed that the presented AlphaD-aviNLM algorithm is better than the other four algorithms, indicating that the presented AlphaD-aviNLM algorithm can obtain the best overall DECT image quality among the five algorithms.

Fig.7

The DECT images reconstructed by the FBP (second row), PWLS-TV (third row), PWLS-aviNLM (fourth row), AlphaD-TV (fifth row), and presented AlphaD-aviNLM (sixth row) algorithms from low-dose measurements. All the images are displayed in the same window.

Table 6

The FSIM measurements of the different SIR algorithms for DECT images

Methods	FBP	PWLS-TV	PWLS-aviNLM	AlphaD-TV	AlphaD-aviNLM
120 kVp	0.9345	0.9891	0.9947	0.9938	0.9956
90 kVp	0.9547	0.9933	0.9974	0.9965	0.9978

Figure. 8 shows the decomposed images reconstructed by different algorithms where the Teflon and water are selected as the basis materials. From the results, we have following observations: (1) Direct decomposition on the FBP reconstructed results in excessive noise in the decomposed images; (2) All the SIR algorithms yield better decomposed images quality than the FBP algorithm, indicating that the SIR algorithms can suppress noise-induced artifacts; (3) The aviNLM based algorithms performs better than the TV-based algorithms, especially in the regions indicated by the red arrows, implying that the patch-based reconstruction algorithms are superior to the pixel-based reconstruction algorithm. In addition, it can be observed that some shading artifacts occur in the basic material images, especially in the “water” images, which may be caused by the scattering effects. Table 7 summarizes the mean and STD measurements of ROIs indicated by the blue boxes as shown in Fig. 8. From the table, it can be observed that the presented AlphaD-aviNLM algorithm obtains mean values of the basic materials with negligible difference from the normal dose ones, and the smallest noise STDs on the basic materials among the SIR algorithms, respectively.

Fig.8

The Teflon and water images reconstructed by different algorithms. All the basic material images are displayed in the same window.

Table 7

The mean and STD measurements of ROIs indicated by the blue boxes as shown in Fig. 8

Methods	Normal dose	PWLS-TV	PWLS-aviNLM	AlphaD-TV	AlphaD-aviNLM
Teflon	1.864±0.114	1.859±0.184	1.862±0.155	1.861±0.129	1.863±0.117
Water	0.996±0.104	0.989±0.193	0.993±0.156	0.992±0.132	0.995±0.125

3.3 Real patient data study

To further demonstrate the performance of the presented AlphaD-aviNLM algorithm, Fig. 9 shows the patient 1 images reconstructed by the FBP, PWLS-TV, PWLS-aviNLM, AlphaD-TV, and presented AlphaD-aviNLM algorithms from the low-dose DECT sinogram data, respectively. As visualized in the results, the TV-based SIR algorithms can obtain a non-uniform intensity distribution in the homogeneous area in the reconstructed DECT images, but smooth out some detail as indicated by the red arrows. On the contrary, the aviNLM-based SIR reconstructed images are less noisy and more realistic than the TV-based ones, especially in the pulmonary lobe as indicated by the yellow boxes. The presented AlphaD-aviNLM algorithm can produce slightly sharper features and higher spatial resolution results. Additionally, Table 8 lists the corresponding noise level of the two ROIs, and also demonstrates the AlphaD-aviNLM algorithm can yield high quality DECT images. In addition, the VMS images (65 keV) generated from the DECT images are shown in Fig. 10. It can be seen that the presented AlphaD-aviNLM algorithm can yield the noticeable gains over the other three SIR algorithms in terms of noise suppression in VMS images. The real patient data study exhibits the same trend as the digital XCAT phantom and physical phantom studies through visual inspection. To qualitatively demonstrate the performance of the presented AlphaD-aviNLM algorithm in a subjective manner, five experienced physicians with at least five years in CT imaging scored the VMS images. The physicians were asked to score the VMS images from 0 (worst) to 5 (best) in terms of image noise, artifacts, edge and structures and overall image quality. In the experiments, the VMS images reconstructed by the different algorithms were displayed on the screen randomly, which is completely blind to the physicians. Table 9 lists subjective scores of VMS images (65 keV) from different algorithms. It is evident that the VMS image reconstructed by the presented AlphaD-aviNLM algorithm generally has the highest scores, indicating the presented AlphaD-aviNLM algorithm outperforms the other four algorithms from the physicians point of view. Figure. 11 shows the patient 2 images reconstructed by the different algorithms, and we can observe that the presented AlphaD-aviNLM algorithm is superior to the other four algorithms on the reconstruction of the lung region.

Fig.9

The patient 1 images reconstructed by the FBP (first row), PWLS-TV (second row), PWLS-aviNLM (third row), AlphaD-TV (fourth row), and presented AlphaD-aviNLM (fifth row) algorithms from low-dose measurements. All the images are displayed in the same window.

Fig.10

The virtual monochromatic spectral images (65 keV) generated from the DECT images. (a) is from the FBP algorithm; (b) is from the PWLS-TV algorithm; (c) is from the PWLS-aviNLM algorithm; (d) is from the AlphaD-TV algorithm; and (e) is from the presented AlphaD-aviNLM algorithm. All the images are displayed in the same window.

Fig.11

The patient 2 images reconstructed by the FBP, PWLS-TV, PWLS-aviNLM, AlphaD-TV and presented AlphaD-aviNLM algorithms from low-dose measurements. All the images are displayed in the same window.

Table 8

Material types in the digital XCAT phantom shown in Fig. 10

Methods	ROI1	ROI2
PWLS-TV	0.1103	0.1091
PWLS-aviNLM	0.0987	0.0952
AlphaD-TV	0.1009	0.0916
AlphaD-aviNLM	0.0893	0.0822

Table 9

Subjective scores of VMS images (65 keV) from different algorithms

Methods	FBP	PWLS-TV	PWLS-aviNLM	AlphaD-TV	AlphaD-aviNLM
Scores	3.032±0.652	3.538±0.564	3.902±0.539	3.865±0.588	4.053±0.524

4 Discussion and conclusion

In this study, we proposed and validated a SIR reconstruction algorithm for low-dose DECT data acquisitions. One motivation of this study is that the data fidelity term in the objection function is characterized by the Alpha-divergence criterion which is different from the conventional ML and MAP criteria. Similar to the ML and MAP criteria, the Alpha-divergence criterion is also used to measure the discrepancy between probability distributions [21]. In particular, the WLS criterion can be represented by the Alpha-divergence criterion in the case of α = 2 as well as α = -1 [21]. Therefore, the Alpha-divergence criterion can be of broad prospects in CT imaging applications. Another motivation of this study is that a NLM-based regularization term is introduced. Unlike the conventional quadratic and TV regularizations, the NLM-based regularizations adequately incorporate image nonolocal self-similarity, which is based on the fact that a nonlocal patch often has many nonlocal similar patches to it across the image [8, 27]. Many studies have proved that NLM-based regularizations have achieved a great success in image processing, especially in CT imaging [8 , 27]. Considering strong similarity between DECT images [8, 11], in this study, the aviNLM regularization utilizes data redundant information within DECT images. Moreover, the aviNLM regularization takes into account both the geometrical similarity and spatial distance while the TV regularization is only restricted by the spatial distance. This is also the reason why the patch-based models outperform the pixel-based models. The experiments were conducted with the digital phantom, physical phantom and real patient data. Qualitative and quantitative results in Section 3 demonstrated the feasibility and efficacy of the presented AlphaD-aviNLM algorithm in terms of quality-measure-utility metrics used.

The study suggests a number of interesting points meriting future study. The first is the parameter tuning. For the presented AlphaD-aviNLM algorithm, three main parameters need to be selected manually, namely the parameter p and the hyper-parameters β_h and β_l. In the study, the sizes of patch-window and patch-window do not show noticeable effects on the reconstructed DECT images when they are set in a reasonable range. The parameter p is determined through a broad range of values manually in terms of visual inspection and quantitative measurements. Once the parameter p is fixed, the hyper-parameters β_h and β_l are tuned manually to achieve a good compromise between a low noise level and a high level of detail. It should be argued that it would have been possible to yield better images using the PWLS-aviNLM algorithm than using the presented AlphaD-aviNLM algorithm with different parameters. Because there are indeed many parameters in the four SIR algorithms, we cannot guarantee that the ones we chose are the most suitable. Optimal parameters selection is the subject of future investigation. Second, the mismatch between the DECT images should be taken into account. In the physical phantom experiment, the physical phantom was fixed to be scanned with a single-source CT scanner. Therefore the consecutive acquisition DECT data can be obtained by using sequential scanning. The motion artifacts can be of a minimal problem. However, in the clinic, the slow-kVp DECT negatively impact the accuracy and precision of DE analysis, particularly when imaging moving organs. Therefore, to obtain accurate DE analysis, the misalignments due to motion or other effects should be minimized before the DECT reconstruction, which would be another topic in our future research plan. Third, the data were required from limited sample of patients and potential selection bias is unknown. Review of a clinical study with a variety of different lesion types would be beneficial to confirm if the presented AlphaD-aviNLM algorithms can be extended to a broader population, which is a worthy of future investigation. Finally, in this study, only the simplest basic material decomposition algorithm, i.e., direct matrix inversion operation, is utilized to obtain final basic material image. To further improve the performance, the statistical material decomposition algorithms can be incorporated into the presented AlphaD algorithm framework to obtain better basic material images, which are inspiring and worth investigation in future work.

In this work, our efforts were focused on the robust DECT images reconstruction via divergence constrained spectral redundancy information. In clinics, the presented algorithm can be used not only in DECT reconstruction but also in multi-energy CT imaging [24], perfusion imaging [27] and 4D CT imaging [25], which may be one of our interests for future research.

Footnotes

Acknowledgments

This work was supported in part by the China Postdoctoral Science Foundation funded project under Grants No. 2016M602489, the National Natural Science Foundation of China under Grant Nos. 81371544, 61571214, and 81501466, the Guangdong Natural Science Foundation under Grant Nos. 2015A030313271, and 2015A030310018, the Science and Technology Program of Guangdong, China under Grant No. 2015B020233008, the Science and Technology Program of Guangzhou, China under Grant No. 201510010039.

References

Heismann

B.J.

, Schmidt

B.T.

and Flohr

, Spectral Computed Tomography, Bellingham, WA: SPIE, (2012).

Silva

A.C.

, Morse

B.G.

, Hara

A.K.

, Paden

R.G.

, Hongo

and Pavlicek

, Dual-energy (spectral) CT: Applications in abdominal imaging, Radiographics 31 (2011), 1031–1046.

Mccollough

C.H.

, Leng

, Yu

and Fletcher

J.G.

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

Noh

, Fessler

J.A.

and Kinahan

P.E.

, Statistical sinogram restoration in dual-energy CT for PET attenuation correction, IEEE Trans Med Imaging 28 (2009), 1688–1702.

Harms

, Wang

, Petrongolo

, Niu

and Zhu

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

Petrongolo

and Zhu

, Noise suppression for dual-energy CT through entropy minimization, IEEE Trans Med Imaging 34 (2015), 2286–2297.

Zeng

, Huang

, Bian

, Zhang

, Niu

, Fan

, Liang

and Ma

, Dual energy CT images restoration via dictionary learning and spectral gradient Modeling, in Fully Three-Dimensional Image Reconstruction in Radiology and Nuclear Medicine Newport (2015), 659–662.

Zeng

, Huang

, Zhang

, Bian

, Niu

, Zhang

, Feng

, Chen

and Ma

, Spectral CT image restoration via an average image-induced nonlocal means filter, IEEE Trans Biomed Eng 63 (2016), 1044–1057.

Zhao

, Ding

, Lu

, Wang

, Zhao

and Molloi

, Dual-dictionary learning-based iterative image reconstruction for spectral computed tomography application, Phys Med Biol 57 (2012), 8217–8229.

10.

Zeng

, Huang

, Bian

, Zhang

, Niu

, Fan

, Liang

and Ma

, Dual energy CT image iterative reconstruction using an average image induced edge-preserving prior, in Fully Three-Dimensional Image Reconstruction in Radiology and Nuclear Medicine Newport (2015), 610–613.

11.

Szczykutowicz

T.P.

and Chen

G.H.

, Dual energy CT using slow kVp switching acquisition and prior image constrained compressed sensing, Phys Med Biol 55 (2010), 6411–6429.

12.

Zhang

, Thibault

J.B.

, Bouman

, Sauer

and Hsieh

, Model-based iterative reconstruction for dual-energy X-ray CT using a joint quadratic likelihood model, IEEE Trans Med Imaging 33 (2014), 117–134.

13.

Cai

, Rodet

, Legoupil

and Mohammaddjafari

, A full-spectral Bayesian reconstruction approach based on the material decomposition model applied in dual-energy computed tomography, Med Phys 40 (2013), 111916.

14.

Niu

, Gao

, Bian

, Huang

, Chen

, Yu

, Liang

and Ma

, Sparse-view X-ray CT reconstruction via total generalized variation regularization, Phys Med Biol 59 (2014), 2297–3017.

15.

Gao

, Bian

, Huang

, Zhang

, Niu

, Feng

, Chen

, Liang

and Ma

, Low-dose x-ray computed tomography image reconstruction with a combined low-mas and sparse-view protocol, Optics Express 22 (2014), 15190–15210.

16.

Zhang

, Han

, Liang

, Hu

, Liu

, Moore

, Ma

and Lu

, Extracting information from previous full-dose ct scan for knowledge-based bayesian reconstruction of current low-dose ct images, IEEE Transactions on Medical Imaging 35 (2016), 860–870.

17.

Zhang

, Ouyang

, Ma

, Huang

, Chen

and Wang

, Noise correlation in CBCT projection data and its application for noise reduction in low-dose CBCT, Medical Physics 41 (2014), 031906.

18.

Zhang

, Wang

, Ma

, Lu

and Liang

, Statistical models and regularization strategies in statistical image reconstruction of low-dose X-ray CT: A survey, Läkartidningen 108 (2014), 2660–2661.

19.

La Riviere

P.J.

, Penalized-likelihood sinogram smoothing for low-dose CT, Med Phys 32 (2005), 1676–1683.

20.

, Liang

, Fan

, Liu

, Huang

, Chen

and Lu

, Variance analysis of X-ray CT sinograms in the presence of electronic noise background, Med Phys 39 (2012), 4051–4065.

21.

, Liang

, Fan

, Liu

, Huang

, Lu

and Chen

, A study on CT sinogram statistical distribution by information divergence theory, In IEEE Nuclear Science Symposium and Medical Imaging Conference (2011), 3191–3196.

22.

Dong

, Niu

and Zhu

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

23.

Rigieand

D.S.

and La Riviere

P.J.

, Joint reconstruction of multi-channel, spectral CT data via constrained total nuclear variation minimization, Phys Med Biol 60 (2015), 1741–1762.

24.

Zeng

, Gao

, Huang

, Bian

, Zhang

, Lu

and Ma

, 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.

Zhang

, Ouyang

, Huang

, Ma

, Chen

and Wang

, Few-view cone-beam CT reconstruction with deformed prior image, Med Phys 41 (2014), 121905.

26.

Zhang

, Mou

, Wang

and Yu

, Tensor-based dictionary learning for spectral CT reconstruction, IEEE Trans Med Imaging 36 (2016), 142–154.

27.

, Zhang

, Gao

, Huang

, Liang

, Feng

and Chen

, Iterative image reconstruction for cerebral perfusion CT using pre-contrast scan induced edge-preserving prior, Phys Med Biol 57 (2012), 7519–7542.

28.

Zhang

, Ma

, Wang

, Liu

, Han

, Lu

, Moore

and Liang

, Statistical image reconstruction for low-dose CT using nonlocal means-based regularization. Part II: An adaptive approach, Computerized Medical Imaging and Graphics 43 (2015), 26–35.

29.

Huang

, Ma

, Liu

, Feng

and Chen

, Projection data restoration guided non-local means for low-dose computed tomography reconstruction, IEEE International Symposium on Biomedical Imaging: From Nano To Macro (2011), 1167–1170.

30.

Tian

, Ma

, Liang

, Huang

and Chen

, Information divergence constrained total variation minimization for positron emission tomography image reconstruction, In IEEE Nuclear Science Symposium and Medical Imaging Conference (2011), 2587–2592.

31.

Combettes

P.L.

and Wajs

V.R.

, Signal recovery by proximal forward-backward splitting, Siam Journal on Multiscale Modeling & Simulation 4 (2006), 1168–1200.

32.

Minka

, Divergence measures and message passing, Microsoft Research Ltd (2005), MSR-TR-2005-173.

33.

Segars

W.P.

, Sturgeon

G.S.

, Grimes

and Tsui

B.M.

, 4D XCAT phantom for multimodality imaging research, Med Phys 37 (2010), 4902–4915.

34.

Dong

, Huang

, Bian

, Niu

, Zhang

, Feng

, Liang

and Ma

, A simple low-dose X-ray CT simulation from high-dose scan, IEEE Trans Nucl Sci 62 (2015), 2226–2233.

Iterative reconstruction for low dose dual energy CT using information-divergence constrained spectral redundancy information

Abstract

Keywords

1 Introduction

2 Methods and materials

2.1 Statistical model in DECT sinogram domain

2.4.2 Selection of the parameters in the aviNLM model

2.4.3 Selection of the parameters β ⇀

2.5 Experimental data acquisition

2.5.1 XCAT phantom data

2.5.3 Patient data

2.6 Comparison Methods

3.1 XCAT Phantom Study

Footnotes

Acknowledgments

References

2.4.3 Selection of the parameters $\overset{⇀}{β}$