1.

INTRODUCTION

Recently, photon-counting computed tomography (PCCT) has been developed in clinics. It utilizes the photon counting detectors (PCDs) to simultaneously count photons in multi-energy bins. By obtaining abundant energydependent material-specific information and high contrast-to-noise ratio (CNR) for soft materials, PCCT has the advantages on material decomposition and lesion diagnostic.^{1, 2}

However, the narrow energy bin receives decreasing photons and the spectral data suffers from serious quantum noise.³ This challenge significantly degrades the quality of the reconstructed PCCT images in multi-energy bins and effects the imaging performance of PCCT imaging. To obtain high quality PCCT images, several statistical iteration reconstruction (SIR) methods have been exploited. It incorporates the statistical property of X-ray photons and prior information of the desired PCCT images to build reconstruction model with data fidelity and regularization terms, respectively. For examples, Xu et al. introduced a statistical interior tomography method with TV regularization to reconstruction PCCT images,⁴ Kim et al. stacked similar image patches at the same position and utilized a low-rank regularization to suppress image noise.⁵ Zhang et al. utilized a nonlocal mean regularization of the full-spectrum image to maintain the image details.⁶ Yao et al. proposed to improve the reconstruction performance by utilizing a nonlocal spectral similarity of a weighted image.⁷ Liu et al. developed a nonlocal total variation (TV) regularization term to construct the weights from the full-spectrum images.⁸ Tao et al. utilized the structural redundancy between the base materials and the spectral images to establish a prior-knowledge-aware material decomposition method.⁹ Zeng et al. analyzed the intrinsic tensor properties of the PCCT images and constructed a full-spectrum-knowledge-aware tensor model for PCCT imaging.¹⁰ These methods have been shown great potential in preserving image structures and suppressing noise. Moreover, deep learning (DL) based methods have been utilized in spectral CT. For examples, Li et al. constructed a cascade DNN to estimate the high-energy image from low-energy image,¹¹ and Cong et al. estimated the PCCT images from the data of energy-integrating detectors by using DNN based method.¹²

However, most of the current SIR methods assume the noise in PCCT images follows a Gaussian distribution, which deviates the real distribution of the noise in PCCT images. Because the artifacts, which maybe induced by beam hardening effect in lower energy bins or photon starvation effect of high density materials, would complicate the noise distribution in the image and invalidate the performances of the aforementioned SIR methods. In addition, the DL-based methods need a quantity of paired data to train a desired network, and the collection of training data is time-consuming and the clinical PCCT data is hard to be obtained. To solve it, we propose to utilize Gaussian Mixture Model (GMM),¹³ which is a universal approximator for any continuous density function, to model the noise in the PCCT images. Moreover, we treat the PCCT images as a 3-order tensor and serve the three dimensional total variation (3DTV)¹⁴ as image prior. Finally, we construct a statistical iteration reconstruction method based on GMM and 3DTV. For shorten, we call the presented reconstruction method as “GMM-3DTV”.

In summary, the main contributions of this work are:

2.

METHODS

PCCT receives spectral data among multi-energy bins. Considering the spatial and energy dimensions of the measurements, the 3D PCCT imaging model can be expressed as follows:

where y = {y_n, n ≤ N} and X = {x_n, n ≤ N} denote the measurements and reconstructed PCCT images among multi-energy bins, respectively. N is the total number of the multi-energy bins. A is a linear projection operator. It should be note that additional constraints are incorporated to stable the solution of the above model.

Due to the variance of the material attenuations for X-ray along the whole spectrum, the feature of the noise is different among multi-energy bins. Specifically, the reconstructed images in lower bins suffer from strip artifacts due to the beam hardening, and the ones in higher bins are corrupted by the photon starvation induced noise. Therefore, a single Gaussian is no sufficient to approximate the complex noise distribution in the PCCT images.

Fig. 1 shows the noisy PCCT images, noise and ground truth. It can be seen that the PCCT images have multiple modalities of noise, as illustrated by the three Gaussian components in Bin 1, and a simple Gaussian noise model may deviate the real cases. In order to solve the complex noise distribution, we model the noise in the PCCT images as a parametric probability distribution by GMM for more flexibly adapting different cases. Specifically, the noise term ε_n in the nth energy bin is expressed as follows:

Figure 1.

Illustration of the approximated results by the GMM for complex noise in multi-energy bin PCCT images. First row (from left to right): Noisy PCCT images, noise and ground truth. Second row: three noise components of GMM in Bin 1.

where π_nk, μ_nk and denote the mixture rate, mean and variance values of the kth Gaussian compound in nth energy bin, respectively. In this work, the mean values are set to be zero. K is the total number of Gaussian compounds, and .

Considering the sparsity structures and low rank property of the PCCT images, a general 3D total variation (3DTV) term is utilized as the prior, as follows:

where ∇_m is different operations along spatial height, width and spectrum modes of X, || · ||₁ is the L₁ norm. Therefore, coupling the GMM approximation for noise and sparse regularization for PCCT images, we can formulate a robust penalized weighted least squares method for PCCT imaging as follows:

