# Adaptively Tuned Iterative Low Dose CT Image Denoising

Authors: SayedMasoud Hashemi; Narinder S. Paul; Soosan Beheshti; Richard S. C. Cobbold
Journal: Computational and Mathematical Methods in Medicine
 (2015)

---

## Abstract

Improving image quality is a critical objective in low dose computed tomography (CT) imaging and is the primary focus of CT image denoising. State-of-the-art CT denoising algorithms are mainly based on iterative minimization of an objective function, in which the performance is controlled by regularization parameters. To achieve the best results, these should be chosen carefully. However, the parameter selection is typically performed in an ad hoc manner, which can cause the algorithms to converge slowly or become trapped in a local minimum. To overcome these issues a noise confidence region evaluation (NCRE) method is used, which evaluates the denoising residuals iteratively and compares their statistics with those produced by additive noise. It then updates the parameters at the end of each iteration to achieve a better match to the noise statistics. By combining NCRE with the fundamentals of block matching and 3D filtering (BM3D) approach, a new iterative CT image denoising method is proposed. It is shown that this new denoising method improves the BM3D performance in terms of both the mean square error and a structural similarity index. Moreover, simulations and patient results show that this method preserves the clinically important details of low dose CT images together with a substantial noise reduction.


---

## Body

## 1. Introduction

While X-ray computed tomography (CT) enables ultrafast acquisition of patient images obtained with excellent spatial resolution, the dose needed to achieve diagnostic image quality can result in a significant increase in the risk of developing cancer [1]. Consequently, low-dose CT imaging is clinically desired and has been under investigation for several years. Lowering the radiation dose may seriously degrade diagnostic performance or undermine physician confidence by producing noisier images [2, 3]. Several different algorithmic approaches have been proposed to reduce the effect of noise in the CT images, including projection data denoising [4–6], optimizing the reconstruction algorithms to include the noise statistics [7–9], and CT image denoising [10–12]. The latter is the focus of this paper, where an adaptively tuned iterative CT image denoising algorithm is presented.The main source of noise in X-ray projection data is quantum noise caused by statistical fluctuations of X-ray quanta reaching the detectors, so that the CT projection noise follows the Poisson distribution [3]. However, because of the use of different reconstruction algorithms and signal processing steps in CT reconstruction, the noise statistics of the processed CT images are usually unknown and hard to model and are spatially changing. Moreover, directional noise in the form of streak artifacts is present in many CT images. As a result, incorporating accurate noise statistics into image-based CT denoising can be very challenging. When denoising is based on the projection data and its statistics, other difficulties arise. Specifically, such denoising methods and the associated iterative reconstructions require access to the CT raw data, which is often unavailable. Furthermore, these methods have a high computational complexity making it challenging to obtain a final image in a reasonable length of time, depending on the available computational resources. On the other hand, image-based denoising methods are fast and can be applied directly on the CT images without changing the clinical workflow.A simplified noise model is usually used in image based denoising algorithms, in which, following the Central Limit Theorem (CLT) [13], the final noise in each voxel follows a Gaussian distribution [11, 14–17]. The CLT can be used since each voxel in CT images is computed by adding values from many different projections. With this assumption, a noisy CT image y can be modeled by(1)y=x+n,where x is the noiseless image and n is a zero mean additive anisotropic Gaussian noise with variance of σn2, which varies with the pixel location and its value.Different image based denoising algorithms have been used to estimate the noiseless CT images, such as anisotropic diffusion [10], total variation (TV) [18], bilateral filtering [19], or wavelet-based techniques [11, 12, 20]. These methods can usually be formulated as an unconstrained Lagrangian multiplier optimization problem [18, 21–23], that is,(2)x^=argminx⁡12x-y22+λhx,in which λ is a regularization parameter that controls the tradeoff between the data fidelity and the regularization term h(x) and xq=∑i=1n‍xiq. Different regularization terms, h(x), lead to different denoising methods. For example, in TV-based methods [24–26], h(x)=δhx2+δvx2 where δh and δv are the gradient in horizontal and vertical directions, and in wavelet soft thresholding methods h(x)=∑i‍|Ψ2Dx)i, where Ψ2D(x) is the 2D wavelet transform [27].There is a strong dependence of the quality of the result on the regularization parameter. It is a challenging task to find the regularization parameterλ that provides the best balance between signal smoothing and feature preservation [26]. Specifically, if λ is not appropriately adjusted, the optimization is trapped in a local minimum; that is, if λ is too small, noise is only partially removed and, if it is too large, the image may be oversmoothed [25]. Some methods have been proposed to update the regularization parameters iteratively, such as use of the discrepancy principle [28], generalized cross-validation [26], and L-curve [29]. These methods fail in certain situations, are problem specific, and generally increase the computational complexity of the algorithms.One straight forward approach used in many algorithms is to use a heuristicλ value combined with a criterion to stop the algorithm before the estimated signal is oversmoothed. Different stopping criteria have been proposed for iterative denoising problems. For instance, Akkoul et al. [30] used a switching median filter algorithm to stop the iterative process when the number of changed pixels in the denoising iterations is a minimum. In [20, 31] the statistical properties of high frequency wavelet subbands were used to stop the TV iterations. However, such methods are unable to differentiate oversmoothed data from well-denoised data. As a result, to avoid oversmoothing, the updating steps are typically chosen to be small, which decreases the convergence speed.In this paper, a noise confidence region evaluation (NCRE) method is used to address the regularization selection and the algorithm stopping problems. It adaptively updates the regularization parameters at the end of each iteration by validating the result of that iteration. The algorithm stops when the statistical properties of the denoising residual resemble those of the additive white Gaussian noise. Using NCRE, a new iterative block matching and 3D filtering (BM3D) method is proposed, which outperforms BM3D [32] itself. The proposed method is compared with anisotropic diffusion denoising, which is generally regarded as a standard denoising method in CT imaging [10], nonlocal mean [33–35], and BM3D [32]. In addition, we study the noise properties of CT images and show that the noise in small image blocks has an additive white Gaussian model, which justifies the success of the nonlocal based denoising algorithms in image based CT denoising [33, 36–39].


