Projection data smoothing through noise-level weighted total variation regularization for low-dose computed tomography

Abstract

Reducing radiation dose while maintaining the quality of the reconstructed images is a major challenge in the computed tomography (CT) community. In light of the non-stationary Gaussian noise distribution, we developed a model that incorporates a noise-level weighted total variation (NWTV) regularization term for denoising the projection data. Contrary to the well-known edge-weighted total variation method, which aims for better edge preserving, the proposed NWTV tries to adapt the regularization with the spatially varying noise levels. Experiments on simulated data as well as the real imaging data suggest that the proposed NWTV regularization could achieve quite competitive results. For sinograms with sharp edges, the NWTV could do a better job at balancing noise reduction and edge preserving, such that noise is removed in a more uniform manner. Another conclusion from our experiments is that the well-recognized stair-casing artifacts of TV regularization play little role in the reconstructed images when the NWTV method is applied to low-dose CT imaging data.

Keywords

Denoising weighted total variation low-dose ct non-stationary gaussian noise

1 Introduction

Currently, X-ray computed tomography (CT) is widely applied in clinical diagnosis and industry for nondestructive testing. One particular concern regarding to its application is the X-ray dose. In medical applications, the X-ray radiation may increase the risk of somatic damage [1, 2], which calls for reducing the radiation dose. In industrial CT and Micro-CT applications, low-dose X-ray projection data are commonly seen due to efficiency requirement or in time critical applications like online test. Another scenario is that when scanning large or high-density object, the counts of detected photons would be low due to heavy attenuation. Low-dose leads to increased level of noise in the projection data, and the images reconstructed via conventional analytical reconstruction methods will suffer from obvious noise and artifacts. Thus, low-dose CT has always been an active research field in the CT community.

Many techniques have been developed to improve the quality of low-dose reconstructions, which can be roughly divided into three categories: (A) projection data pre-processing approaches; (B) iterative reconstruction with proper regularizations; and (C) post-processing approaches. The methods in category (A) try to estimate the noise-free projection data based on some statistical properties of the noisy data and then do reconstruction from the recovered sinogram by the filtered back projection (FBP) algorithm [3]. In [4], the Rudin-Osher-Ratemi (ROF) model [5] was modified to deal with Poisson noise and applied to raw data smoothing. In [6], the Anscombe transform is applied on raw data and then the bilateral filtering method is utilized to smooth the resulted projection data. To explore local similarity, the patch-based dictionary learning approach was adopted for sinogram smoothing in [7]. According to experimental results, the calibrated and log-transformed projection data is believed to follow the so-called non-stationary Gaussian distribution [8]. Please also refer to [9] for a more complex model which considers additional electronic background noise. Based on this noise model, many methods have been proposed in the literature for projection data smoothing [8 , 11]. An obvious advantage of the methods in this category is that both the denoising and reconstruction processes could be performed quite efficiently.

Iterative reconstruction methods in category (B) are gaining popularity due to their better noise suppression property and the increasing computation power. Two iterative reconstruction methods are commonly used: statistical iterative reconstruction and (simultaneous) algebraic reconstruction (ART or SART) [12, 13]. To better suppress noise, ART and SART are usually combined with regularizers, such as total variation (TV) [14, 15], Non-local means (NLM) filtering [16, 17] and sparsity constraint based on dictionary learning [18, 19]. Iterative reconstruction can produce better CT images, at the cost of higher computational burden.

In contrast to the first two categories, the post-processing methods in category (C) can be directly applied to the reconstructed CT images, and easily integrated into the current CT workflow. Traditional methods include NLM [20, 21] and block-matching 3D (BM3D) [22, 23] as well as dictionary learning based algorithms [24, 25]. Deep learning techniques have also been applied for CT image processing which could achieve very competitive results [26, 27].

Generally speaking, to achieve better result, the noise properties shall be given fully consideration and properly incorporated into the workflow. In this regard, the post-processing methods in Category (C) suffer a common difficulty that the statistical properties of the noise in the reconstructed images are very complex to analysis, due to various data pre-processing procedures like system calibration and log-transformation, as well as inexact image reconstruction model. Some efforts dedicated to noise analysis in the image domain could be found in [21].

On the other hand, in the projection space, the noise model is relatively easier to establish. In fact, it has been verified that the projection data suffers from non-stationary Gaussian noise [8]. This noise model is based on the assumption that the raw data, i.e. the collected data by the detector before calibration and log-transform, follows the Poisson distribution. In fact, there exists studies on denoising the raw data directly [4]. However, many reconstruction algorithms are applied on the projection data rather than the raw data.

Except the noise model, the smoothness properties of the projection data should also be taken into consideration for successful denoising. The ideal image is usually modeled as a piecewise constant function, which lives in the space of bounded variation (BV). It has been shown that the projection data is half-order smoother than the ideal reconstructed image in terms of some Sobolev space [28]. Another smoothness characterization can be found in [29]. So, BV is not the proper space for modeling the projection data, and piecewise smooth might be a better choice This seems to rule out TV based denoising methods, such as the famous ROF model which favors piecewise constant approximations.

In this paper, however, we will show that TV based methods could still perform well for projection data denoising. We propose a weighted TV regularization term (NWTV), whose smoothing strength is proportional to the spatially varying noise levels of the projection data, i.e. the non-stationary Gaussian noise. As a regularization term, NWTV is designed to achieve two goals: edge-preserving smoothing and adaptive (to spatially varying noise levels) smoothing. A related idea can be found in [30], where the TV term incorporating adaptive (to local noise levels) weighting regularization parameters is proposed to remove Rician noise from MRI images. In the NWTV case, the adaptivity is built into the TV term and a constant regularization parameter is adopted. As stated before, TV might not be a good choice since it prefers piecewise constant approximations even for smooth signals, which leads to stair-casing artifacts, also known as blocky effect [31]. This concern is actually not very serious when TV is applied for projection data smoothing. The reason behind is that projection data smoothing, whether TV based or not, should be gentle, i.e. the regularization strength should be small. Otherwise, heavy streak artifacts shall be introduced. With dominating data term, the stair-casing artifacts of NWTV shall be reduced. On the other hand, even for denoising piecewise smooth images, edge-preserving is an essential requirement, which fits well to TV regularization.

The remainder of this paper is organized as follows. The NWTV model will be introduced in Section 2. The efficient split-Bregman algorithm for solving the NWTV model will be developed in Section 3. Section 4 dedicates to experiments with both simulated and real projection data, and conclusions and acknowledgments will be given in Section 5.

2 The NWTV regularization

As stated previously, the noise in the projection data follows approximately the non-stationary Gaussian distribution. To successfully suppress the noise, a denoising method should be aware of this noise model. In this section, we will introduce and verify the non-stationary Gaussian noise distribution for the projection data and then introduce the NWTV regularization term to fully utilize this distribution.

