Computational and Mathematical Methods in Medicine

Volume 2015, Article ID 638568, 12 pages

http://dx.doi.org/10.1155/2015/638568

## Adaptively Tuned Iterative Low Dose CT Image Denoising

^{1}Institute of Biomaterials and Biomedical Engineering, University of Toronto, Toronto, ON, Canada M5S 3G9^{2}Joint Department of Medical Imaging, Toronto General Hospital, University Health Network, Toronto, ON, Canada M5G 2N2^{3}Department of Electrical and Computer Engineering, Ryerson University, Toronto, ON, Canada M5B 2K3

Received 23 January 2015; Revised 2 May 2015; Accepted 3 May 2015

Academic Editor: Hugo Palmans

Copyright © 2015 SayedMasoud Hashemi et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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

#### 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 can be modeled bywhere is the noiseless image and is a zero mean additive anisotropic Gaussian noise with variance of , 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,in which is a regularization parameter that controls the tradeoff between the data fidelity and the regularization term and . Different regularization terms, , lead to different denoising methods. For example, in TV-based methods [24–26], where and are the gradient in horizontal and vertical directions, and in wavelet soft thresholding methods , where 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: grouping similar patches with a reference patch, calculation of the 3D wavelet coefficients of each stack of patches, denoising the wavelet coefficients (thresholding in step or Wiener filtering in step ), and 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 . To find best match to similar patches in the step and to reliably estimate the Wiener coefficients, the method requires a reliable estimate of the noiseless image, which is the main purpose of step . 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 is used in step . 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 onin which is the patch grouping information, is the number of 3D patches, denotes the noisy patches, denotes the noiseless 3D patches, and is the noise at each 3D patch. Conventional BM3D uses an independent identical additive Gaussian noise model in which the noise variances are similar in all patches. Using this assumption, its regularization term is a nonlocal wavelet -norm , where [43] in which is the 3D wavelet transform. Using this regularization in (2), BM3D solves the following optimization problem in its first step: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 , 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:where is the mean and is the variance of the projections at the th projection angle () and the th detector bin whose distance is from the detectors center, is a scaling factor, is the electronic noise variance, and 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,the noise variance in the reconstructed images can be described by [38]where is the distance between the center of two adjacent detectors, is the parallel projection at the th angle and the th detector bin, and 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 is symmetric about the center of rotation and 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 . Using (3), the modified optimization problem used in this paper is given byin which is a set of regularization parameters that are functions of . 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

If denotes the signal recovered at th iteration, the denoising residual at the end of th iteration can be expressed by which, ideally, should be the noise (3), . Here, we provide a quantitative measure that verifies the similarity between the structure of and that of .

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

In [47] it was shown that the following function of zero mean white Gaussian noise, , with length , and for any given scalar value of :is equivalent to sorting the absolute value of the noise elements . The expected value of this function is and its variance is , where and is the cumulative distribution function (CDF) of a Gaussian distribution. Therefore, is bounded by the following lower () and upper () values: with probability of . 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 Figure 1, these boundaries divide the space into three regions. At the end of each iteration 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 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 th iteration has removed not only the noise but also parts of the noiseless data itself, will have some of the image information making its samples larger than Gaussian noise. Therefore, for a specific value of , (average number of s with absolute values smaller than ) is smaller than and falls in Region III. This will enforce continuation of the denoising to the th step with changing of the regularization parameters such that the denoising algorithm extracts less noise in the next iteration, that is,* decreasing *, . If falls in Region I, it indicates that the noise is partially removed. In this case the algorithm continues to the th step and changes the regularization parameters such that more noise is extracted by the denoising algorithm, that is,* increasing *, . In summary, at each iteration when is in either Region I or III, the regularization parameter is updated such that it moves toward Region II. The value of can be tuned as a fixed or an adaptively changing variable based on the euclidean distance between and , . In our proposed method to update a global 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 , 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.