## 2. Problem Formulation and CT Noise Properties

Recently, it has been shown that nonlocal patch based algorithms outperform others in CT image denoising [33, 34, 36–41]. For example, in [41] a nonlocal means (NLM) based method, which takes advantage of the presence of repeating structures in a given image, was compared with a principle component analysis based denoising method and a highly constrained backprojection method. It was shown that the NLM method outperformed both of the methods in terms of the contrast to noise ratio, noise standard deviation, and squared error. Another class of algorithms looks for similar blocks in the whole 2D image and stacks them together in 3D arrays. Denoising is then performed through transform domain shrinkage of the 3D arrays. An algorithm called K-SVD [37, 38, 42] uses these 3D patches to train an optimum dictionary. This method, which assumes that each 3D block can sparsely be represented by the trained dictionary atoms, uses shrinkage algorithms to denoise the patches.Our proposed algorithm makes use of the block matching and 3D filtering (BM3D) technique [32, 36]. This is a noniterative denoising method that currently outperforms many newer algorithms [40]. It is composed of two major filtering steps. In both stages collaborative filtering is utilized, which itself has four stages: (1) grouping similar patches with a reference patch, (2) calculation of the 3D wavelet coefficients of each stack of patches, (3) denoising the wavelet coefficients (thresholding in step (1) or Wiener filtering in step (2)), and (4) recovering the denoised image by calculating the inverse 3D wavelet transformation. The BM3D approach aims to denoise the patches by Wiener filtering, which is done in step (2). To find best match to similar patches in the step (2) and to reliably estimate the Wiener coefficients, the method requires a reliable estimate of the noiseless image, which is the main purpose of step (1). The input for this step, which is a hard thresholding block, is the 3D noisy wavelet coefficients of similar patches located by block matching applied to the available noisy image. A hard threshold with a heuristically determined value of 2.7σn is used in step (1). The resulting denoised coefficients are then transformed back to a spatial domain to be used as the initial estimate of the noiseless data used for calculating the Wiener filter coefficients.To denoise the images, patch based methods generally use a model based on(3)yGj=xGj+nGj,j∈1,…,m,in which G is the patch grouping information, m is the number of 3D patches, yGj denotes the noisy patches, xGj denotes the noiseless 3D patches, and nGj~N(0,σGi2) is the noise at each 3D patch. Conventional BM3D uses an independent identical additive Gaussian noise model in which the noise variances σGi2 are similar in all patches. Using this assumption, its regularization term is a nonlocal wavelet l1-norm h(x)=NLW(x,G), where [43] (4)NLWx,G=∑j=1mΨ3DxGj1,in which Ψ3D is the 3D wavelet transform. Using this regularization in (2), BM3D solves the following optimization problem in its first step:(5)x^=argminx⁡12x-y22+λ∑j=1mΨ3DxGj1.In our approach, we modify the BM3D formulation for CT image denoising by incorporating a more realistic noise model, to include the nonstationarity of the noise and its dependence on the position and value of the pixels. Our proposed method uses the noise properties of the patches xGj, studied in the Appendix, to improve the performance of BM3D for CT image denoising.