where α is a hyper-parameter of the image prior term. Simply, we call the above method as “GMM-3DTV” method. Finally, by imposing the negative form of likelihood function of the GMM, Eq. (4) is rewritten as follows:

where β is the hyper-parameters of GMM likelihood term. Π = {π_nk n ≤ N,k ≤ K} and K} are the sets of mixture coefficients and variance values, respectively. Moreover, we adopt an Expectation Maximization (EM) algorithm ¹⁵ to iteratively optimize the presented GMM-3DTV method.

The whole EM algorithm for optimizing Eq. (5) can be summarized in Algorithm 1.

3.

RESULTS

3.1

Implementation Details

In this work, we compare the presented GMM-3DTV method with the filtered back projection (FBP) method using a ramp filter and a tensor-based dictionary learning regularization (TDL) method. In addition, two more methods are implemented to investigate different parts of the presented GMM-3DTV method. The one is GMM-based reconstruction model, and the other is 3DTV-based reconstruction model. The images at normal dose serve as the ground truth. Two numerical phantoms are utilized to validate the performance of the presented method. Specifically, XCAT phantom¹⁶ contains alcohol, water, bone and gadolinium, and the synthesized clinical phantom includes water, bone and iodine. The phantoms are scanned under the simulated 120 kVp X-ray spectrum with 1.6 mm Al filtration by SPEKTR toolbox.¹⁷ Five energy bins with equal photon counts are determined by the thresholds: 25, 50, 60, 70 and 85 keV. The PCCT imaging parameters are set as follows: (1) source-to-detector and source-to-center distances are 1040.0 and 570.0 mm, respectively. (2) 1160 projection views are evenly scanned the objections over 360^o. (3) 816 detector channels are placed along the parallel X-ray beam. To generate the noisy projections, Poisson noise is applied into the simulated noise-free projections.

3.2

Qualitative Analysis

Fig. 2 shows visual comparisons of the presented and compared methods on XCAT phantom. It can be observed that: 1) the FBP method suffers from noise-induced artifacts, and the other methods outperform the FBP method in terms of improving the image quality; 2) the 3DTV-based method suffers from blocky artifacts and the GMM-based method remains producing noisy images; 3) the TDL method smooths the images and enhances the texture of the noise-induced artifacts; 4) the presented GMM-3DTV method better handles the noise-induced artifacts and preserves the image details. Moreover, zoomed in regions-of-interest (ROIs) indicated by the red boxes in each image are illustrated. It can be observed that the presented GMM-3DTV maintains the details of the anatomic structures.

Figure 2.

Reconstructed images of the presented and compared methods on XCAT phantom. The display windows from Bin 1 to 5 are [0.00, 0.02], [0.00, 0.010], [0.00, 0.017], [0.00, 0.015] and [0.00, 0.010] mm^–1, respectively. Zoomed ROIs indicated by the red boxes are displayed for better visualization.

Fig. 3 illustrates the results of different methods on synthesized clinical phantom. Similar with the results on XCAT phantom, the TDL method fails to denoise among the multi-energy bins, the 3DTV induces additional block artifacts for the denoising results, and the GMM-based method hardly removes the noise. On the contrary, the presented GMM-3DTV method can effectively remove the noise-induced artifacts and recover the images details. For better visual inspection, zoomed in ROIs indicated by the blue boxes are shown, which also demonstrates the same conclusion.

Figure 3.

Reconstructed images of the presented and compared methods on synthesized clinical phantom. The display windows from Bin 1 to 5 are [0.007, 0.0105], [0.0032, 0.0072], [0.003, 0.004], [0.0022, 0.0035] and [0.0018, 0.0027] mm^–1, respectively. Zoomed ROIs indicated by the blue boxes are displayed for better visualization.

Fig. 4 shows the profiles of different methods on synthesized clinical phantom indicated by the green line in Fig 3. From the results, we can seen that the TDL and 3DTV-based methods produce blurry results, and the GMM-based method fails to suppress the noise. In contrast, the results of the presented GMM-3DTV method are the closest to the ground truth.

Figure 4.

Profiles of the results from different methods on synthesized clinical phantom. The location of the profiles is indicated by the green line in Fig. 3 for each reconstructed images in Bin 1.

4.

DISCUSSION AND CONCLUSION

Due to the beam-hardening and photon starvation effects, the reconstructed PCCT images suffer from noise-induced artifacts with complex noise distribution in image domain. Most of the SIR methods hardly handle the artifacts and produce suboptimal results. To address this issue, in this work, we presented a novel SIR PCCT reconstruction method. Specifically, we utilized GMM to approximate the complex noise distribution in the PCCT image domain. Moreover, a 3DTV term, which served as the image prior, was also incorporated into the reconstruction model to encourage structural similarity of the PCCT images along the multi-energy bins. Experiments were conducted to demonstrate the effectiveness of the presented GMM-3DTV method. In the future, clinical and more scene studies would be included to further demonstrate the reconstruction performance of the presented GMM-3DTV method.