2.1 The non-stationary Gaussian noise model

As far as noise properties are concerned, the low-dose CT measurement data (raw data) has been proved to be consistent with the Poisson distribution [32], which is based on the fact that the number of detected photons emitted from X-ray source conforms to the Poisson distribution. It was shown in [33] that $σ_{i}^{2} \approx 1 / - {\bar{I}}_{i}$ , where ${\bar{I}}_{i}$ and $σ_{i}^{2}$ are the mean and variance of the ith raw data, i.e. the detected photon numbers I_i, respectively.

However, the CT projection data (sinogram) used for reconstruction is obtained by performing several operations on the raw data, such as system calibration and log-transformation. Li et al. [8] proposed a non-stationary Gaussian noise model for the projection data, which relates its mean and variance by

$σ_{i}^{2} (p_{i}) = l_{i} e^{α {\bar{p}}_{i}},$ (1) where p _i is the ith value of the sinogram, ${\bar{p}}_{i}$ denotes the mean of p_i, α and l_i are object-independent parameters, which are determined by the scanning system.

To calibrate the parameters l_i and α in Equation (1), we did a similar experiment as described in [8]. An aluminum sample was repeatedly scanned with an industrial CT system located in our lab. The geometrical settings of our system are: the distance from the X-ray source to the turntable center (SOD) is 700.00 mm, the distance between the X-ray source and the detector center (SDD) is 1060.78 mm and the length of each detector cell is 0.083 mm. The scanning was repeated 2500 times at one fixed angle with 350 kV and 0.55 mA. Then the raw data was transformed by the system calibration (background noise subtraction) and log-transformation to get the projection data. Figure 1(a) shows the obtained projection data. The relationship between the sample mean and variance is plotted in Fig. 1(b), while the log-variance vs. mean is plotted in Fig. 1(c), where the straight lines clearly show the exponential relationship between the mean and the variance.

Fig.1

Verification of the non-stationary Gaussian noise model. (a) the sinogram; (b) variance vs. mean; (c) log-variance vs. mean.

In our simulations, for simplicity, the parameter l_i was treated as a constant, and the mean and variance relationship was modeled as $σ_{i}^{2} (p_{i}) = l e^{α {\bar{p}}_{i}}$ . After taking the Log-transform, a linear equation $ln σ_{i} = ln l + α {\bar{p}}_{i}$ can be derived, from which ln l (so the parameter l) and α can be computed by fitting the data shown in Fig. 1(c) with a linear model. Based on our experimental data, the parameters l and α for our CT system were computed approximately as l = 1.6e– 4 and α= 1.1.

2.2 NWTV in light of non-stationary Gaussian noise

A denoising method or model usually assumes certain noise distributions. For instance, the ROF model introduced by Rudin, Osher and Fatemi in 1992 [5] assumes Gaussian noise, which has been proven to be very effective for suppressing noise while preserving edges. Other popular conventional denoising methods, such as NLM and BM3D, also assume Gaussian noise. Another easily seen noise model in the literature is Poisson distribution. Please refer to [4, 34].

Compared to methods like NLM and BM3D, optimization based denoising models enjoy the property that the solution, i.e. the denoised image could be well characterized in terms of some functional spaces. So, we would like to incorporate the noise property of the projection data into some energy functional, so as to build an optimization model to recover the noise-free projection data. This idea has been utilized in our previous work for removing Poisson noise for the raw data [4]. In this section, we will generalize this idea to deal with the projection data. This is meaningful since conventional CT reconstruction algorithms use directly the projection data rather than the raw data. So, if one wants to combine the denoising procedure into some reconstruction algorithm, it’s much easier to do it with the projection data.

We start from the ROF model to derive our NWTV model. To make our formulations easier, here we adopt a special notation for digital images. We would like to model a digital image as a 2-dimensional piecewise constant function rather than a discrete 2-dimensional vector. For a digital image of size M×N, its pixels are indexed by (x, y) , x ∈ [0, M] , y ∈ [0, N], such that (x₀, y₀) refers to the first pixel whenever x₀ ∈ (0, 1] and y₀ ∈ (0, 1]. When we say two pixels (x₁, y₁) and (x₂, y₂) are different, it actually means that $(floor (x_{1}), ceil (x_{1})] \cap (floor (x_{2}), ceil (x_{2})] = φ,$ or $(floor (y_{1}), ceil (y_{1})] \cap (floor (y_{2}), ceil (y_{2})] = φ,$ where the function floor takes as input a real number x and outputs the greatest integer less than or equal to x. Similarly, the function ceil maps x to the least integer greater than or equal to x. Define Ω to be the set that hosts all the pixels, i.e. Ω = (0, M] × (0, N]. By the above notations, formulations for continuous signals could be easily adapted to describe discrete digital images. Assuming that f is a noisy image with the Gaussian noise n, and u is the ideal noise-free image, we have f (x, y) = u (x, y) + n (x, y). For Gaussian noise, the digital image f (x, y) can be regarded as samplings of independent and identical random variables, so we have

$\int_{Ω} fd Ω = \int_{Ω} ud Ω, \frac{1}{| Ω |} 1 pt \int_{Ω} {(u - f)}^{2} d Ω = σ^{2},$ (2) and

$\frac{1}{| Ω |} 1 pt \int_{Ω} \frac{{(u - f)}^{2}}{σ^{2}} 1 pt 1 pt d Ω = 1 .$ (3)

For Gaussian noise, Equation (3) provides a data fidelity term that might be utilized in denoising models, like the ROF model $\hat{u} = arg min_{u} 1 pt \int_{Ω} | \nabla u | 1 pt d Ω 1 pt + \frac{μ}{2} 1 pt \int_{Ω} {(u - f)}^{2} d Ω,$ (4) where μ is a weighted parameter to balance the data term and the TV term, and Δ denotes the gradient operator. For continuous image u, $\nabla u = {(\frac{\partial u}{\partial x}, \frac{\partial u}{\partial y})}^{T}$ , where the superscript T denotes matrix transposition.

In the following, the Equation (3) will be generalized to deal with non-Gaussian noise models, e.g. the non-stationary Gaussian noise model. As mentioned previously, each pixel of the projection data p (x, y) can be regarded as a random variable whose variance and mean possess the exponential relationship. Let u (x, y) denote the mean of p (x, y) that needs to be recovered, and $σ_{p}^{2}$ denote the variance of p. Consider the random variable p’s standardization $\tilde{p} = (p - u) / - σ_{p}$ . It’s well known that for any $(x, y), \tilde{p} (x, y)$ follows a distribution with mean zero and variance one. As always done in the literature, we could further assume that any two variables p(x₁, y₁), p(x₂, y₂) should be independent of each other if (x₁, y₁) (x₂, y₂). If we further assume that ${\tilde{p} (x, y), (x, y) \in Ω}$ are actually following the same distribution (with mean zero and variance one), then we have