### 2.1. Noise in CT Images

Although a reasonable statistical model for the CT projection data is the independent Poisson distributions [3], it has been shown that the corrected polyenergetic X-ray projections can be modeled more accurately by a Gaussian distribution with the following relationship between its mean and variance:(6)σi,n2=Ji×exp⁡ξ¯i,ns2+σe2i,j,where ξ¯ is the mean and σi,n2 is the variance of the projections at the ith projection angle (φi) and the nth detector bin whose distance is ln from the detectors center, s is a scaling factor, σe2 is the electronic noise variance, and Ji is a parameter adaptive to different detector channels [6]. During the reconstruction process the noise distribution is changed by the reconstruction algorithm and filters. As a result, due to the complicated dependencies of noise on scan parameters and on spatial position, the noise distribution in the reconstructed CT images is usually very difficult to determine (a detailed study of the noise model in CT can be found in [44, 45]). Using the discrete filtered back projection relation,(7)fx,y=πΔtNπg∑i=1Nπg∑n=1Ngcxcos⁡φi-ysin⁡φi-lngi,n,the noise variance in the reconstructed images can be described by [38](8)σn2x,y=πΔtNπg2∑i=1Nπg∑n=1Ngc2xcos⁡φi-ysin⁡φi-lnσi,n2,where Δt is the distance between the center of two adjacent detectors, gi,n is the parallel projection at the ith angle and the nth detector bin, and c(·) is the ramp filter in the spatial domain. This can be interpreted as the backprojection of the projection noise variances, making it nonstationary, object dependent, and correlated. Moreover, since the variance of each voxel is the summation of the variances from many angles, if the variance in one direction is significantly larger than that at another direction, then the variances along that direction will be more correlated than that for other directions [46] producing what is known as a streak artifact in the reconstructed images. It should be noted that this effect is not included in our model. In (8), the ramp filter c(·) is symmetric about the center of rotation and σi,n2 depends on the attenuation of the media through which the X-ray beams pass. Therefore, we assert that the noise variances of small neighborhoods with similar attenuations and similar radial distances from the center of rotation should be very similar. The results of experimental tests of this assertion are presented in the Appendix.


### 2.2. Modified Formulation

Based on the above discussion and our experimental results, presented in the Appendix, it can be assumed that the noise in 3D similar patches of the CT images follows a white additive Gaussian distribution, with different variances for different 3D patchesσGi2. Using (3), the modified optimization problem used in this paper is given by(9)Jy,G,Λ=x^=argminx⁡12x-y22+∑j=1mλjΨ3DxGj1,in which Λ=[λ1,λ2,…,λm] is a set of regularization parameters that are functions of σGi. To improve the regularization parameter selection, we used an adaptive updating method based on an evaluation of the noise statistics. It automatically avoids the oversmoothing without lowering the convergence speed. The NCRE method validates the statistical properties of the error residuals at the end of each iteration and categorizes the result as well denoised, partially denoised, or oversmoothed. This information is then used to update the parameters in the next iteration or to stop the algorithm at the end of the iteration for which the similarity between error residual and Gaussian noise is satisfied.


## 2.1. Noise in CT Images

Although a reasonable statistical model for the CT projection data is the independent Poisson distributions [3], it has been shown that the corrected polyenergetic X-ray projections can be modeled more accurately by a Gaussian distribution with the following relationship between its mean and variance:(6)σi,n2=Ji×exp⁡ξ¯i,ns2+σe2i,j,where ξ¯ is the mean and σi,n2 is the variance of the projections at the ith projection angle (φi) and the nth detector bin whose distance is ln from the detectors center, s is a scaling factor, σe2 is the electronic noise variance, and Ji is a parameter adaptive to different detector channels [6]. During the reconstruction process the noise distribution is changed by the reconstruction algorithm and filters. As a result, due to the complicated dependencies of noise on scan parameters and on spatial position, the noise distribution in the reconstructed CT images is usually very difficult to determine (a detailed study of the noise model in CT can be found in [44, 45]). Using the discrete filtered back projection relation,(7)fx,y=πΔtNπg∑i=1Nπg∑n=1Ngcxcos⁡φi-ysin⁡φi-lngi,n,the noise variance in the reconstructed images can be described by [38](8)σn2x,y=πΔtNπg2∑i=1Nπg∑n=1Ngc2xcos⁡φi-ysin⁡φi-lnσi,n2,where Δt is the distance between the center of two adjacent detectors, gi,n is the parallel projection at the ith angle and the nth detector bin, and c(·) is the ramp filter in the spatial domain. This can be interpreted as the backprojection of the projection noise variances, making it nonstationary, object dependent, and correlated. Moreover, since the variance of each voxel is the summation of the variances from many angles, if the variance in one direction is significantly larger than that at another direction, then the variances along that direction will be more correlated than that for other directions [46] producing what is known as a streak artifact in the reconstructed images. It should be noted that this effect is not included in our model. In (8), the ramp filter c(·) is symmetric about the center of rotation and σi,n2 depends on the attenuation of the media through which the X-ray beams pass. Therefore, we assert that the noise variances of small neighborhoods with similar attenuations and similar radial distances from the center of rotation should be very similar. The results of experimental tests of this assertion are presented in the Appendix.


## 2.2. Modified Formulation

Based on the above discussion and our experimental results, presented in the Appendix, it can be assumed that the noise in 3D similar patches of the CT images follows a white additive Gaussian distribution, with different variances for different 3D patchesσGi2. Using (3), the modified optimization problem used in this paper is given by(9)Jy,G,Λ=x^=argminx⁡12x-y22+∑j=1mλjΨ3DxGj1,in which Λ=[λ1,λ2,…,λm] is a set of regularization parameters that are functions of σGi. To improve the regularization parameter selection, we used an adaptive updating method based on an evaluation of the noise statistics. It automatically avoids the oversmoothing without lowering the convergence speed. The NCRE method validates the statistical properties of the error residuals at the end of each iteration and categorizes the result as well denoised, partially denoised, or oversmoothed. This information is then used to update the parameters in the next iteration or to stop the algorithm at the end of the iteration for which the similarity between error residual and Gaussian noise is satisfied.


## 3. Proposed Denoising Method: BM3D-NCRE

Ifx^i denotes the signal recovered at ith iteration, the denoising residual at the end of ith iteration can be expressed by Δyi=y-x^i which, ideally, should be the noise (3), n={nGi,∀i=1,…,m}. Here, we provide a quantitative measure that verifies the similarity between the structure of Δyi and that of n.


### 3.1. Noise Confidence Region Evaluation (NCRE)