$\frac{1}{| Ω |} 1 pt \int_{Ω} {\tilde{p}}^{2} d Ω = 1, i . e . \frac{1}{| Ω |} 1 pt \int_{Ω} \frac{1}{σ_{p}^{2}} 1 pt {(u - p)}^{2} 1 pt d Ω = 1 .$ (5)

It should be pointed out that, the assumption for the standardization variables, i.e. ${\tilde{p} (x, y), (x, y) \in Ω}$ collectively follow the same distribution, can not be true in general. However, in the non-stationary Gaussian case, i.e. p follows the Gaussian distribution [8] $PDF (p) = \frac{1}{\sqrt{2 π σ_{p}^{2}}} 1 pt exp (- \frac{{(p - u)}^{2}}{2 σ_{p}^{2}}), σ_{p}^{2} = l e^{α u},$ (6)

Then it’s easy to verify that $\tilde{p}$ really follows the standard Gaussian distribution. Thus, we have

$\frac{1}{| Ω |} 1 pt \int_{Ω} \frac{1}{l e^{α u}} 1 pt {(u - p)}^{2} 1 pt d Ω = 1 .$ (7)

By assuming a prior distribution related to total variation, we arrive at the following optimization problem $\hat{u} = arg 1 pt min_{u} 1 pt \int_{Ω} | \nabla u | 1 pt d Ω + \frac{μ}{2} 1 pt 1 pt \int_{Ω} \frac{1}{e^{α u}} 1 pt {(u - p)}^{2} 1 pt d Ω .$ (8)

Here the system parameter l has been absorbed by the parameter μ.

We acknowledge that the above model could have been derived simply from the Bayes point of view. However, we think that the above derivation is valuable in the sense that, with proper approximations, it might still be applicable even the assumption with the standardization variables fails to hold. For example, in the Poisson distribution case, each of the standardization variables follows a different distribution [4], even though all of them have mean zero and variance one.

Compared to the ROF model (4), the data term is highly nonlinear due to the extra exponential function, which leads to difficulties in devising fast solvers. So as done in [35] for Poisson noise removal, we move the exponential function to the regularization term, and arrive at $\hat{u} = arg 1 pt min_{u} 1 pt \int_{Ω} e^{α u} 1 pt | \nabla u | 1 pt d Ω + \frac{μ}{2} 1 pt \int_{Ω} {(u - p)}^{2} 1 pt d Ω .$ (9)

The regularization term ∫_Ωe^αu1pt| ∇ u|1ptdΩ is the so called NWTV, which encodes the noise level as well as the gradient sparsity prior for the noise-free data. As explained before, the piecewise constant assumption might not be appropriate for the projection data. However, the TV regularization enjoys the well-known edge-preserving property, which is also needed for projection data smoothing. Another reason behind for employing the TV regularization is that the projection data smoothing must be gentle, i.e. one can only expect to reduce part of the noise with small regularization parameters, so as to avoid streak artifacts. That’s to say, the TV regularization should also be incorporated with a small weighting parameter, in which case, the data term dominates the energy and the stair-casing artifacts should be slight. The unique attribute of the NWTV is that besides edge preserving, it also adapts to spatially varying noise levels.

The NWTV model (9) corresponds to the isotropic ROF model. One can also consider its anisotropic version

$\hat{u} = arg 1 pt min_{u} 1 pt \int_{Ω} e^{α u} 1 pt | \nabla_{x} u | 1 pt d Ω + \int_{Ω} e^{α u} | \nabla_{y} u | 1 pt d Ω + \frac{μ}{2} 1 pt \int_{Ω} {(u - p)}^{2} 1 pt d Ω .$ (10)

Both models could be solved by the split-Bregman algorithm [36] efficiently. Since the formulations are quite similar, and in applications, both isotropic NWTV and anisotropic NWTV achieve similar results, we will consider only the anisotropic version.

3 Algorithm

The split-Bregman algorithm is a very efficient and popular solver for the ROF model. The basic idea is to introduce additional variables to split the TV term to remove its nonlinearity, and then apply the Bregman iteration to deal with the constraints. It has been proved that split-Bregman is equal to the augmented Lagrangian algorithm in certain situations [37], and it also has close connections to ADMM (Alternating Direction Method of Multipliers) [38, 39]. Here we adopt the same idea to generalize it to solve the NWTV model (10). Let us introduce variables d_x = ∇ _xu, d_y = ∇ _yu, and the relaxation variables b_x, b_y. In the discrete setting, the split-Bregman algorithm applied to our model (10) reads ${\begin{matrix} (u^{k + 1}, d_{x}^{k + 1}, d_{y}^{k + 1}) = arg 1 pt min_{u, d_{x}, d_{y}} 1 pt e^{α u} {∥ d_{x} ∥}_{1} + e^{α u} {∥ d_{y} ∥}_{1} + \frac{μ}{2} 1 pt {∥ u - p ∥}_{2}^{2} \\ + \frac{θ}{2} 1 pt {∥ u - u^{k} ∥}_{2}^{2} + \frac{λ}{2} 1 pt {∥ d_{x}^{k} - \nabla_{x} u - b_{x}^{k} ∥}_{2}^{2} + \frac{λ}{2} 1 pt {∥ d_{y}^{k} - \nabla_{y} u - b_{y}^{k} ∥}_{2}^{2} \\ b_{x}^{k + 1} = b_{x}^{k} + \nabla_{x} u^{k + 1} - d_{x}^{k + 1} \\ b_{y}^{k + 1} = b_{y}^{k} + \nabla_{y} u^{k + 1} - d_{y}^{k + 1} \end{matrix}$ (11)

Note that the new term $\frac{θ}{2} 1 pt {∥ u - u^{k} ∥}_{2}^{2}$ is added to increase the stability of the solution process of u. In addition, to further easy the optimization with respect to u, the term e^αu will be linearized and approximated by e^{αu
^k} during the iterations, i.e.

$\begin{matrix} (u^{k + 1}, d_{x}^{k + 1}, d_{y}^{k + 1}) = arg 1 pt min_{u, d_{x}, d_{y}} 1 pt e^{α u^{k}} {∥ d_{x} ∥}_{1} + e^{α u^{k}} {∥ d_{y} ∥}_{1} + \frac{μ}{2} 1 pt {∥ u - p ∥}_{2}^{2} \\ + \frac{θ}{2} 1 pt {∥ u - u^{k} ∥}_{2}^{2} + \frac{λ}{2} 1 pt {∥ d_{x}^{k} - \nabla_{x} u - b_{x}^{k} ∥}_{2}^{2} + \frac{λ}{2} 1 pt {∥ d_{y}^{k} - \nabla_{y} u - b_{y}^{k} ∥}_{2}^{2} . \end{matrix}$ (12)

The above optimization problem (12) can be further split into three simpler subproblems by the alternating optimization strategy, which read

u-subproblem

$\begin{matrix} u^{k + 1} = arg 1 pt min_{u} 1 pt \frac{μ}{2} 1 pt {∥ u - p ∥}_{2}^{2} + \frac{θ}{2} 1 pt {∥ u - u^{k} ∥}_{2}^{2} + \frac{λ}{2} 1 pt {∥ d_{x}^{k} - \nabla_{x} u - b_{x}^{k} ∥}_{2}^{2} \\ + \frac{λ}{2} 1 pt {∥ d_{y}^{k} - \nabla_{y} u - b_{y}^{k} ∥}_{2}^{2} \end{matrix}$ (13)d-subproblem ${\begin{matrix} d_{x}^{k + 1} = arg 1 pt min_{d_{x}} e^{α u^{k}} 1 pt {∥ d_{x} ∥}_{1} + \frac{λ}{2} 1 pt {∥ d_{x}^{k} - \nabla_{x} u - b_{x}^{k} ∥}_{2}^{2} \\ d_{y}^{k + 1} = arg 1 pt min_{d_{y}} e^{α u^{k}} 1 pt {∥ d_{y} ∥}_{1} + \frac{λ}{2} 1 pt {∥ d_{y}^{k} - \nabla_{y} u - b_{y}^{k} ∥}_{2}^{2} \end{matrix}$ (14)

For the u-subproblem, its optimality condition requires $[(μ + θ) I - λ Δ] u^{k + 1} = μ p + θ u^{k} + λ \nabla_{x}^{T} 1 pt (d_{x}^{k} - b_{x}^{k}) + λ \nabla_{y}^{T} 1 pt (d_{y}^{k} - b_{y}^{k}) .$ (15)

The above equation could have been solved by applying FFT efficiently, at the cost of possible boundary artifacts due to the periodic boundary conditions presumed by FFT. To avoid this problem, we use the Gauss-Seidel algorithm to solve u.

$\begin{matrix} G_{i, j}^{k} = u_{i, j}^{k + 1} = \frac{λ}{μ + θ + 4 λ} 1 pt u_{i + 1, j}^{k} + u_{i - 1, j}^{k} + u_{i, j + 1}^{k} + u_{i, j - 1}^{k} + d_{x, i, j - 1}^{k} \\ - d_{x, i, j}^{k} + d_{y, i, j - 1}^{k} - d_{y, i, j}^{k} - b_{x, i - 1, j}^{k} + b_{x, i, j}^{k} - b_{y, i, j - 1}^{k} \\ + b_{y, i, j + 1}^{k} + \frac{λ}{μ + θ + 4 λ} 1 pt p_{i, j} + \frac{λ}{μ + θ + 4 λ} 1 pt u_{i, j}^{k}, \end{matrix}$ (16) where i, j refer to the pixel (x_i, y_j). The d-subproblems could be solved very efficiently by the shrinkage operator, i.e.

${\begin{matrix} d_{x}^{k + 1} = shrink (\nabla_{x} u^{k + 1} + b_{x}^{k}, \frac{e^{α u^{k + 1}}}{λ}), \\ d_{x}^{k + 1} = shrink (\nabla_{y} u^{k + 1} + b_{y}^{k}, \frac{e^{α u^{k + 1}}}{λ}), \end{matrix}$ (17) where the function shrink is defined as $shrink (x, a) = {\begin{matrix} x - a & x \in (a, + \infty), \\ 0 & x \in [- a, + a], \\ x + a & x \in (- \infty, - a), \end{matrix}$ and when the inputs are vector or matrices, it performs its operation in an elementwise manner.

The whole algorithm to solve the NWTV model can be summarized as Algorithm 1.

4 Experiment

To verify the effectiveness of the NWTV model and its fast solution algorithm, experiments against popular algorithms specialized for CT image denoising will be performed on both synthesized images and real data. Two popular algorithms: the penalized-weighted least-squares (KL_PWLS) algorithm [10] and adaptive non-local means (ANLM) filter [21] are selected as the competing algorithms. Comparisons will be performed on both the projection domain and the CT image domain.

In light of the correlations between the neighboring views, the KL_PWLS algorithm transforms the sinograms into the KL domain (i.e. Karhunen- Loéve transform domain, for details, see [10]), which are then denoised and transformed back to the projection domain. The final CT images are then reconstructed by the FBP algorithm. The ANLM filter works on the image domain which incorporates the spatially varying noise levels of the reconstructed images into constructing the weights of the non-local means algorithm. The noise levels are estimated from the projection data by tracking and analyzing the noise propagation behavior during the reconstruction procedure.

In the projection domain, the NWTV will be compared with the KL_PWLS algorithm in terms of the quality of the recovered sinograms. In the CT image domain, the quality of the reconstructed images (by the FBP algorithm) from the denoised sinograms will be compared.

Postprocessing denoising algorithms could be still applied on the reconstructed images for further improvement. In this regard, the sinogram restoration induced non-local means (SR_NLM) filter [20] provides a strategy to combine the reconstructed images with the non-local means algorithm. So, we will propose another algorithm, named NWTV_NLM. Following the idea of SR_NLM, NWTV_NLM will first denoise the sinogram by the NWTV algorithm, then images are reconstructed by the FBP algorithm. The last step is to apply the non-local means filter on the reconstructed images with weighting parameters computed by taking the reconstructed image as mean values. To compare NWTV and KL_PWLS further, one could easily construct another algorithm, named KLPWLS_NLM, which just replaces the NWTV denoising part with the KL_PWLS algorithm. The ANLM filter will also be implemented as one of the competing algorithms.

To draw some conclusions in a quantitative way, quality measures including the peak signal to noise ratio (PSNR) [40], normal root mean square error (NRMSE), structural similarity (SSIM) [41] and universal quality index (UQI) [42] are computed and compared respectively. The simplest and most widely used quality index is the PSNR which reflects the noise level in the signal, and higher PSNR value generally indicates higher quality of the recovered images. NRMSE shows the divergence of two images with regard to the human vision, and a lower NRMSE value generally means better similarity. The SSIM index is to measure the structural similarity of two images. Compared to PSNR, the UQI is a better measure for image distortions like loss of correlation, luminance distortion and contrast distortion.

4.1 Validation of the NWTV regularization

In this subsection, we will demonstrate how the NWTV algorithm performs on balancing the spatially varying noise levels of the non-stationary Gaussian distribution. The test image consists of piecewise constant and piecewise smooth blocks of different shapes. The gray value of the image ranges from 0 to 3. According to the relationship between noise variance and mean in Equation (1), we first calculate the noise variances with the original image as the means. Then a random field is generated according to the Gaussian distribution with the calculated variances, which acts as the noisy image. The added noise (non-stationary Gaussian noise) can be identified as the difference image between the noisy image and the mean image (ideal image). Note that in all the experiments, when needed, noisy image is always generated in the same way described above. The non-stationary Gaussian noise with l = 0.001, α = 2 is added to the test image. The ideal image, noisy image and the added noise are shown in Fig. 2(a), (b) and (c), respectively. It can be clearly seen from Fig. 2(c) that, the noise level of each part is different, and higher gray value shows higher noise level.

Fig.2

(a) the ideal image; (b) the noisy image; (c) the added non-stationary Gaussian noise.

Four denoising methods including NLM, BM3D, TV (i.e. ROF model) and the proposed NWTV are used to process the noisy image. The parameters of each method have been tuned to achieve best balance between noise removal (high PSNR values) and edge-preserving by sampling the parameter space. The denoising results are shown in Fig. 3. The line profiles of the denoised images show that NLM, BM3D and the TV methods remove the noise in some regions while fail to do so in some other regions, which indicates that these three methods could not suppress the non-stationary noise very well. This is in accordance with our intuition, since they don’t consider the spatially varying noise levels. The last row of Fig. 3 shows the results of the NWTV model. Compared to the other three methods, the NWTV model could suppress the noise quite well in all the regions of the noisy image, which indicates that NWTV adapts well to the spatially varying noise levels. Another observation is that, for the TV and NWTV methods, the stair-casing artifacts are observable in continuous areas, as shown by the line profiles.

Fig.3

The first column shows the noisy image and the denoised images by NLM, BM3D, ROF and NWTV, respectively. The second and third columns demonstrate the line profiles of the noise-free image and the denoised images. Note that the 180th row covers the upper regions of the images, while the 300th row covers the bottom blocks. Results are obtained with best tuned parameters (in terms of both noise suppression and edge-preserving): μ= 2, λ= 2μfor the ROF model and μ= 20, λ= 10μ, θ= 2.5μ, α= 2for the NWTV model.

The quantitative measures in Table 1 also demonstrate the effectiveness of the proposed NWTV model. By checking the computed quality measures, i.e. PSNR, NRMSE, UQI and SSIM, one can see that our NWTV model outperforms all the competing algorithms.

Table 1

Accuracy comparison between NLM, BM3D, ROF and the NWTV for simple synthesized image

	Noisy	NLM	BM3D	TV	NWTV
PSNR	20.8918	30.4836	35.0074	36.6259	37.6741
NRMSE	0.111710	0.012272	0.004330	0.002983	0.002343
UQI	0.947080	0.993849	0.997826	0.998486	0.998818
SSIM	0.172772	0.733970	0.925919	0.905788	0.934088

Algorithm 1 split-Bregman NWTV Denoising
Initialize: $u^{0} = p, and d_{x}^{0} = d_{y}^{0} = b_{x}^{0} = b_{y}^{0} = 0$
While ∥u^k - u^k - 1 ∥ ₂ > tol do
u^k +1 = G^k
$d_{x}^{k + 1} = shrink (\nabla_{x} u^{k + 1} + b_{x}^{k}, \frac{1}{λ} 1 pt e^{α u^{k + 1}})$
$d_{y}^{k + 1} = shrink (\nabla_{y} u^{k + 1} + b_{y}^{k}, \frac{1}{λ} 1 pt e^{α u^{k + 1}})$
$b_{x}^{k + 1} = b_{x}^{k} + \nabla_{x} u^{k + 1} - d_{x}^{k + 1}$
$b_{y}^{k + 1} = b_{y}^{k} + \nabla_{y} u^{k + 1} - d_{y}^{k + 1}$
k = k + 1
end while

4.2 Performance tests on simulated data

In this subsection, the performance of NWTV will be tested against the KL_PWLS algorithm in the projection domain, and the NWTV_NLM algorithm will be tested against the ANLM and KLPWLS_NLM algorithms in the image domain. When applying the NLM method, the sizes of search- and patch-windows are fixed to 21×21 and 5×5, respectively, just as done in [20]. The sinograms are produced by applying the fan-beam forward projection on two simulated phantoms. The number of detector cells is 700, and the length of each detector cell is set to 0.5 mm. The total number of views is set to 720, which are evenly spanned on a circular track of 360 degrees. The SOD and SDD are set to 550 mm and 800 mm, respectively. The non-stationary Gaussian noise is added to the sinograms, in the same way as described in Fig. 2.

4.2.1 Resolution test

In this experiment, the structure-preserving capability of NWTV will be tested. The ideal image is shown in Fig. 4(a), which consists of four line-pairs and two rings. The size of the ideal image is 512×512, the width of line-pairs in positions “1” and “2” are three pixels and the width of line-pairs in positions “3” and “4” are two pixels. Note that the line-pairs in positions “1” and “3” are approximately along the radial direction of the image, while the line-pairs in positions “2” and “4” are approximately along the angular direction. The non-stationary Gaussian noise with l = 0.002, and α = 1.1 is added to the projection data. Figure 4(b) and (c) show the noise-free and the noisy projection data, respectively. To illustrate the noise better, the line profiles of their 344th rows are plotted in Fig. 4(d). The display windows for both the ideal image and the reconstructed CT image are set to [0, 0.06], while the display windows for the projections are set to [0, 4.5].

Fig.4

The simulated projection data. (a) the ideal image; (b) the projection data of (a); (c) the noisy projection data; (d) line profiles for the 344th rows of the noisy and noise-free projection data.

The results of sinogram restoration by KL_PWLS and the proposed NWTV model are shown in Fig. 5 (first and second rows). At first glance, it seems that both algorithms produce quite similar results. However, when checking the large oscillation parts located at the right part, one could see that NWTV recovers the amplitudes of the oscillation better than KL_PWLS. This is also verified by the quantitative indices computed and illustrated in Table 2 (the rows indicated by “Projection”). All the indices win a bit over KL_PWLS. Another observation is that, NWTV does exhibit observable stair-casing artifacts in this experiment.

Fig.5

The first row shows the noisy and denoised projection data. From left to right: Noisy, NWTV and KL_PWLS respectively. The second row shows the line profiles “344th column” of the projection data. The third row shows the reconstructed CT images, the fourth row and the fifth raw show the zoomed-in images of the red- and blue-framed subregions, respectively. Results are obtained with best tuned parameters (in terms of both noise suppression and edge-preserving): μ= 130, λ= 10μ, θ= 0.2μ, α= 1.1 for the NWTV model and β= 30 for the KL_PWLS model.

Table 2

Accuracy comparison between NWTV, KL_PWLS, ANLM, NWTV_NLM and the KLPWLS_NLM for simulated line-pairs data

		Noisy	NWTV	KL_PWLS
Projection	PSNR	30.7751	38.9912	37.1524
	NRMSE	0.036968	0.005575	0.008513
	UQI	0.981772	0.997196	0.995719
	SSIM	0.288634	0.686740	0.590407
Reconstruction	PSNR	21.6094	25.4410	23.8508
	NRMSE	0.198680	0.082223	0.118581
	UQI	0.897866	0.954492	0.933919
	SSIM	0.400502	0.530847	0.517702
		ANLM	NWTV_NLM	KLPWLS_NLM
Reconstruction	PSNR	26.7891	26.9747	25.7632
	NRMSE	0.060157	0.057760	0.076345
	UQI	0.964783	0.967272	0.956197
	SSIM	0.621822	0.612425	0.607518

In the image domain, the difference of the smoothing behavior of the two algorithms becomes apparent. The third row shows the images reconstructed by applying the FBP algorithm on the noisy as well as the denoised projection data, and two of the line pairs are then zoomed in for better display in the fourth and fifth rows. An interesting observation is that while for the line pair positioned along the angular direction, both the NWTV and KL_PWLS produce satisfying results, as shown in the fifth row. For the line pair positioned along the radial direction, however, the NWTV produces much better result and heavy blurring have been introduced by the KL_PWLS algorithm, as shown in the fourth row. This phenomenon could be explained off by the KL transform. The KL transform is along the angular direction, and the principle component analysis in the KL domain would lead to information losing between adjacent angles, which shall introduce blurring along the angular direction.

The quantitative indices computed and illustrated in Table 2 further verify the advantages of NWTV over the competing algorithm, since all the measures of NWTV in the row indicated by “reconstruction” win over the KL_PWLS algorithm, since all the measures of NWTV in the row indicated by “reconstruction” win over the KL_PWLS algorithm.

Another comparison is performed among the three NLM-based methods: ANLM, NWTV_NLM and KLPWLS_NLM. The results are shown in Fig. 6. A quick observation is that pure post-processing algorithm, i.e. ANLM, could improve the reconstructions to a great deal. Indeed, the line pairs are reconstructed quite well, as shown in the first column. Even so, as shown in the second column, improvement still be observable in the results of NWTV_NLM. With additional NLM denoising (guided by the reconstructed image), both the results of NWTV and KL_PWLS are improved. However, for the KL_PWLS algorithm, due to severe blurring along the angular direction, the line pair framed by the red rectangle still suffers from blurring. As a comparison, just the remaining noise in the results of NWTV reconstruction is removed by the NWTV_NLM algorithm, and all the line pairs are well reconstructed. The quantitative indices computed and illustrated in Table 2 also clearly show the advantages of NWTV_NLM algorithm.

Fig.6

The first row shows the images recovered by the ANLM, NWTV_NLM and KLPWLS_NLM algorithms, respectively. The second row shows the zoomed-in images of the red-framed subregions and the third row shows the zoomed-in images of the blue-framed subregions. Results are obtained with fine tuning parameters of h = 0.6 for ANLM, τ= 0.025 for NWTV_NLM and τ= 0.035 for KLPWLS_NLM.

It needs to be pointed out that ANLM assumes parallel-beams rather than fan-beams. So, to apply ANLM, a transform is needed to rearrange the fan-beams to parallel-beams. This transform, however, introduces interpolation errors, which might lead to blurring for the reconstructed images. This blurring effect can actually be identified in the middle image of the first column of Fig. 6.

4.2.2 Streak artifacts suppression test

In low-dose image reconstructions, streak artifacts are commonly seen due to photon starvation [43, 44]. The photon starvation artifacts are partially caused by the noise in the sinograms, and denoising the sinograms should help to suppress them.

In this experiment, we will try to demonstrate the streaks suppressing capabilities of the proposed NWTV algorithm. The ideal image is shown in Fig. 7(a), which is one slice of the thorax phantom (http://www.imp.uni-erlangen.de/phantoms/thorax/thora-x.htm). The non-stationary Gaussian noise with and α= 1.2 is added to the projection data. Figure 7(b) and (c) show the noise-free and noisy projection data respectively, and Fig. 7(d) shows the line profiles of their 176th rows, which clearly show the signal dependent noise levels. The display windows of the ideal image as well as the reconstructed CT images are [0, 0.09], and the display windows of the corresponding projections are [0, 8.5].

Fig.7

Acquisition of the simulated projection data. (a) the ideal image; (b) the projection data of (a); (c) the noisy projection data; (d) line profiles for the 176th rows of the noisy and noise-free projection data.

The results are shown in Fig. 8. Due to the existence of high absorption objects (i.e. bones) as well as the added noise, streak artifacts have been introduced in the reconstructed image, as shown in the last image of the first column. By applying the two sinogram smoothing methods (i.e. NWTV and KL_PWLS), these artifacts have been effectively removed, without sacrificing the spatial resolution, see Fig. 8 (the last two reconstructions in the third row). It’s worthy to point out that, as shown in the line profiles for the NWTV algorithm, the stair-casing artifacts are observable in this experiment. However, it doesn’t lead to unfavorable artifacts in the reconstructed images. This is understandable if one checks the noisy line profile. It’s hard to say that the denoised line profile for the KL_PWLS is better than the one for NWTV.

Fig.8

The first row shows the noisy and denoised projection data. From left to right: Noisy, NWTV and KL_PWLS respectively. The second row shows the line profiles (176th column) of the projection data. The third row shows the reconstructed CT images. Results are obtained with best tuned parameters (in terms of both streaks suppression and edge-preserving): for the NWTV model and for the KL_PWLS model.

The three NLM-based methods for the reconstructed CT images are also implemented and the results are shown in Fig. 9. Compared to ANLM, NWTV_NLM and KLPWLS_NLM perform better in artifacts suppression as well as structure preserving. Please check the location indicated by the red arrow in Fig. 9(b) and (c). For the result of ANLM, blurring is introduced and the spatial resolution is reduced. Besides, the streak artifacts are still clearly observable. In contrast, NWTV_NLM and KLPWLS_NLM could effectively remove the streaks while leaves the spatial resolution almost the same.

Fig.9

These images are resulted from NLM-based methods. (a) the result of ANLM; (b) the result of NWTV_NLM; (c) the result of KLPWLS_NLM. Results are obtained with best tuned parameters of h = 0.8 for ANLM, τ= 0.025 for NWTV_NLM and τ= 0.035 for KLPWLS_NLM.

4.3 Performance test on real data

The real projection data is obtained through scanning a line-pairs phantom by an industrial CT system. The number of detector cells is 700, and the length of each detector cell is 0.083 mm. The total number of views is set to 720, which are evenly spanned on a circular track of 360 degrees. The SOD and SDD are set to 450 mm and 1246.05 mm, respectively.

A normal dose projection data is obtained as follows. The X-ray source is tuned to have voltage 300 kV and the current intensity is set to 3 mA. Five frames are merged into one frame in each view and the integration time for each frame is set to 0.02 ms. Then the projection data is obtained by applying system calibration and the Log-transform to the acquired raw data. The image reconstructed directly from this projection data, which is shown in Fig. 10 (first image in the third row), will act as the reference image to compare the reconstruction results from the tested algorithms. We decrease the current to 0.75 mA and scan this phantom again to get the low-dose data (about 25% of normal dose). The reconstructed image from this low-dose data is shown in Fig. 10 (first image in the first row). Note that the line pairs are zoomed in and displayed for better illustration. The size of the reconstructions is fixed to 1024×1024 and the display windows are set to [0, 0.18]. In this experiment, the ANLM algorithm is excluded out of the comparison. This is because the transform from fan-beam to parallel-beam is problematic. For the equipment in our lab, the unavoidable nonzero offset and tilt angle of the detector as well as other errors in the geometrical parameters bring difficulties to the rearrangement transform. As a result, heavy blurring (non-uniformly distributed) would be introduced in the ANLM reconstructions.

Fig.10

The first and third rows show the reconstructed images, while the second and fourth rows show the zoomed-in images for the red-framed subregions marked in the first and third rows, respectively. Results are obtained with best tuned parameters (in terms of noise suppression and edge-preserving): μ= 250, λ= 10μ, θ= 0.2μ, α= 1.1 for the NWTV model and β= 10 for the KL_PWLS model, τ= 0.03 for NWTV_NLM and KLPWLS_NLM.

The reconstructed images are shown in Fig. 10. From the zoomed-in images of the red-framed subregions, one could find out that the reconstruction from the NWTV algorithm, as shown in the second image of first row, is really close to the normal dose reconstruction, and the noise has been efficiently reduced compared with the direct FBP reconstruction. In contrast, the reconstruction from the KL_PWLS algorithm, as shown in the last image of first row, is slightly blurred and the amount of remained noise is bigger. When combining with the NLM filtering, i.e. applying the NWTV_NLM and KLPWLS_NLM algorithms, both algorithms produce almost equally good results, as shown in the second and third images in the third row of Fig. 10.

5 Discussion

In this study, we proposed and tested the NWTV regularizer for smoothing the low-dose CT projection data, which incorporates the prior knowledge that the projection data suffers from the non-stationary Gaussian noise. Fast split-Bregman algorithm is formulated and generalized to solve the NWTV model efficiently, and experiments against other popular methods, e.g. KL_PWLS, KLPWLS_NLM and ANLM method, are performed to verify the effectiveness and advantages of the NWTV model. Generally speaking, the performance of NWTV is at least comparable to the competing algorithms. It could demonstrate advantages at balancing the spatially varying noise levels with the smoothing strength, such that the noise is removed in a more uniform manner and blurring artifacts is minimized.

Our experiments suggest that proper post-processing of the reconstructed images enables to further improve their quality. Especially, when employing the idea of the SR_NLM algorithm, the performance difference between KL_PWLS and NWTV is much reduced. Indeed, the difference between the reconstructed images between NWTV_NLM and KLPWLS_NLM is negligible in some of our tests. Another observation is that the ANLM might not work well for very noisy data, which indicates that pre-processing is valuable and indispensable.

TV regularization is usually thought of being pursuing piecewise constant approximations and its drawback of stair-casing artifacts when applied on smooth or piecewise smooth data is well recognized. When dealing with the projection data, however, it seems that the stair-casing artifacts don’t create obvious artifacts in the reconstructed images. The reason behind might be that a small regularization parameter is used such that the smoothing is gentle. Indeed, strong smoothing would introduce streaks to the reconstructed images whether TV or other regularizers are used. On the other hand, it might be interesting to see if NWTV could be improved by replacing TV with other regularizers which prefer smoother functions. For example, in [45], the Huber function is adopted to build a regularization term, which could be thought of approximating TV, while behaves like a quadratic norm (rather than the ℓ₁ norm) for small gradients. Huber regularizer could be thought of pursuing piecewise smooth functions since the quadratic norm prefers smooth approximations. We would like to reserve it for future research.

Regularization is usually done in the CT image domain to bring in the prior information for the reconstructed CT images, since one could have certain knowledge for the scanning samples. This idea has been extensively exploited in the literature, and TV regularization is among the most popular ones. Recent research work [46, 47] suggest that the image quality could be further improved if regularization was done on both image domain and the projection domain. It could be possible to incorporate the NWTV term into some double domain regularization framework to deal with more challenging problems like sparse view or limited-angle reconstruction. This will be investigated in our future work.

Footnotes

Acknowledgments

Thanks for the support of the National Natural Science Foundation of China (NSFC) (61871275 and 61771324). And the authors are grateful to Beijing Higher Institution Engineering Research Center of Testing and Imaging as well as Beijing Advanced Innovation Center for Imaging Technology for funding this research work.

References

D.J.

Brenner and

E. J. J. N. E. J. o. M.

Hall , Computed tomography—an increasing source of radiation exposure 357(22) (2007), 2277–2284.

Sodickson , et al., Recurrent CT, cumulative radiation exposure, and associated radiation-induced cancer risks from CT of adults, Radiology 251(1) (2009), 175–184.

A.C.

Kak and

Slaney , Principles of computerized tomographic imaging, SIAM, 2001.

Zhu ,

Zhao ,

Li and

Zhang , Noise reduction with low dose CT data based on a modified ROF model, Optics Express 20(16) (2012), 17987–18004.

L.I.

Rudin ,

Osher and

Fatemi , Nonlinear total variation based noise removal algorithms, Physica D: Nonlinear Phenomena 60(1-4) (1992), 259–268.

Manduca , et al., Projection space denoising with bilateral filtering and CT noise modeling for dose reduction in CT, Medical Physics 36(11) (2009), 4911–4919.

Zhang ,

Sun and

Zhang , Projection domain denoising method based on dictionary learning for low-dose CT image reconstruction, Journal of X-ray Science and Technology 23(5) (2015), 567–578.

Li , et al., Nonlinear sinogram smoothing for low-dose X-ray CT, IEEE Transactions on Nuclear Science 51(5) (2004), 2505–2513.

Ma , et al., Variance analysis of x-ray CT sinograms in the presence of electronic noise background, Medical Physics 39(7) (2012), 4051–4065.

10.

Wang ,

Li ,

Lu and

Liang , Penalized weighted least-squares approach to sinogram noise reduction and image reconstruction for low-dose X-ray computed tomography, IEEE Transactions on Medical Imaging 25(10) (2006), 1272–1283.

11.

Zhang ,

Zhang and

Lu , Statistical sinogram smoothing for low-dose CT with segmentation-based adaptive filtering, IEEE Transactions on Nuclear Science 57(5) (2010), 2587–2598.

12.

Gordon ,

Bender and

G.T.

Herman , Algebraic reconstruction techniques (ART) for three-dimensional electron microscopy and X-ray photography, Journal of Theoretical Biology 29(3) (1970), 471IN1477–476IN2481.

13.

A.H.

Andersen and

A.C.

Kak , Simultaneous algebraic reconstruction technique (SART): A superior implementation of the ART algorithm, Ultrasonic Imaging 6(1) (1984), 81–94.

14.

E.Y.

Sidky and

Pan , Image reconstruction in circular cone-beam computed tomography by constrained, total-variation minimization, Physics in Medicine and Biology 53(17) (2008), 4777.

15.

Sun ,

Li ,

Cong ,

Wang and

Guo , Sparse-view image reconstruction via total absolute curvature combining total variation for x-ray computed tomography, Journal of X-ray Science and Technology 25(6) (2017), 857–868.

16.

Zhang ,

Ma ,

Wang ,

Liu ,

Lu and

Liang , Statistical image reconstruction for low-dose CT using nonlocal means-based regularization, Computerized Medical Imaging and Graphics 38(6) (2014), 423–435.

17.

M.A.

Al-antari ,

M.A.

Al-masni ,

M.K.

Metwally , et al., Denoising images of dual energy X-ray absorptiometry using non-local means filters, Journal of X-ray Science and Technology 26(3) (2018), 395–412.

18.

Xu ,

Yu ,

Mou ,

Zhang ,

Hsieh and

Wang , Low-dose X-ray CT reconstruction via dictionary learning, IEEE Transactions on Medical Imaging 31(9) (2012), 1682–1697.

19.

Tan , et al., Tensor-based dictionary learning for dynamic tomographic reconstruction, Physics in Medicine and Biology 60(7) (2015), 2803.

20.

Bian , et al., SR-NLM: A sinogram restoration induced non-local means image filtering for low-dose computed tomography, Computerized Medical Imaging and Graphics 37(4) (2013), 293–303.

21.

Li , et al., Adaptive nonlocal means filtering based on local noise level for CT denoising, Medical Physics 41(1) 2014.

22.

P.F.

Feruglio ,

Vinegoni ,

Gros ,

Sbarbati and

Weissleder , Block matching 3D random noise filtering for absorption optical projection tomography, Physics in Medicine and Biology 55(18) (2010), 5401.

23.

Kang , et al. , Image denoising of low-radiation dose coronary CT angiography by an adaptive block-matching 3D algorithm 8669 (2013), 86692G.

24.

Chen , et al., Improving abdomen tumor low-dose CT images using a fast dictionary learning based processing, Physics in Medicine and Biology 58(16) (2013), 5803.

25.

Chen , et al., Artifact suppressed dictionary learning for low-dose CT image processing, IEEE Transactions on Medical Imaging 33(12) (2014), 2271–2292.

26.

Chen , et al., Low-dose CT denoising with convolutional neural network (2017), 143–146.

27.

J.M.

Wolterink ,

Leiner ,

M.A.

Viergever and

I. J. I. t. o. m. i.

Išgum , Generative adversarial networks for noise reduction in low-dose CT, IEEE Transactions on Medical Imaging 36(12) (2017), 2536–2545.

28.

Natterer and

Wubbeling , Mathematical Methods in Image Reconstruction, SIAM, 2001.

29.

Wu , On the continuity of images by transmission imaging, Communications in Mathematical Sciences 11(2) (2013), 573–595.

30.

V.N.

Varghees ,

M.S.

Manikandan and

Gini , Adaptive MRI image denoising using total-variation and local noise estimation, 2012, pp. 506–511.

31.

Y.L.

You ,

Xu ,

Tannenbaum and

Kaveh , Behavioral analysis of anisotropic diffusion in image processing, IEEE Transactions on Image Processing 5(11) (1996), 1539–1553.

32.

Hsieh , Computed tomography: Principles, design, artifacts, and recent advances. SPIE press, 2003.

33.

Macovski , Medical Imaging Systems. Prentice Hall, 1983.

34.

Salmon ,

Harmany ,

C.-A.

Deledalle and

Willett , Poisson noise reduction with non-local PCA, Journal of Mathematical Imaging and Vision 48(2) (2014), 279–294.

35.

Wang ,

Zhu and

Li , Imaging model for the scintillator and its application to digital radiography image enhancement, Optics Express 23(26) (2015), 33753–33776.

36.

Goldstein and

Osher , The split Bregman method for L1-regularized problems, SIAM Journal on Imaging Sciences 2(2) (2009), 323–343.

37.

Wu and

X.-C.

Tai , Augmented Lagrangian method, dual methods, and split Bregman iteration for ROF, vectorial TV, and high order models, SIAM Journal on Imaging Sciences 3(3) (2010), 300–339.

38.

Chen ,

Lan and

Ouyang , Optimal primal-dual methods for a class of saddle point problems, SIAM Journal on Optimization 24(4) (2014), 1779–1814.

39.

Ouyang ,

Chen ,

Lan and

Pasiliao Jr, An accelerated linearized alternating direction method of multipliers, SIAM Journal on Imaging Sciences 8(1) (2015), 644–681.

40.

Huynh-Thu and

Ghanbari , Scope of validity of PSNR in image/video quality assessment, Electronics Letters 44(13) (2008), 800–801.

41.

Wang ,

A.C.

Bovik ,

H.R.

Sheikh and

E.P.

Simoncelli , Image quality assessment: From error visibility to structural similarity, IEEE Transactions on Image Processing 13(4) (2004), 600–612.

42.

Wang and

A.C.

Bovik , A universal image quality index, IEEE Signal Processing Letters 9(3) (2002), 81–84.

43.

J. F.

Barrett and

Keat , Artifacts in CT: Recognition and avoidance, Radiographics 24(6) (2004), 1679–1691.

44.

Mori ,

Machida ,

Osanai and

Iinuma , Photon starvation artifacts of X-ray CT: Their true cause and a solution, Radiological Physics and Technology 6(1) (2013), 130–141.

45.

K. J.

Little and L. R. P. J., Sinogram restoration in computed tomography with an edge-preserving penalty, Medical Physics 42(3) (2015), 1307–1320.

46.

Zhan and

Dong , CT Image Reconstruction by Spatial-Radon Domain Data-Driven Tight Frame Regularization, SIAM Journal on Imaging Sciences 9(3) (2016), 1063–1083.

47.

Zhang ,

Dong and

Liu , A Reweighted Joint Spatial-Radon Domain CT Image Reconstruction Model for Metal Artifact Reduction, SIAM Journal on Imaging Sciences 11(1) (2018), 707–733.