In [47] it was shown that the following function of zero mean white Gaussian noise, n, with length n, and for any given scalar value of z:(10)Gz,n=1n∑j=1nGz,nj,Gz,nj=1ifnj≤z0ifnj>z;is equivalent to sorting the absolute value of the noise elements nj. The expected value of this function is E(G(z,n))=F(z) and its variance is var(G(z,n))=F(z)(1-F(z))/n, where F(z)=2ϕ(z/σn)-1 and ϕ(·) is the cumulative distribution function (CDF) of a Gaussian distribution. Therefore, G(z,n) is bounded by the following lower (Ln) and upper (Un) values: (11)Lnz=Fz-ζ1nFz1-Fz,Unz=Fz+ζ1nFz1-Fzwith probability of Pr{|G(z,n)-F(z)|≤ζ1/mF(z)(1-F(z))}≈2ϕ(ζ)-1. If the sorted absolute values of a signal lie between these two boundaries for a large enough ζ, that signal will follow a white Gaussian distribution with a confidence probability close to one.As shown in Figure1, these boundaries divide the z,G(z,n) space into three regions. At the end of each iteration G(z,Δyi) the sorted absolute value of the residual is calculated. If this sequence falls into Region II (in our proposed algorithm being a subset of a region is evaluated by having a high fraction of G(z,Δyi) in that region; for example, in our simulations this fraction is 90%), it means that the residual has a Gaussian-like structure and denoising stops. On the other hand, if the denoising at the ith iteration has removed not only the noise but also parts of the noiseless data itself, Δyi will have some of the image information making its samples larger than Gaussian noise. Therefore, for a specific value of z, G(z,Δyi) (average number of Δyis with absolute values smaller than z) is smaller than G(z,n) and falls in Region III. This will enforce continuation of the denoising to the (i+1)th step with changing of the regularization parameters λj such that the denoising algorithm extracts less noise in the next iteration, that is,decreasing∀λj=λj/s, s>1. If G(z,Δyi) falls in Region I, it indicates that the noise is partially removed. In this case the algorithm continues to the (i+1)th step and changes the regularization parameters λj such that more noise is extracted by the denoising algorithm, that is,increasing∀λj=s×λj, s>1. In summary, at each iteration when G(z,Δyi) is in either Region I or III, the regularization parameter is updated such that it moves toward Region II. The value of s can be tuned as a fixed or an adaptively changing variable based on the euclidean distance between G(z,Δyi) and F(z), di=Gz,Δyi-Fz2. In our proposed method to update Λ=[λ1,…,λm] a global s is used (similar to [33, 41]), which is updated based on the placement of the denoising residual of the patches in different Regions I–III. The algorithm is stopped when the denoising residual of the soft tissue around the lung is placed in Region II. An example of a soft tissue mask for a thoracic phantom, denoted by M, is shown in Figure 2. The pixels in this region have very similar CT# and have almost the same radial distances from the center of the rotation. Therefore, it could be assumed that the noise in this region has a white Gaussian distribution.Figure 1
Three possible regions for the residual at the end of each iteration. If it lies in Region II (noise confidence region), denoising is stopped.Figure 2
Noise statistics of the soft tissue region surrounding the lung. (a) Top: the thoracic phantom which is used to evaluate the noise characteristics and bottom: the soft tissue region of the phantom (denoted byM in Algorithm 1). (b) The statistical distribution of the noise in the soft tissue region is shown by the blue experimental values and these are compared with Gaussian distribution with the same variance (red dashed line).
(a)
(b)Algorithm 1:Proposed Iterative Regularization Parameter Updating.
Initialize: i=1, M: Soft tissue mask,
x
^
0
=
J
(
y
,
G
i
,
Λ
i
), Δy0=y-x^0, G0=sort(M·Δy0),
while  M·Gi  ⊈  Region II and  i < maxiter  do,
if  Gi⊂RegionI  then
increase  ∀λj∈Λ,j=[1,…,m]  towards Region II,
x^i=J(yi,Gi,Λi),
yi=(1-α1)y+β1x^i+(α1-β1)x^i-1,
end if
if  Gi ⊂ Region III then
decrease  ∀λj∈Λi,j=[1,…,m]  towards Region II,
x^i=J(yi,Gi,Λi),
yi=(1-α2)y+β2x^i+(α2-β2)x^i-1,
end if
Δyi=y-x^i, Gi=sort(M·Δyi), i←i+1,
end while


### 3.2. Summary of the Proposed Method: BM3D-NCRE

Algorithm1 shows the proposed iterative Λ updating scheme, in which · denotes an element-wise multiplication. The updating method uses a memory strategy for recovery of possible lost edges and fine details. This process is represented by(12)yi=1-αky+βkx^i+αk-βkx^i-1,in which k∈[1,2] and αk and βk are positive scalars, chosen based on the conditions given in [48]. This stage of the algorithm was inspired by the second order iterative methods [48, 49] that improve the convergence rate of the iterative methods.To denoise CT images the fundamentals of BM3D were used in an iterative scheme: the output of the Wiener filter is a better estimate of the original image than the input of Wiener filter from the first step. Therefore, this output can be fed into the first step to provide better Wiener coefficients in the second iteration. The modified BM3D formulation (9) is used iteratively in Algorithm 1 to denoise the CT images, where NCRE adjusts the threshold values applied on 3D wavelet coefficients of the first step. The parameters of the BM3D algorithm are chosen based on the ones that resulted in the best performance in [50]. The initial value for Λ is 2.7σGi with noise variances estimated independently for each 3D stack, using the median-absolute deviation method as described in [51]. In each iteration, if the sorted absolute value of the denoising residual of the soft tissues around the lung, M·Δyi, falls into Region I, the threshold value will be increased and, if it falls into Region III, the threshold values are decreased, so that in the next iteration the residual moves towards Region II.


## 3.1. Noise Confidence Region Evaluation (NCRE)

In [47] it was shown that the following function of zero mean white Gaussian noise, n, with length n, and for any given scalar value of z:(10)Gz,n=1n∑j=1nGz,nj,Gz,nj=1ifnj≤z0ifnj>z;is equivalent to sorting the absolute value of the noise elements nj. The expected value of this function is E(G(z,n))=F(z) and its variance is var(G(z,n))=F(z)(1-F(z))/n, where F(z)=2ϕ(z/σn)-1 and ϕ(·) is the cumulative distribution function (CDF) of a Gaussian distribution. Therefore, G(z,n) is bounded by the following lower (Ln) and upper (Un) values: (11)Lnz=Fz-ζ1nFz1-Fz,Unz=Fz+ζ1nFz1-Fzwith probability of Pr{|G(z,n)-F(z)|≤ζ1/mF(z)(1-F(z))}≈2ϕ(ζ)-1. If the sorted absolute values of a signal lie between these two boundaries for a large enough ζ, that signal will follow a white Gaussian distribution with a confidence probability close to one.As shown in Figure1, these boundaries divide the z,G(z,n) space into three regions. At the end of each iteration G(z,Δyi) the sorted absolute value of the residual is calculated. If this sequence falls into Region II (in our proposed algorithm being a subset of a region is evaluated by having a high fraction of G(z,Δyi) in that region; for example, in our simulations this fraction is 90%), it means that the residual has a Gaussian-like structure and denoising stops. On the other hand, if the denoising at the ith iteration has removed not only the noise but also parts of the noiseless data itself, Δyi will have some of the image information making its samples larger than Gaussian noise. Therefore, for a specific value of z, G(z,Δyi) (average number of Δyis with absolute values smaller than z) is smaller than G(z,n) and falls in Region III. This will enforce continuation of the denoising to the (i+1)th step with changing of the regularization parameters λj such that the denoising algorithm extracts less noise in the next iteration, that is,decreasing∀λj=λj/s, s>1. If G(z,Δyi) falls in Region I, it indicates that the noise is partially removed. In this case the algorithm continues to the (i+1)th step and changes the regularization parameters λj such that more noise is extracted by the denoising algorithm, that is,increasing∀λj=s×λj, s>1. In summary, at each iteration when G(z,Δyi) is in either Region I or III, the regularization parameter is updated such that it moves toward Region II. The value of s can be tuned as a fixed or an adaptively changing variable based on the euclidean distance between G(z,Δyi) and F(z), di=Gz,Δyi-Fz2. In our proposed method to update Λ=[λ1,…,λm] a global s is used (similar to [33, 41]), which is updated based on the placement of the denoising residual of the patches in different Regions I–III. The algorithm is stopped when the denoising residual of the soft tissue around the lung is placed in Region II. An example of a soft tissue mask for a thoracic phantom, denoted by M, is shown in Figure 2. The pixels in this region have very similar CT# and have almost the same radial distances from the center of the rotation. Therefore, it could be assumed that the noise in this region has a white Gaussian distribution.Figure 1
Three possible regions for the residual at the end of each iteration. If it lies in Region II (noise confidence region), denoising is stopped.Figure 2
Noise statistics of the soft tissue region surrounding the lung. (a) Top: the thoracic phantom which is used to evaluate the noise characteristics and bottom: the soft tissue region of the phantom (denoted byM in Algorithm 1). (b) The statistical distribution of the noise in the soft tissue region is shown by the blue experimental values and these are compared with Gaussian distribution with the same variance (red dashed line).
(a)
(b)Algorithm 1:Proposed Iterative Regularization Parameter Updating.
Initialize: i=1, M: Soft tissue mask,
x
^
0
=
J
(
y
,
G
i
,
Λ
i
), Δy0=y-x^0, G0=sort(M·Δy0),
while  M·Gi  ⊈  Region II and  i < maxiter  do,
if  Gi⊂RegionI  then
increase  ∀λj∈Λ,j=[1,…,m]  towards Region II,
x^i=J(yi,Gi,Λi),
yi=(1-α1)y+β1x^i+(α1-β1)x^i-1,
end if
if  Gi ⊂ Region III then
decrease  ∀λj∈Λi,j=[1,…,m]  towards Region II,
x^i=J(yi,Gi,Λi),
yi=(1-α2)y+β2x^i+(α2-β2)x^i-1,
end if
Δyi=y-x^i, Gi=sort(M·Δyi), i←i+1,
end while


## 3.2. Summary of the Proposed Method: BM3D-NCRE

Algorithm1 shows the proposed iterative Λ updating scheme, in which · denotes an element-wise multiplication. The updating method uses a memory strategy for recovery of possible lost edges and fine details. This process is represented by(12)yi=1-αky+βkx^i+αk-βkx^i-1,in which k∈[1,2] and αk and βk are positive scalars, chosen based on the conditions given in [48]. This stage of the algorithm was inspired by the second order iterative methods [48, 49] that improve the convergence rate of the iterative methods.To denoise CT images the fundamentals of BM3D were used in an iterative scheme: the output of the Wiener filter is a better estimate of the original image than the input of Wiener filter from the first step. Therefore, this output can be fed into the first step to provide better Wiener coefficients in the second iteration. The modified BM3D formulation (9) is used iteratively in Algorithm 1 to denoise the CT images, where NCRE adjusts the threshold values applied on 3D wavelet coefficients of the first step. The parameters of the BM3D algorithm are chosen based on the ones that resulted in the best performance in [50]. The initial value for Λ is 2.7σGi with noise variances estimated independently for each 3D stack, using the median-absolute deviation method as described in [51]. In each iteration, if the sorted absolute value of the denoising residual of the soft tissues around the lung, M·Δyi, falls into Region I, the threshold value will be increased and, if it falls into Region III, the threshold values are decreased, so that in the next iteration the residual moves towards Region II.


## 4. Results and Discussion

Three test methods were used to evaluate the performance of the proposed algorithm. The first method consisted of a simulated Shepp-Logan phantom corrupted by adding Poisson noise to the fan beam X-ray projections. The images were reconstructed using theifanbeam command in Matlab. The number of unattenuated photons in each projection was taken to be 1×1013, which led to a noise variance similar to the images reconstructed from the tube current of 50 mAs and peak voltage of 120 kVp. White Gaussian noise with standard deviation of 100 was added to all the projections to simulate the presence of electronic noise. The second method used a CATPHAN phantom (Phantom Laboratory, Greenwich, NY, USA). This is a standard phantom widely used for CT image quality evaluation and contains spheres of differing contrast as well as line pairs with differing spacing that can be used to test the spatial resolution. Ideally, the denoising algorithms should enable us to distinguish between smaller spheres with lower contrasts in the low contrast slice and should keep the line pair resolution unchanged. The third method uses axial chest CT images from a clinical patient. All three test methods used the following parameters in the Algorithm 1: ζ=6, α1=0.9, β1=0.8, α2=0.4, and β2=0.3. The parameters of BM3D are chosen similar to the ones proposed in [50]: the size of the blocks is 8×8, sliding step to process every next reference block is 3, maximum number of similar blocks is 16, and the size of the search neighborhood for full-search blockmatching is 39×39.Figure3 shows the mean square error (MSE) and the structural similarity index (SSIM) [52] resulting from successive iterations of BM3D-NCRE applied to the reconstructed noisy Shepp-Logan phantom images. The first iteration is equivalent to the result of BM3D. As shown, the MSE decreases and the SSIM increases in successive iterations to a point where the algorithm is stopped by falling into Region II. The results of denoising the Shepp-Logan phantom with BM3D and BM3D-NCRE are shown in Figure 4. As can be seen, the noise is removed more effectively by BM3D-NCRE. However, the streak artifacts are still visible in both denoised images.Figure 3
Squared error (blue line) and the Structural Similarity Index (SSIM) (green dashed line) showing the changes with each iteration when using BM3D-NCRE. The shading colors show the region in whichΔyi’s are placed after each iteration. The algorithm stops when Region II is reached.Figure 4
Shepp-Logan phantom simulations: (a) original phantom, (b) noisy reconstructed phantom, (c) denoised by BM3D, and (d) denoised by BM3D-NCRE.
(a)
(b)
(c)
(d)In the second test, the CATPHAN phantom was scanned using a low dose (50 mAs, 120 kVp) and a high dose (300 mAs, 120 kVp) protocol. Image reconstructions were performed with a Toshiba Aquilion One CT using the proprietary lung kernel FC52 and the proprietary iterative reconstruction algorithm Adaptive Iterative Dose Reduction 3D (AIDR3D) [53]. The latter uses anisotropic diffusion denoising as its base to improve the image quality at each iteration. Our proposed denoising method was applied to the images reconstructed with the high spatial resolution filter algorithm, FC52. These are compared to the images reconstructed with AIDR3D, and the FC52 reconstructed images denoised by nonlocal mean and BM3D. The nonlocal mean package provided by Gabriel Peyre on Mathwork File Exchange was used here (http://www.mathworks.com/matlabcentral/fileexchange/13619-toolbox-non-local-means). This package is based on the method described in [35]. The BM3D code is also based on the package provided by Alessandro Foi on his homepage (http://www.cs.tut.fi/~foi/GCF-BM3D/). It should be noted that the parameters of nonlocal mean and BM3D were heuristically adjusted to achieve the best performance, based on visual inspection of the results. In addition, the parameters are adjusted to keep the spatial resolution in the line pair resolution slice the same. Figure 5 shows the line pair slice reconstructed by FC52, AIDR3D, and FC52 denoised by BM3D-NCRE, nonlocal mean, and BM3D. As can be seen, all these methods have the same spatial resolution as the original image whose resolution was not improved by the proprietary iterative reconstruction method used (AIDR3D).Figure 5
Top: line pair slice of the CATPHAN phantom scanned with 50 mAs and 120 kVp (window width/window level = 400/60 HU). Bottom: red rectangular ROI of the images with the red line showing the cut-off line pair resolution (window width/window level = 400/500 HU). Left to right: image reconstructed with FC52 (STD = 64 HU) and reconstructed with AIDR3D (STD = 41 HU), FC52 denoised with the proposed method (STD = 22 HU), FC52 denoised with nonlocal mean (STD = 34 HU), and FC52 denoised with BM3D method (STD = 27 HU).Figure6 shows that the detectability of low contrast objects is improved with our method and outperforms that achieved by AIDR3D, nonlocal mean, and BM3D. Visibility of the spheres with higher contrasts in the images denoised by BM3D is very close to the ones denoised by BM3D-NCRE. However, the visibility is significantly less for the spheres with lower contrasts. Both BM3D and BM3D-NCRE outperform AIDR3D and nonlocal mean, while the number of visible spheres is almost the same as in the images denoised by nonlocal mean and the images reconstructed by AIDR3D.Figure 6
Low contrast study using CATPHAN low contrast slice. Top: images scanned with 50 mAs/120 kVp and bottom: scanned with 300 mAs/120 kVp. Left to right: image reconstructed with FC52, reconstructed with AIDR3D, reconstructed with FC52 and denoised by BM3D-NCRE, reconstructed with FC52 and denoised by nonlocal mean, and reconstructed with FC52 and denoised by BM3D. In all images window width/window level = 100/70 HU.The final comparison was made using a low dose (50 mAs, 120 kVp) lung CT of a patient reconstructed using FC52 and processed by anisotropic diffusion denoising, BM3D-NCRE, nonlocal mean, and BM3D. A single axial slice of the images is shown in Figure7. As can be seen, anisotropic diffusion removes some fine details and reduces the contrast of the small features. Nonlocal mean, BM3D, and BM3D-NCRE outperform anisotropic diffusion denoising in sense of preserving the small structures. Comparing the results of nonlocal mean with BM3D-NCRE, it can be seen that the low contrast features are kept perfectly unchanged in BM3D-NCRE, while they are removed or their contrasts are reduced in the image denoised by nonlocal mean. In addition, the noise is not homogeneously removed from the image. Comparing BM3D with BM3D-NCRE, it can be seen that the image denoised by BM3D-NCRE has less noise and the removed noise is more homogeneous. It should be noted that the parameters of the anisotropic diffusion denoising are adjusted in such a way that the noise variance in the reconstructed images is the same as the images denoised by BM3D-NCRE.Figure 7
Comparison of the effects of anisotropic diffusion, nonlocal mean, BM3D, and BM3D-NCRE. (a) In the original image, the circular regions show the area from which the noise variance is measured, the dashed rectangular region is the area which is shown in (d)-(e), and average noise of the three regions is around 55 HU. (b) Denoised by anisotropic diffusion, noise is around 25 HU. (c) Denoised by BM3D-NCRE, noise is around 25 HU. (d) Denoised by nonlocal mean, noise is around 21 HU. (e) Denoised by BM3D, noise is around 28 HU. ((f)–(j)) The zoomed-in region shown by dashed rectangle in image (a). The difference between original image and the image denoised by (k) anisotropic diffusion, (l) BM3D-NCRE, (m) nonlocal mean, and (n) BM3d are shown. The window width/window level in (a)–(j) is 1600/−300 HU and is 100/0 HU in (k)–(n).
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)
(k)
(l)
(m)
(n)


## 5. Conclusions

An iterative denoising scheme was proposed for low dose CT images, which adjusts the denoising parameters at each iteration based on the effect of the denoising method in the previous iteration. Noise confidence region evaluation (NCRE) was used to compare the Gaussian noise with denoising residual to determine if the denoising was effectively, weakly, or strongly executed. Based on this information the denoising parameters were adjusted for the next iteration. BM3D was used in the new proposed iterative scheme. The phantom study showed that our proposed method improved low contrast detectability. The patient study demonstrated that the image was efficiently denoised and the visibility of small objects was preserved. However, it should be noted that the modified optimization model is not accurate when the electronic noise dominates the photon fluctuations as could occur for very low doses. In addition, the streak artifacts would still be present in images denoised by the proposed method.



---

*Source: 638568-2015-05-24.xml*