Adaptive nonlocal means filtering based on local noise level for CT denoising Zhoubo Li, Lifeng Yu, Joshua D. Trzasko, David S. Lake, Daniel J. Blezek, Joel G. Fletcher, Cynthia H. McCollough, and Armando Manduca Citation: Medical Physics 41, 011908 (2014); doi: 10.1118/1.4851635 View online: http://dx.doi.org/10.1118/1.4851635 View Table of Contents: http://scitation.aip.org/content/aapm/journal/medphys/41/1?ver=pdfcov Published by the American Association of Physicists in Medicine Articles you may be interested in Noise, sampling, and the number of projections in cone-beam CT with a flat-panel detector Med. Phys. 41, 061909 (2014); 10.1118/1.4875688 A simple approach to measure computed tomography (CT) modulation transfer function (MTF) and noise-power spectrum (NPS) using the American College of Radiology (ACR) accreditation phantom Med. Phys. 40, 051907 (2013); 10.1118/1.4800795 Three-dimensional anisotropic adaptive filtering of projection data for noise reduction in cone beam CT Med. Phys. 38, 5896 (2011); 10.1118/1.3633901 Projection space denoising with bilateral filtering and CT noise modeling for dose reduction in CT Med. Phys. 36, 4911 (2009); 10.1118/1.3232004 Adaptive mean filtering for noise reduction in CT polymer gel dosimetry Med. Phys. 35, 344 (2008); 10.1118/1.2818742

Adaptive nonlocal means filtering based on local noise level for CT denoising Zhoubo Li Department of Physiology and Biomedical Engineering, Mayo Clinic, Rochester, Minnesota 55905

Lifeng Yu Department of Radiology, Mayo Clinic, Rochester, Minnesota 55905

Joshua D. Trzasko, David S. Lake, and Daniel J. Blezek Department of Physiology and Biomedical Engineering, Mayo Clinic, Rochester, Minnesota 55905

Joel G. Fletcher and Cynthia H. McCollough Department of Radiology, Mayo Clinic, Rochester, Minnesota 55905

Armando Manducaa) Department of Physiology and Biomedical Engineering, Mayo Clinic, Rochester, Minnesota 55905

(Received 19 June 2013; revised 22 October 2013; accepted for publication 4 December 2013; published 31 December 2013) Purpose: To develop and evaluate an image-domain noise reduction method based on a modified nonlocal means (NLM) algorithm that is adaptive to local noise level of CT images and to implement this method in a time frame consistent with clinical workflow. Methods: A computationally efficient technique for local noise estimation directly from CT images was developed. A forward projection, based on a 2D fan-beam approximation, was used to generate the projection data, with a noise model incorporating the effects of the bowtie filter and automatic exposure control. The noise propagation from projection data to images was analytically derived. The analytical noise map was validated using repeated scans of a phantom. A 3D NLM denoising algorithm was modified to adapt its denoising strength locally based on this noise map. The performance of this adaptive NLM filter was evaluated in phantom studies in terms of in-plane and crossplane high-contrast spatial resolution, noise power spectrum (NPS), subjective low-contrast spatial resolution using the American College of Radiology (ACR) accreditation phantom, and objective low-contrast spatial resolution using a channelized Hotelling model observer (CHO). Graphical processing units (GPU) implementation of this noise map calculation and the adaptive NLM filtering were developed to meet demands of clinical workflow. Adaptive NLM was piloted on lower dose scans in clinical practice. Results: The local noise level estimation matches the noise distribution determined from multiple repetitive scans of a phantom, demonstrated by small variations in the ratio map between the analytical noise map and the one calculated from repeated scans. The phantom studies demonstrated that the adaptive NLM filter can reduce noise substantially without degrading the high-contrast spatial resolution, as illustrated by modulation transfer function and slice sensitivity profile results. The NPS results show that adaptive NLM denoising preserves the shape and peak frequency of the noise power spectrum better than commercial smoothing kernels, and indicate that the spatial resolution at low contrast levels is not significantly degraded. Both the subjective evaluation using the ACR phantom and the objective evaluation on a low-contrast detection task using a CHO model observer demonstrate an improvement on low-contrast performance. The GPU implementation can process and transfer 300 slice images within 5 min. On patient data, the adaptive NLM algorithm provides more effective denoising of CT data throughout a volume than standard NLM, and may allow significant lowering of radiation dose. After a two week pilot study of lower dose CT urography and CT enterography exams, both GI and GU radiology groups elected to proceed with permanent implementation of adaptive NLM in their GI and GU CT practices. Conclusions: This work describes and validates a computationally efficient technique for noise map estimation directly from CT images, and an adaptive NLM filtering based on this noise map, on phantom and patient data. Both the noise map calculation and the adaptive NLM filtering can be performed in times that allow integration with clinical workflow. The adaptive NLM algorithm provides effective denoising of CT data throughout a volume, and may allow significant lowering of radiation dose. © 2014 Author(s). All article content, except where otherwise noted, is licensed under a Creative Commons Attribution 3.0 Unported License. [http://dx.doi.org/10.1118/1.4851635] Key words: CT dose reduction, image denoising, nonlocal means filtering, adaptive denoising, noise estimation 011908-1

Med. Phys. 41 (1), January 2014

0094-2405/2014/41(1)/011908/16

© Author(s) 2014

011908-1

011908-2

Li et al.: Adaptive nonlocal means filtering for CT denoising

1. INTRODUCTION Radiation dose from clinical CT scanning is an increasing health concern worldwide.1 The current guiding principle in CT clinical practice is to use radiation dose levels as low as reasonably achievable while maintaining acceptable diagnostic accuracy. However, lowering radiation dose alone generally produces a noisier image and may seriously degrade diagnostic performance or undermine physician confidence. Many techniques have been proposed for controlling noise in CT, and these can be broadly categorized into three major types: projection space denoising, image space denoising, and iterative reconstruction (IR). Projection space techniques, which work on either the raw projection data or the log-transformed sinogram, attempt to reduce noise in the projection data domain prior to image reconstruction.2–8 In general, these techniques have the advantage that noise properties in projection space are fairly well understood. One potential drawback of projection-based methods is that they result in some loss of image sharpness due to the fact that edges in projection data are not welldefined. Image-space denoising involves applying linear or nonlinear filters directly to the reconstructed images. Most such techniques [e.g., bilateral filtering,9 total variation denoising,10 nonlocal means (NLM) denoising (Fig. 1),11 and k-SVD denoising12 ] take advantage of the strong structural and statistical properties of objects in image space (e.g., sharp edges, similarities between neighboring pixels). In CT, they can be implemented directly and without access to the raw data. However, CT noise in image space is difficult to model accurately and has strong spatial variations and correlations. It can therefore be more difficult for such techniques to achieve an optimal tradeoff between denoising and blurring or artifacts, or to get consistent performance across an entire scan volume. IR techniques are more accurately considered reconstruction rather than denoising techniques, and take advantage of statistical assumptions about the noise properties in projection space, prior information in image space, and CT scanner system models.13–15 IR techniques require access to the raw data and accurate knowledge of the details of the scanner geometry, photon statistics, data-acquisition, and correction physics, thus highly dependent on specific scanner models. True IR is very computationally intensive (e.g., up to several hours per data set), which has prevented fast clinical application to date, although software methods16, 17 and hardware methods18–21 have been investigated to accelerate the iterative procedure. Due to the extremely high computational load of true IR, hybrid techniques have recently been developed that attempt to gain many of the benefits of true IR with much lower computational load (e.g., Sinogram Affirmed Iterative Reconstruction, SAFIRE from Siemens). Some of these are now available commercially, but there is little data to demonstrate differences in performance between approaches. Despite the development of various IR methods, they are only available on some of the latest scanner models from each manufacturer, and are implemented on a CT scanner image reconstruction Medical Physics, Vol. 41, No. 1, January 2014

011908-2

system—so purchase of this hardware and software is on a per scanner basis and very expensive. Most hospitals and departments of radiology have a heterogeneous mix of CT scanner system makes and models, and many older generation CT systems do not have commercial alternatives for upgrading their image reconstruction system with IR methods. Our goal was to develop a denoising strategy that can be broadly used across the CT community over a heterogeneous scanner fleet that encompasses different manufacturers as well as different models of varying age and software revision. These requirements lead us to consider image space denoising techniques, which are relatively simple to implement, work on the image data alone, and can be applied retrospectively. As mentioned above, it is difficult for image space techniques to model CT noise or scanner details accurately, and thus they may appear to necessarily be at a disadvantage with respect to projection space or IR methods. However, the spatial structure models in modern image denoising algorithms are significantly more advanced than the spatial regularization terms that are incorporated in current IR methods. It is thus not at all clear that image space results will necessarily be inferior in terms of noise reduction. NLM denoising11 is an effective image denoising strategy that exploits the inherent redundant information present in most images. NLM generalizes the notion of finite spatial differences and utilizes a measure of similarity between nearby image patches to estimate underlying image structure. This allows NLM to preserve a high degree of image texture and fine detail. We previously incorporated NLM denoising into the clinical CT workflow using both algorithmic and hardware speedups, such that images are returned to viewing workstations in under 5 min.22 However, the standard NLM algorithm uses a uniform filtering strength to denoise the image, while in CT images the noise level varies significantly within and across slices. Therefore, applying NLM filtering to CT images using a single global filtering strength cannot achieve optimal denoising performance (Fig. 2). In this work, we have developed a technique for efficiently estimating the local noise level for CT images by creating a map of local noise distribution in the image space, and have modified the NLM algorithm to adapt to local variations in noise levels. This paper is organized as follows. First, we describe the principles of adapting NLM denoising to local noise level. Second, we present a simple method that can efficiently and accurately estimate the local noise (i.e., calculate the noise map). Third, an accelerated implementation of the method based on graphical processing units (GPU) is described. Finally, we present results from phantom and patient evaluations.

2. METHODS 2.A. Adapting NLM to local noise level

NLM assumes that images contain a substantial amount of redundant local structure, and this property can be exploited

011908-3

Li et al.: Adaptive nonlocal means filtering for CT denoising

011908-3

F IG . 1. Closeup of an original CT slice of the abdomen (left) and the image after NLM denoising (right). Note improvement in conspicuity of hepatic metastases and appearance of liver parenchyma.

F IG . 2. Images of a CT slice of the abdomen (first column), after standard NLM denoising with low strength (second column), high strength (third column), and with adaptive NLM (fourth column). Top row: A zoomed version of one region of the images. The standard deviations inside ROI 1 are 37.57, 20.38, 14.49, and 20.12 HU, respectively. Middle row: A zoomed version of another region of the images. The standard deviations inside ROI 2 are 48.08, 31.44, 20.73, and 20.50 HU, respectively. Bottom row: Difference between filtered results and original images (zoomed version of the same region as in the second row). Note standard NLM denoising based on a single noise level may be too strong in some regions (ROI 1, third column) or too weak in other regions (ROI 2, second column). Conversely, adaptive NLM denoising based on an estimated noise map can achieve more uniform noise reduction in different regions with both high and low noise levels. No obvious local structures exist in the difference images. Medical Physics, Vol. 41, No. 1, January 2014

011908-4

Li et al.: Adaptive nonlocal means filtering for CT denoising

to reduce noise by performing appropriately weighted averages of pixel intensities. The NLM filtering process can be expressed as  w(s, t)v(t) t∈W , (1) u(s) =  w(s, t) t∈W

where u(s) is the filtered pixel intensity at pixel s and v(t) is the pixel intensity at pixel t in the original image. W represents a window (termed a search window) centered at pixel s, and w(s,t) is the weight assigned to an individual pixel t within this window when filtering pixel s. In NLM, this weight is

 u(x0 , y 0 , z0 ) =

v(x, y, z)e

(x,y,z)∈W



e

−

1 Nh2

−

1 Nh2



011908-4

⎛  ⎜ δ∈P w(s, t) = exp ⎝−

Gσ (δ) [v(s + δ) − v(t + δ)]2 h2

⎞ ⎟ ⎠

(2)

representing the similarity between two patches of the same size P centered at pixels s and t. Gσ is a Gaussian kernel of variance σ 2 , and δ represents the relative offsets of pixels inside a patch. h is a smoothing parameter used to control the amount of denoising and is usually taken to be proportional to the assumed or known noise level. Standard NLM assumes that noise in the image is Gaussian distributed and spatially invariant and thus h is a global constant. In our case, we extend the NLM concept to 3D images and the filtering process can be expressed as

[v(x0 +δx ,y 0 +δy ,z0 +δz )−v(x+δx ,y+δy ,z+δz )]2

(δx ,δx ,δx )∈P



[v(x0 +δx ,y 0 +δy ,z0 +δz )−v(x+δx ,y+δy ,z+δz )]2

,

(3)

(δx ,δx ,δx )∈P

(x,y,z)∈W

where pixels are specified by their 3D coordinates, (x, y, z). Both the search window and the patch are 3D blocks. N is the number of pixels in the patch and we have added this to the formula to normalize the sum squared difference calculation to express average sum squared difference per pixel, allowing a more meaningful expression of h in terms of the noise level. As suggested in Ref. 23, the Gaussian kernel in the original formula can be removed without noticeable differences in the filtered result, allowing for a more efficient precalculation of sum squared differences. The original NLM, although noniterative, usually has a high computational cost due to the computation of similarities between neighborhoods in a large search window region. In our case, we extend the NLM concept to full 3D blocks in the data volume, further increasing the computational load. To meet the demand of clinical flow on the processing speed, we adopt a 15 × 15 × 15 search window and a 3 × 3 × 3 patch size for all examples in this work. We have implemented and improved the precalculation of sum squared differences described in Refs. 23 and 24, increasing the computational speed 20-fold. We have also reduced computation time by an additional 35% by taking advantage of symmetry between weighting factors when two specific pixels switch roles of s and t in Eq. (2), and 30% by precalculating the exponential function in a lookup table. Finally, porting the algorithm to GPUs (Ref. 22) gave an additional 35× speed gain. In CT volumes, the noise level varies within and across slices, often by 2× within a slice and as much as 3× across slices. This implies that NLM denoising based on a single noise level may be too weak in some places (accomplishing little), too strong in others (blurring fine detail), or both (Fig. 2). It is therefore desirable to modify the NLM algorithm to adapt to the local noise level. Here, we have adjusted the Medical Physics, Vol. 41, No. 1, January 2014

strength of h = k*SD locally, where SD is the estimated noise level of the pixel being denoised and k is a proportionality factor that serves to tune the denoising strength. This can be easily integrated with the speedup technique described above and involves little additional computational effort. We maintain a single k throughout the whole volume so that h varies proportionally to the local SD. We also introduce a blending factor B, representing an averaging of the denoised image (with fraction B) with the original image (with fraction 1-B). The Appendix gives pseudocode of our algorithm for fast 3D NLM filtering. In all the examples in this work, k was determined empirically so that the noise level in the filtered image achieved a predetermined value (e.g., to match the noise level of a given commercial kernel). However, this adaptive NLM filter requires a map of the local noise level, which in turn requires developing a way to efficiently estimate such a map, from image data alone, in a time compatible with the clinical workflow. 2.B. Noise map estimation

The CT image noise distribution can be estimated or calculated using many different approaches. Repeating scans multiple times for the same object is ideal but impractical to implement on patients. Another approach is through Monte Carlo simulation, adding simulated noise to raw data and reconstructing multiple realizations of CT images. However, this method requires access to the raw data, and a large number of repeated noise simulations and reconstructions, which is time consuming and makes it difficult to meet the clinical workflow requirement. A more elegant approach is to derive the noise distribution analytically by propagating a noise

011908-5

Li et al.: Adaptive nonlocal means filtering for CT denoising

model through the reconstruction equations. This has been implemented in simple fan-beam CT,25, 26 where variance and covariance at each location of image space can be analytically derived, assuming a simple CT noise model. In modern CT scanners, helical scans with multidetector rows are typically employed rather than a simple 2D sequential fan-beam mode. The derivation of an analytical formula of noise (variance and covariance) in reconstructed images requires an accurate knowledge of reconstruction algorithms, which in modern scanners typically involve a rebinning process to convert cone-beam data to quasiparallel-beam data and a weighted 3D filtered backprojection (FBP) process.27 Due to the complicated numerical operations in the reconstruction process, an accurate and fast derivation of noise in the final image may be difficult. However, for the purpose of image-based NLM denoising adaptive to local noise level, only an approximate noise level is needed, and it may be suffice to assume a simple CT geometry and reconstruction process. In this paper, we propose to use a simple technique, based on estimating local noise pixelwise in each 2D slice separately with a fan beam approximation, which does not require access to the raw data, is computationally efficient, and is easily parallelizable for speed. The basic process is described below:

011908-5

We therefore use a 2D fan-beam geometry, with fan-angle and focal length consistent with the clinical scanners tested in this study. A standard ray-driven forward projection method was employed to generate the CT sinogram. 2.B.3. Noise modeling in sinogram data

Although modern CT detectors are not photon-counting elements, but rather energy integrators that generate a signal proportional to the total energy deposited in the detector, a photon-counting model is still a good approximation of quantum noise and is widely used for characterizing noise properties of CT data.15, 26, 29 In our simulation, after obtaining the CT sinogram, we calculated the transmitted photon map by dividing the incident photon map by the exponential of the sinogram. Assuming a Poisson distribution, the variance of the measured integral attenuation map, p, can be calculated by30 Var{p(kτ, θ )} =

1 . ¯ N (kτ, θ )

(5)

2.B.1. Calculation of linear attenuation coefficient

Here, Var{p(kτ , θ )} is the variance of the integral attenuation map at the detector index of k when the projection angle is θ . τ is the sampling interval and is equal to the interval between adjacent detectors and N (kτ, θ ) is the mean value of the detected photon numbers. The incident number of photons varies for each projection angle due to the use of the AEC technique,31 and is also nonuniform across the radiation field of the x-ray beam, mainly due to the use of the beam-shaping bowtie filter.32 The effect of the bowtie filter can be characterized by measuring a map of noise-equivalent number of photons along the detector row from a set of air scans.33 The tube current modulation can be estimated based on the attenuation level along each projection angle and modulation strategy described in Ref. 34. Therefore, the incident number of photons is a function of both detector bin index and projection angle. Details on how to incorporate the effects of AEC and the bowtie filter, and how to quantify the incident number of photons have been described previously.8

Pixel intensities in CT images are designated in Hounsfield units (HU), which can be converted into the linear attenuation coefficient of corresponding tissues by the equation

2.B.4. Analytical calculation of noise in reconstructed images

i. Calculate the linear attenuation coefficient from the CT image. ii. Generate CT sinogram data using a fan-beam CT geometry. iii. Estimate the noise distribution of the sinogram data, incorporating the effect of the bowtie filter and automatic exposure control (AEC). iv. Apply the analytical formula to reconstruct the noise map in the final reconstructed images. v. Scale the obtained noise map to match the actual local noise level in CT images. Below we explain each step in more detail.

μ=

(CT + 1000) × μw . 1000

(4)

Here, μ refers to the linear attenuation of material with units of cm−1 . The value of μw for different x-ray photon energies can be obtained from the NIST X-ray Mass Attenuation Coefficients database.28 Here, we use a monoenergetic approximation of the tube spectrum. 2.B.2. Generation of approximate CT sinogram

In principle, calculation of the sinogram requires an accurate knowledge of the CT acquisition geometry. Here, the noise map estimation must be fast, but can be approximate. Medical Physics, Vol. 41, No. 1, January 2014

In our current implementation, we calculated the analytical noise map based on a simple 2D fan-beam geometry and a rebinning FBP reconstruction. For further simplification and efficiency, the correlation introduced in the rebinning step was also neglected. As will be demonstrated below, these simplifications still yield a reasonably accurate noise map estimate, which can be implemented very efficiently, as required for this technique to be clinically viable. Once the noise variance in the projection space is determined, the noise distribution of each pixel in the CT image can be obtained by analytically propagating the noise variance through the reconstruction algorithms. The detailed theoretical description of the parallel-beam filtered backprojection reconstruction can be found in Ref. 35.

Li et al.: Adaptive nonlocal means filtering for CT denoising

011908-6

In order to obtain the variance expression of the reconstructed image in a discrete format, a discrete version of the FBP reconstruction has to be used. Here, we assume a linear interpolation in the backprojection, which is given by π [(1 − η(x, y, θn )Qθn (kτ ) Np n=1 +η(x, y, θn )Qθn (kτ + τ )],

I /2

h(kτ − k  τ )pθ (k  τ ),

(7)

k  =−I /2

where I+1 denote the number of detector bins and h(kτ ) represents the discrete and spatial version of the ramp or Ram-Lak filter and can be expressed as ⎧ 2 ⎪ ⎨ 1/4τ , k = 0 0, k even . h(kτ ) = (8) ⎪ ⎩ − 1 , k odd k2 π 2 τ 2 The above filtered backprojection reconstruction process can be divided into three steps: filtering, interpolation, and backprojection. All three processes can be treated as linear processes, and the variance for a linear combination of random variables Xi can be expressed as36   Var ak Xk = ak2 Var(Xk ) k

k

+

k

2ak ak1 Cov(Xk , Xk1 ).

(9)

k1 =k

After ramp filtering, the variance of the filtered projections can thus be written as Var(Qθ (kτ )) = τ 2

I /2 h2 (kτ − k  τ ) . N¯ (kτ, θ ) k  =I /2

(10)

After filtering, the projection data are interpolated and backprojected into image space and the variance can be expressed as Np πτ2 Var(f (x, y)) = [((1 − η(x, y, θn ))2 Var(Qθn (kτ )) Np n=1

+ η2 (x, y, θn )Var(Qθn (kτ + τ )) + 2η(x, y, θn )(1 − η(x, y, θn ) × Cov(Qθn (kτ ), Qθn (kτ + τ ))], Medical Physics, Vol. 41, No. 1, January 2014

Cov(Qθn (kτ ), Qθn (kτ + τ )) =

I /2

h(kτ )h(kτ + k  τ )Var2 (Qθn (kτ )).

(12)

2.B.5. Scaling to match the noise level in CT images

(6)

where η(x, y, θ n ) ∈ [0, 1] defines the linear interpolation weight between “filtered” projection data, Qθ (kτ ), and the adjacent detector, Np number of projection angles, θ n ∈ [0, π ], and Qθ (kτ ) = τ

where the covariance between the adjacent detectors can be estimated by

k  =−I /2

Np

f (x, y) =

011908-6

(11)

For proper denoising, the raw analytical noise map, which has arbitrary scaling due to the lack of knowledge of the actual number of input photons, must be scaled to match the correct noise level in the CT image. To do this, the mean and standard deviation within 5 × 5 neighborhoods are calculated throughout the CT image. The standard deviations and the corresponding noise map values are compared for those points whose mean value falls within the range of −200 to +400 HU, and the relative scaling between the two that minimizes the summed-squared-difference between the two is determined by a golden-section-search optimizer.37 In homogeneous tissue, the measured standard deviations would be a noisy but accurate estimate of local noise levels. However, since real CT images may contain tissue heterogeneity or image structures, the calculated standard deviation (i.e., “noise”) values are often too high, reflecting both actual noise and underlying heterogeneity. To minimize this bias, only scaled noise map values that fall within ±10 HU of the corresponding sampled standard deviations are used in the summed-squared-difference calculations. In this way, erroneously high noise estimates are rejected from the difference calculation during the optimization. The resulting scale value yields the best overall match between the (scaled) noise map and the estimated noise level in the CT image, without manual intervention. 2.C. Validation of analytical noise estimate by multiple repetitive scan of a phantom

We validated the analytical noise estimate by comparing the results with a reference standard. For determination of a reference noise level distribution in CT images, we repeatedly scanned an anthropomorphic thoracic phantom 100 times with both an axial scan mode and a helical scan mode (Definition-Flash, Siemens Healthcare, Forchheim, Germany). Only a single source was used for both scan modes. A routine abdominal protocol was used. For axial scan, the parameter setting was: 120 kV, 0.5 s rotation time, 64 × 0.6 mm detector collimation, 237 mAs. For helical scan, the parameter setting was: 120 kV, 0.5 s rotation time, 64 × 0.6 mm detector collimation with z-flying focal spot, effective mAs (mAs/pitch) of 118, and a helical pitch of 0.8. The corresponding half dose scans were also acquired to investigate the noise properties of projection data. All images were reconstructed using a B40 kernel with a slice thickness of 1 mm and a FOV of 50 cm. The reference noise map was determined by calculating the standard deviation of the CT number of the same pixel from 100 repeat scans.

011908-7

Li et al.: Adaptive nonlocal means filtering for CT denoising

2.D. Image quality evaluation of the adaptive NLM method

We performed phantom studies to evaluate the in-plane and cross-plane high-contrast spatial resolution, noise power spectrum (NPS), and low-contrast resolution of the adaptive NLM filter. The in-plane high-contrast spatial resolution was evaluated using the modulation transfer function (MTF), which was measured using a 0.2 mm thin wire. The scanning parameters were: 120 kVp, 0.5 s rotation time, 64 × 0.6 mm detector collimation with a z-flying focal spot, helical pitch of 0.8, and an effective mAs of 120. Images were reconstructed using a B40 kernel with a slice thickness of 0.6 mm and a FOV size of 50 mm, and also denoised using the adaptive NLM filter, with strength set to achieve 30% noise reduction. The MTF for both original and denoised images was calculated using a previously described method.38 The cross-plane high-contrast spatial resolution was evaluated using the slice sensitivity profile (SSP), which was measured using a thin gold foil phantom (25 μm thickness) that was embedded inside tissue-equivalent plastic cylinder with a diameter of 23 mm (QRM, Moehrendorf, Germany). The scanning parameters are: 120 kVp, 0.5 s rotation time, 64 × 0.6 mm detector collimation with a z-flying focal spot, helical pitch of 0.8, a mAs value of 240. Images were reconstructed with a kernel of B40 and a nominal slice thickness of 1 mm with a slice interval of 0.1 mm, and also denoised using the adaptive NLM filter, again with strength set to achieve 30% noise reduction. For each image, we subtracted the background mean and recorded the maximum CT number within a region of interest (ROI) centered over the gold foil. The SSP was plotted as the normalized CT numbers as a function of slice location. The NPS was measured using a 30 cm diameter cylindrical phantom filled with water. The phantom was scanned using following parameters: 120 kV, 64 × 0.6 mm detector collimation, 0.5 s rotation time, a helical pitch of 0.8. Images were reconstructed with commercial kernels of B40 and B20, a FOV size of 307 mm, a slice thickness of 0.6 mm, and an interval of 0.6 mm. 100 images were obtained for each reconstruction. The B40 image was denoised using the adaptive NLM filter. The strength setting of the adaptive NLM was adjusted, both with and without a blending factor of 0.5, such that the noise level measured in the center of the image matched that in the B20 image. A procedure following the general framework outlined in Ref. 39 was used to calculate a fully 3D NPS for each of the kernel setting and the NLM denoising. A circular average over a single 2D slice of the 3D NPS in the central axial plane yields a 1D NPS profile. The low-contrast resolution was evaluated both subjectively using the low-contrast resolution module of a phantom used for American College of Radiology (ACR) accreditation and objectively using a validated channelized Hotelling observer (CHO).40 For the subjective evaluation, the low-contrast module of the ACR phantom was scanned with three dose levels at mAs values of 240, 180, and 120. Other scanning parameMedical Physics, Vol. 41, No. 1, January 2014

011908-7

ters were 120 kVp, 0.5 s rotation time, 64 × 0.6 mm detector collimation with a z-flying focal spot, helical pitch of 0.8. Images were reconstructed using a B40 kernel with a slice thickness of 5 mm. For denoising, a thin slice thickness of 1 mm was first reconstructed. After denoising (both with and without blending, with strength set to achieve 30% noise reduction), images with a slice thickness of 5 mm were created. The original and the denoised images were compared to evaluate the effect of the filter on low contrast resolution. For the objective evaluation of low-contrast performance, we used the same experimental setup and CHO configurations as in our previous study, where a CHO model observer was demonstrated to be able to predict human observer performance in a 2-alternative forced choice (2AFC) lesion detection task.40 In that study, three rods (3, 5, 9 mm in diameters) with a CT number of −15 HU at 120 kV were placed in a 35 × 26 cm water tank to simulate lesions with three different sizes. The phantom was scanned 100 times at each of five dose levels (CTDIvol = 2.8, 5.7, 11.4, 17.1, and 22.8 mGy). Twenty-one 2AFC studies were created, including 15 for FBP (five dose levels × three lesion sizes) and six for IR (two dose levels × three sizes). A CHO with Gabor channels was used to predict the percent correct for each 2AFC task. The performance predicted by the CHO was compared with that obtained by four medical physicists for both FBP and IR for the same tasks. The results showed an excellent agreement between human and CHO model observers. In the current evaluation study, we used the same experimental setup and CHO configurations to predict the performance of adaptive NLM filtering (with strength set to achieve 30% noise reduction) for six 2AFC tasks (two dose levels × three sizes). The performance measure, percent correct in the 2AFC task, was compared between the original FBP-based reconstruction and adaptive NLM filtering.

2.E. Computation time

Noise map estimation and NLM filtering (particularly in 3D) are computationally intensive, normally requiring a computation time far beyond clinically acceptable timeframes. Using general purpose GPU, algorithm performance was reduced from more than 60 min for a 300 slice series to less than 5 min.22 The GPU hardware consists of a server chassis with eight Nvidia GTX 680 graphics cards. The server has eight Intel Xeon CPU cores running at 2.4 GHz with 24 G of RAM. Each 680 card has 1536 CUDA cores operating at 1 GHz, with 2 G of GDDR5 RAM. The implementation of adaptive NLM filtering is broken into two logical steps: noise map estimation and NLM filtering. Both steps operate in a hybrid mode, utilizing both CPU and GPU resources to best leverage the relative strengths of each platform. Noise map estimation is slice-by-slice and naturally lends itself to a parallel implementation with a final serial step of scaling. First the image volume is divided into sections and distributed onto each GPU on the server. The images are loaded into fast texture memory. Second, a back-

011908-8

Li et al.: Adaptive nonlocal means filtering for CT denoising

projection step is performed using the texture map hardware to speed up processing. Next, the forward projection step is performed using a combination of GPU and CPU processing. Finally, the resulting noise map estimate from each GPU card is recombined into a single volume. Scaling is performed in the CPU and the final result is transferred to the NLM filtering. NLM filtering is performed entirely on the GPU. The only responsibilities of the CPU are to divide the volume and noise map into sections, configure the GPU processing kernels, and collect the final output. The first step is to divide the image and noise map volumes into sections sized for processing on each GPU. Sizing is critical because GPUs do not gracefully handle out-of-memory situations, and the algorithm uses many double precision buffers to cache intermediate results. Second, the volume data are transferred to the GPU. Next, the NLM algorithm is implemented in a “brute-force” manner. Each GPU core is responsible for processing a single pixel at a time. The CPU orchestrates the GPU to accumulate intermediate results for each search location. Finally, the intermediate results are collected and scaled according to the noise map estimation. Image data are recombined by the CPU and saved in DICOM format. We performed a speed test by processing a series of CT images (n = 64) on an eight-GPU server and recording the total computing and transfer time for each case. The mean and standard deviation of this total processing time per slice were calculated.

2.F. Clinical assessment and integration into CT practice

Radiologist feedback was obtained by comparing routine dose CT images to lower dose CT images (obtained using a validated noise-insertion algorithm),41 which were reconstructed with routine FBP and sequentially denoised using adaptive NLM. After positive evaluation for a variety of focal liver, renal, and bowel pathologies, the decision was made to undertake a pilot study in the clinical practice. Lower dose CT enterography and CT urography protocols were developed based on our experience with image quality improvements with commercial iterative reconstruction on late-model CT scanners, and were developed for two 64 slice MDCT scanners (i.e., GE VCT 64, General Electric Healthcare; Sensation 64; Siemens Healthcare) without denoising options in our heterogeneous CT practice. After selection of denoising strength (k = 0.6, B = 0.5) by our radiologist coauthor and feedback from GI and GU radiologists, a server containing GPUs with the adaptive NLM algorithm was connected to each CT scanner for potential clinical use. During a two week trial period, lower dose CT enterography and urography images were sent to this computer server and denoised with the adaptive NLM algorithm. At the end of the two week trial period, the GI and GU radiology groups were asked to decide as a group if they wished to integrate denoised images with adaptive NLM into their clinical workflow. Medical Physics, Vol. 41, No. 1, January 2014

011908-8

3. RESULTS 3.A. Validation of analytical noise estimate

The estimated noise map based on the CT image data is shown and compared with the results obtained from the reference noise map from 100 repetitive scans of the same object in Fig. 3. For both axial and helical mode CT images, simply performing a 2D forward projection on individual slices of the image (ignoring the true 3D cone-beam nature of the acquisition), determining the noise variance in the sinogram based on photon statistics that incorporate some of the important physical effects (in particular, the bowtie filter and automatic tube current modulation), and analytically propagating the noise through the reconstruction algorithm can lead to a fairly close approximation to the true noise map. The black curves at the edge of the body and streaks outside the body are artifacts of the clipping of values below −1024 HU in the scanner reconstruction and are not relevant. The bright structure in the reference noise map of the helical scan may be due to the nonrepeatability of some factors in helical scans such as starting angles, which yield slightly different interpolation errors and helical and cone-beam artifacts among repeated scans, especially around sharp boundaries. Because of this effect, the noise estimate from repeated helical scans is not perfect around sharp boundaries. However, the errors around sharp boundaries should not affect the noise estimate in smooth areas obtained by using the repeated-scan method. The ratio maps between the reference noise maps and the scaled analytical noise maps showed little noticeable structure and the pixel standard deviations in the center ROIs are very small. The visual comparison of the 1D line profile from the reference and analytical noise maps also showed reasonable agreement in the shape of local noise distribution.

3.B. Image quality evaluation of the adaptive NLM filtering

3.B.1. High-contrast spatial resolution

To study the effects of filtering on spatial resolution, the MTF curves for the line-wire images obtained with the original B40 reconstruction and after denoising of this image with adaptive NLM are shown in Fig. 4. The dotted and solid curves calculated from the baseline and the filtered images have the same frequency content and consequently similar spatial resolution properties. The calculated SSPs for both baseline and filtered images are calculated and shown in Fig. 5. The SSPs before and after the adaptive NLM filtering do not change, demonstrating that the filter has negligible effects on the cross-plane spatial resolution. 3.B.2. Noise and NPS

Figure 6 shows 1D NPS profiles obtained from the 3D NPS calculation for two original FBP kernels (B40 and B20) and the adaptive NLM applied to the B40 reconstruction (strength set to match the B20 noise level, with and without blending)

011908-9

Li et al.: Adaptive nonlocal means filtering for CT denoising

011908-9

F IG . 3. Reference noise map of 100 repetitive CT scans (first column; top: axial mode; bottom: helical mode) and noise maps estimated by analytical noise propagation through the reconstruction algorithm (second column). The third column shows the ratio between the estimated/reference noise map. The pixel values of the center ROIs in the ratio images are 1.004 ± 0.087 (axial) and 1.006 ± 0.092 (helical), respectively. The fourth column shows that the line profiles across the reference (jagged lines) and estimated (smooth lines) noise maps show reasonable agreement (y axes: CT number, x axes: pixel coordinates).

using the analytical noise map calculation. The noise levels measured in a region of interest (200 mm2 ) in the center of the image were 43.9, 30.1, 30.8, and 29.8 HU, respectively. The spatial frequencies corresponding to the maximum magnitude of the NPS were 2.50, 1.92, 2.25, and 2.00 cm−1 , respectively. Thus, with a noise reduction of 35% compared to B40, adaptive NLM denoising (especially with blending) preserves the shape of the noise power spectrum and the peak frequency better than the B20 commercial smoothing kernel.

3.B.3. Subjective evaluation of low-contrast resolution with ACR phantom

Figure 7 shows the images of the ACR phantom with three dose levels [240 (left), 180 (center), and 120 (right) mAs] after B40 reconstruction (top row) and the images after adaptive NLM filtering without (center row) and with (bottom row) blending. The detection of low contrast objects was improved with the noise reduction, and the mean CT numbers were well-maintained for all three dose levels. 3.B.4. Objective evaluation of low-contrast performance using CHO

The performance predicted by the CHO for detecting lowcontrast lesions in six 2AFC tasks (two dose levels × three lesion sizes) is shown in Fig. 8. No improvement was predicted by CHO for the 3 mm lesion. The performance for the 9 mm lesion was already approaching 100%, so no improvement is observed there either. However, the adaptive NLM filter had a significant improvement over FBP for the 5 mm lesion (FBP: 85.0 ± 3.6% and 90.9% ± 2.8%; NLM: 89.7% ± 2.8%; and 96.9% ± 1.7%; p < 0.01). 3.B.5. Computation time trials

F IG . 4. In-plane noise-resolution properties of adaptive NLM. Adaptive NLM filtering can achieve noise reduction without changing the MTF curve. Medical Physics, Vol. 41, No. 1, January 2014

For 64 cases with varying slice numbers (ranging from 171 to 669), the average processing time per slice is 0.79 ± 0.09 s, which demonstrates that a 300 slice data can be processed and returned in 5 min. The individual step processing times are 0.28 ± 0.06, 0.12 ± 0.01, and 0.38 ± 0.08 s per slice

011908-10

Li et al.: Adaptive nonlocal means filtering for CT denoising

011908-10

F IG . 5. Cross-plane noise-resolution properties of adaptive NLM. Adaptive NLM filtering can achieve noise reduction without changing the cross-plane slice sensitivity profile.

for noise estimation, scaling, and NLM filtering, respectively. GPU technology is advancing at a rapid pace, so future generations of GPUs will further reduce this turnaround time.

3.C. Clinical patient examples and integration into clinical practice

Figure 2 (fourth column) shows that the adaptive NLM denoising based on the local noise map is effective in two different areas of the same slice which have markedly different noise levels. In comparison, standard NLM denoising based on a fixed low strength is less effective in some regions (such as liver), while denoising at a fixed high strength blurs inter-

F IG . 6. Noise power spectrum analysis of two commercial reconstruction kernels and adaptive NLM. Compared to the commercial reconstruction kernels, adaptive NLM (especially with blending) can reduce image noise while better preserving the shape of the noise power spectrum and the peak frequency. Medical Physics, Vol. 41, No. 1, January 2014

nal anatomic detail features in other regions. Adaptive NLM changes denoising strength according to local noise characteristics to result in a more uniform appearance. This adaptive noise reduction capability continues to hold true throughout a CT volume and across patients. Figures 9 and 10 display an A/P routine scan using a SIEMENS Definition-Flash scanner with the following exam parameters: 100 kVp, 319 quality reference mAs, 0.5 s rotation time, 64 × 0.6 mm2 detector collimation, and a helical pitch of 1.0. The volume CT dose index (CTDIvol ) is 10.41 mGy. The images were reconstructed using a commercial reconstruction kernel B40 with a slice thickness of 1 mm. Half and quarter dose images were simulated following a noise insertion method described elsewhere.41 The simulated low dose images then were processed by the adaptive NLM filter and compared with the original full dose ones. The images show that the adaptive NLM filtered low dose images have an image quality typical of higher doses and superior to the low dose images without denoising. The improvements in image quality are illustrative in the detection of lesion and vascular structures in the liver, both intraslice (axial view, Fig. 9) and interslice (coronal view, Fig. 10). Figure 11 shows another A/P routine scan example from a Siemens scanner. The patient was scanned with the following parameters: 120 kVp, 240 quality reference mAs, 0.5 s rotation time, 64 × 0.6 mm2 detector collimation, and a helical pitch of 0.8. The volume CT dose index (CTDIvol ) is 17.1 mGy. The images were reconstructed using a commercial reconstruction kernel B40 with a slice thickness of 5 mm. Half dose axial images were simulated using the validated noise insertion technique,41 and then processed by the adaptive NLM filter and compared with the original images. The filtered images showed conspicuity of a small subcentimeter metastasis similar to the routine dose image, with improved

011908-11

Li et al.: Adaptive nonlocal means filtering for CT denoising

011908-11

F IG . 7. Images for the ACR phantom after B40 reconstruction (top row) and B40 reconstruction followed by adaptive NLM filtering without (center row) and with (bottom row) blending. Note the enhancement of low contrast detection for 240 and 180 mAs, the noise suppression in the adaptive NLM filtered images for all three dose levels, and the well-maintained mean CT number in ROI 1 and 2 after filtering for all three dose levels. From high dose to low dose, the mean/std of CT number in ROI 1 are 96.5/5.3, 96.3/5.8, and 96.4/7.0, respectively. For ROI 2, the values are 90.4/6.0, 90.3/6.7, and 90.3/8.5. After adaptive NLM filtering without blending, the CT number values in ROI 1 are 95.3/3.8, 95.1/4.1, and 95.4/5.0 and for ROI 2 are 89.4/4.4, 89.1/5.1, and 89.4/6.6, respectively. After adaptive NLM filtering with blending, the CT number values in ROI 1 are 95.7/3.8, 95.3/4.1, and 95.7/5.0 and for ROI 2 are 89.8/4.4, 89.3/5.0, and 89.6/6.5, respectively.

image quality in the hepatic parenchyma compared to the half dose image reconstructed with FBP.

Following selection of denoising strength for lower dose CT enterography and CT urography and implementation of adaptive NLM on a computer server, lower dose CT enterography and CT urography protocols underwent routine adaptive NLM denoising, with radiologists evaluating FBP and adaptive NLM images in clinical practice. Following this two week trial, both GI and GU radiology groups unanimously elected to continue to use adaptive NLM in their practice. An example of one lower dose CT enterography from this assessment is shown in Fig. 12, illustrating similar image quality and bowel wall conspicuity with 73% dose reduction. Formal observer performance studies are under way for liver and bowel CT, and future studies are planned.

4. DISCUSSION F IG . 8. Objective evaluation of low-contrast performance using channelized Hotelling model observer. Medical Physics, Vol. 41, No. 1, January 2014

Despite the potential benefit of IR on improving image quality and reducing radiation dose in CT, implementation

011908-12

Li et al.: Adaptive nonlocal means filtering for CT denoising

011908-12

F IG . 9. Effect of adaptive NLM filtering and lower dose levels on appearance of a neuroendocrine tumor metastasis in the liver (axial view). Top row: Simulated 1/4 dose image (left) and simulated 1/2 dose image (right). Bottom row: Simulated 1/4 dose after adaptive NLM (left), simulated 1/2 dose after adaptive NLM (center), and original full dose (right).

of a true IR (e.g, MBIR from GE) remains very challenging in clinical practice mainly due to its long computational time.42 To improve computational speed, various hybrid IR methods have been developed and implemented in some of the latest scanner models by all major CT manufacturers, such

as SAFIRE from Siemens.43 However, most IR or noise reduction methods are only available on relatively newer scanner models, serve only the CT system on which they reside, and image quality improvement and the potential radiation dose reduction remain to be evaluated. In this work, we have

F IG . 10. Coronal images from same patient as prior figure, showing benefit of adaptive NLM filtering of clinical A/P scan in a patient with neuroendocrine metastasis to tumor (coronal images shown). Top row: Simulated 1/4 and 1/2 dose images, respectively. Bottom row: Simulated 1/4 dose image after adaptive NLM (left), 1/2 dose image after adaptive NLM (center), routine dose image with FBP. Medical Physics, Vol. 41, No. 1, January 2014

011908-13

Li et al.: Adaptive nonlocal means filtering for CT denoising

011908-13

F IG . 11. Adaptive NLM filtering of another clinical A/P scan showing a subcentimeter 0.7 cm renal cell carcinoma metastasis at portal-phase contrast-enhanced abdominal CT. Arrows indicate this small metastasis on full dose (left) and simulated 50% dose (middle), with corresponding CTDIvol of 17.1 and 8.6 mGy, respectively. Right-most image shows adaptive NLM filtering of the 50% dose image, with improved visualization of metastasis and portal and hepatic veins compared to the half-dose image with FBP.

proposed an image-based noise reduction method based on NLM filtering adapted to the local noise level. This adaptive NLM filtering is performed on a compute server with eight GPUs and can service multiple CT systems. The advantage of the proposed method is twofold. First, it can be implemented on any CT scanner, without the need for access to raw CT data, as a software postprocessing step. Second, the adaptation to the local noise level of CT images provides an effec-

tive way to control noise throughout the imaged volume. We have also proposed an efficient and accurate way to estimate the noise distribution in CT images that can be implemented without exact knowledge of the CT system model and reconstruction algorithms. Ideally, the noise distribution of CT images should be calculated based on the noise propagation derived from the reconstruction algorithms implemented on the scanner. Due

F IG . 12. 53-year-old man with history of Crohn’s disease and prior ileocolic resection. Top row is contrast-enhanced CT enterography using routine dose settings at 120 kV (CTDIvol 24.0 mGy) reconstructed with FBP and showing a stricture with mild inflammation at the anastomosis (arrow). Patient subsequently underwent endoscopic dilation. Bottom row is lower dose CT enterography performed two years later with 100 kV (CTDIvol 6.4 mGy) image with adaptive NLM showing recurrent stricture (arrow) with less inflammation, with similar image quality and bowel wall conspicuity despite 73% dose reduction. Medical Physics, Vol. 41, No. 1, January 2014

011908-14

Li et al.: Adaptive nonlocal means filtering for CT denoising

to the complexity of the reconstruction algorithms for helical cone-beam CT, the analytical form of the noise distribution in final images may not be readily available. Therefore, in the current implementation, the noise distribution of CT images was estimated using an analytical method that calculates the noise propagation based on a simple filtered backprojection reconstruction. This can be implemented very efficiently and is independent of the scanner platform. As demonstrated in Fig. 3, this method still yields an estimate of noise distribution within a reasonable range of accuracy using the noise distribution obtained from repeated CT scans as a reference. In our simulation, we only consider the noise caused by photon counting statistics and ignore other noise sources like electronic system noise and scatter radiation noise. This assumption is valid since it is known that noise in CT, typically above the electronic noise floor when appropriate scanning protocols are used, is dominated by photon counting statistics44, 45 and is also validated by the ratios (∼1.414) between reference noise maps from full dose and half dose phantom scans in both sequential and helical mode (data not shown). The bowtie filter also has an important influence on the noise map. Different CT scanners may have different bowtie filters and different scanning tasks (e.g., head, abdominal, etc.) may use different bowtie filters. Characterization of different bowtie filters and use of this information in the denoising may improve both noise estimation accuracy and denoising performance. AEC modulates the x-ray tube current according to the anatomy information from the tomogram. AEC includes both angular and longitudinal axes and thus the tube current is reduced for different projections and scan positions to achieve desired image quality with lower radiation dosage. Hence, AEC has effects on the noise distribution in CT images. In our estimation process, the longitudinal modulation information, which is stored in the DICOM header as the average tube current value of one slice, was utilized to adjust the incident x-ray photon number. Incorporation of tube current at each projection angle should further improve the noise estimate, but this information is not available from the DICOM header. We have performed extensive phantom studies to evaluate the proposed adaptive NLM denoising method. It is well known that the spatial resolution in an IR or a nonlinear denoising method depends on the contrast level of the object. The proposed method can lower noise and achieve higher low-contrast performance without degrading high-contrast spatial resolution. The MTF and SSP studies showed that the high-contrast spatial resolution remained essentially unchanged in both the axial plane and along the longitudinal direction (Figs. 4 and 5), while the noise level is substantially reduced. The NPS results demonstrated that adaptive NLM preserves the shape of the noise power spectrum and the peak frequency better than commercial smoothing kernels for the same level of denoising, which indicates a better maintenance of low-contrast spatial resolution (Fig. 6). The ACR phantom results (Fig. 7) and CHO model observer study demonstrated the potential improvement in low-contrast performance by usMedical Physics, Vol. 41, No. 1, January 2014

011908-14

ing the adaptive NLM filter (Fig. 8). The preliminary clinical pilot within our radiology department showed that adaptive NLM based on an estimation of local noise level can achieve more consistent noise reduction in different regions within the same CT volume than standard NLM (Fig. 2). The clinical patient examples show improved image quality and lower noise while maintaining resolution and preserving small features (Figs. 9–11). Observer performance trials evaluating adaptive NLM use with lower dose CT images are nearing completion for the detection of hepatic metastases and Crohn’s disease.

5. CONCLUSIONS We have proposed a novel approach to practically and efficiently estimate a local noise map from a CT image and shown that it agrees well with the result from multiple scans of the same object on a CT scanner. We have also modified the NLM algorithm to adaptively denoise CT images based on such a local noise level map with only marginal additional computational cost. The adaptive NLM filter significantly improves both inter- and intraslice (and interpatient) denoising performance. A complete evaluation of the effect of this image-space denoising incorporating the noise map method on observer performance is currently under way.

ACKNOWLEDGMENTS This research was partially supported by a grant from Mayo Discovery Translational Fund and a grant from Thrasher Research Fund. J.G.F. and C.H.M. received research support from Siemens Healthcare. L.Y. received royalties from Siemens Healthcare.

APPENDIX: PSEUDOCODE OF ALGORITHM FOR FAST 3D ADAPTIVE NONLOCAL MEANS The pseudocode below makes use of the sum squared difference image technique,23, 24 which calculates the squared difference between corresponding pixels in two images, and then integrates this quantity from one corner of the image (e.g., x = y = 0) so that SSD(X,Y) contains the total sum squared difference for all pixels with (x ≤ X, y ≤ Y). The sum squared difference for arbitrary rectangular patches can then be quickly calculated by simply adding or subtracting the SSD values at the four corners of the patch.23, 24 We have extended this technique to 3D here (which requires eight corners for a cubic patch). To make best use of the sum squared difference image speedup technique, all pixels are filtered simultaneously. The outer loop in the code below cycles through the various possible offsets within a search window, while the inner loop cycles through the pixels in the image volume. Partial results for each pixel are maintained and updated at each cycle through the outer loop.

011908-15

Li et al.: Adaptive nonlocal means filtering for CT denoising 8 A.

Input: v(x,y,z)—input volume; W (Wx,Wy,Wz)—dimensions of search window; P(Px, Py, Pz)—dimensions of patch; k—denoising strength parameter; SD(x,y,z)—noise map; B—blending parameter Output: u(x,y,z)—output volume Temporary variables: SSD(x,y,z)—sum square difference image; SSDpatch—total sum squared difference between two patches; Np —number of pixels in a patch = Px ∗Py ∗Pz ; Z(x,y,z)—sum of weights from all pixels within a search window to the center pixel (excluding itself); M(x,y,z)—maximum weight between all other pixels and the center pixel, empirically used to represent the weight of the noisy center pixel; LUT—lookup table for the exponential weight function (precalculated). Initialize u, M, and Z to be 0 for all δz ∈ [−(Wz−1)/2, (Wz−1)/2] for all δy ∈ [−(Wy−1)/2, (Wy−1)/2] for all δx ∈ [−(Wx−1)/2, (Wx−1)/2] do If (δx, δy, δz) are all zero, the center pixel is acting as its own neighbor, and all SSDs are zero. Skip to the next set of offsets - the center pixel is handled separately below. Compute a shifted version of the image for this offset, and compute the SSD(x,y,z) map for this offset. for all s(x0 ,y0 ,z0 ) in the image v do t = (x0 +δx,y0 +δy,z0 +δz) Calculate SSDpatch between patches centered at s and t. Divide SSDpatch by Np *([k* SD(x0 ,y0 ,z0 )]**2) (the divisions by Np and k can be handled in the LUT precalculation instead for speed). Calculate weight w(s,t) using LUT[SSDpatch] u(s) ← u(s) + w(s,t) × v(t) M(s) = max(M(s),w(s,t)) Z(s) ← Z(s) + w(s,t) end for end for end for end for for all s(x0 ,y0 ,z0 ) in the image do - (this adds the contribution of each pixel to its own filtered value) u(s) ← u(s) + M(s) × v(s) u(s) ← u(s)/(Z(s) + M(s)) end for u = B × u + (1−B) × v return u.

a) Author

to whom correspondence should be addressed. Electronic mail: [email protected]; Telephone: (507)284-8163. 1 D. J. Brenner and E. J. Hall, “Computed tomography: An increasing source of radiation exposure,” New Engl. J. Med. 357, 2277–2284 (2007). 2 J. Hsieh, “Adaptive streak artifact reduction in computed tomography resulting from excessive x-ray photon noise,” Med. Phys. 25, 2139–2147 (1998). 3 M. Kachelriess, O. Watzke, and W. A. Kalender, “Generalized multidimensional adaptive filtering for conventional and spiral single-slice, multi-slice, and cone-beam CT,” Med. Phys. 28, 475–490 (2001). 4 T. Li, X. Li, J. Wang, J. Wen, H. Lu, J. Hsieh, and Z. Liang, “Nonlinear sinogram smoothing for low-dose X-ray CT,” IEEE Trans. Nucl. Sci. 51, 2505–2513 (2004). 5 J. Wang, T. Li, H. Lu, and Z. Liang, “Penalized weighted least-squares approach to sinogram noise reduction and image reconstruction for lowdose X-ray computed tomography,” IEEE Trans. Med. Imaging 25, 1272– 1283 (2006). 6 P. J. La Rivière, “Penalized-likelihood sinogram smoothing for low-dose CT,” Med. Phys. 32, 1676–1683 (2005). 7 P. J. La Rivière, J. Bian, and P. A. Vargas, “Penalized-likelihood sinogram restoration for computed tomography,” IEEE Trans. Med. Imaging 25, 1022–1036 (2006).

Medical Physics, Vol. 41, No. 1, January 2014

011908-15

Manduca, L. Yu, J. D. Trzasko, N. Khaylova, J. M. Kofler, C. M. McCollough, and J. G. Fletcher, “Projection space denoising with bilateral filtering and CT noise modeling for dose reduction in CT,” Med. Phys. 36, 4911–4919 (2009). 9 C. Tomasi and R. Manduchi, “Bilateral filtering for gray and color images,” Int. Conf. Comput. Vis. 6, 839–846 (1998). 10 L. I. Rudin, S. Osher, and E. Fatemi, “Nonlinear total variation based noise removal algorithms,” Physica D 60, 259–268 (1992). 11 A. Buades, B. Coll, and J. M. Morel, “A review of image denoising algorithms, with a new one,” Multiscale Model. Simul. 4, 490–530 (2005). 12 M. Aharon, M. Elad, and A. Bruckstein, “K-SVD: An algorithm for designing overcomplete dictionaries for sparse representation,” IEEE Trans. Signal Process. 54, 4311–4322 (2006). 13 K. Lange and J. A. Fessler, “Globally convergent algorithms for maximum a posteriori transmission tomography,” IEEE Trans. Image Process. 4, 1430–1438 (1995). 14 J. Nuyts, B. De Man, P. Dupont, M. Defrise, P. Suetens, and L. Mortelmans, “Iterative reconstruction for helical CT: A simulation study,” Phys. Med. Biol. 43, 729–737 (1998). 15 J.-B. Thibault, K. D. Sauer, C. A. Bouman, and J. Hsieh, “A threedimensional statistical approach to improved image quality for multislice helical CT,” Med. Phys. 34, 4526–4544 (2007). 16 F. J. Beekman and C. Kamphuis, “Ordered subset reconstruction for x-ray CT,” Phys. Med. Biol. 46, 1835–1844 (2001). 17 J. A. Fessler, E. P. Ficaro, N. H. Clinthorne, and K. Lange, “Groupedcoordinate ascent algorithms for penalized-likelihood transmission image reconstruction,” IEEE Trans. Med. Imaging 16, 166–175 (1997). 18 M. Kachelriess, M. Knaup, and O. Bockenbach, “Hyperfast parallel-beam and cone-beam backprojection using the cell general purpose hardware,” Med. Phys. 34, 1474–1486 (2007). 19 J. S. Kole and F. J. Beekman, “Evaluation of accelerated iterative x-ray CT image reconstruction using floating point graphics hardware,” Phys. Med. Biol. 51, 875–889 (2006). 20 F. Xu and K. Mueller, “Accelerating popular tomographic reconstruction algorithms on commodity PC graphics hardware,” IEEE Trans. Nucl. Sci. 52, 654–663 (2005). 21 F. Xu and K. Mueller, “Real-time 3D computed tomographic reconstruction using commodity graphics hardware,” Phys. Med. Biol. 52, 3405–3419 (2007). 22 D. J. Blezek, Z. Li, B. J. Bartholmai, A. Manduca, and B. J. Erickson, “Clinically feasible CT denoising using a GPU-based non-local means algorithm,” Society for Imaging Informatics in Medicine, 2011 Annual Meeting, Washington, DC, 2011. 23 J. Darbon, A. Cunha, T. F. Chan, S. Osher, and G. J. Jensen, “Fast nonlocal filtering applied to electron cryomicroscopy,” in Proceedings of the IEEE International Symposium on Biomedical Imaging: From Nano to Macro, Paris, France, 2008 (IEEE, New York, NY, 2008), pp. 1331–1334. 24 Y. L. Liu, J. Wang, X. Chen, Y. W. Guo, and Q. S. Peng, “A robust and fast non-local means algorithm for image denoising,” J. Comput. Sci. Technol. 23(2), 270–279 (2008). 25 X. Pan and L. Yu, “Image reconstruction with shift-variant filtration and its implication for noise and resolution properties in fan-beam computed tomography,” Med. Phys. 30, 590–600 (2003). 26 A. Wunderlich and F. Noo, “Image covariance and lesion detectability in direct fan-beam x-ray computed tomography,” Phys. Med. Biol. 53, 2471– 2493 (2008). 27 K. Stierstorfer, A. Rauscher, J. Boese, H. Bruder, S. Schaller, and T. Flohr, “Weighted FBP—a simple approximate 3D FBP algorithm for multislice spiral CT with good dose usage for arbitrary pitch,” Phys. Med. Biol. 49, 2209–2218 (2004). 28 J. H. Hubbell and S. M. Seltzer, “Tables of X-ray mass attenuation coefficients and mass energy absorption coefficients from 1 keV to 20 MeV for elements Z 51 to 92 and 48 additional substances of dosimetric interest,” NISTIR Report No. 5632 (U.S. Department of Commerce, Washington, DC, 1995). 29 H. H. Barrett and K. Myers, Foundations of Image Science, 1st ed. (WileyInterscience, Hoboken, NJ, 2003). 30 H. H. Barrett and W. Swindell, Radiological Imaging: The Theory of Image Formation, Detection, and Processing (Academic, New York, 1981). 31 W. A. Kalender, H. Wolf, and C. Suess, “Dose reduction in CT by anatomically adapted tube current modulation. II. Phantom measurements,” Med. Phys. 26, 2248–2253 (1999).

011908-16

Li et al.: Adaptive nonlocal means filtering for CT denoising

32 J. Hsieh,

Computed Tomography: Principles, Design, Artifacts, and Recent Advances, 2nd Revised ed. (SPIE, Bellingham, WA, 2009). 33 B. R. Whiting, P. Massoumzadeh, O. A. Earl, J. A. O’Sullivan, D. L. Snyder, and J. F. Williamson, “Properties of preprocessed sinogram data in x-ray computed tomography,” Med. Phys. 33, 3290–3303 (2006). 34 M. Gies, W. A. Kalender, H. Wolf, and C. Suess, “Dose reduction in CT by anatomically adapted tube current modulation. I. Simulation studies,” Med. Phys. 26, 2235–2247 (1999). 35 A. C. Kak and M. Slaney, Principles of Computerized Tomographic Imaging (Society of Industrial and Applied Mathematics, Philadelphia, PA, 2001). 36 R. C. Mittelhammer, Mathematical Statistics for Economics and Business, Corrected ed. (Springer, New York, 1996). 37 W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, Numerical Recipes with Source Code CD-ROM 3rd Edition: The Art of Scientific Computing, 3rd ed. (Cambridge University Press, Cambridge, England, 2007). 38 J. M. Boone, “Determination of the presampled MTF in computed tomography,” Med. Phys. 28, 356–360 (2001). 39 J. H. Siewerdsen, I. A. Cunningham, and D. A. Jaffray, “A framework for noise-power spectrum analysis of multidimensional images,” Med. Phys. 29, 2655–2671 (2002).

Medical Physics, Vol. 41, No. 1, January 2014

40 L.

011908-16

Yu, S. Leng, L. Chen, J. M. Kofler, R. E. Carter, and C. H. McCollough, “Prediction of human observer performance in a 2-alternative forced choice low-contrast detection task using channelized Hotelling observer: Impact of radiation dose and reconstruction algorithms,” Med. Phys. 40, 041908 (9pp.) (2013). 41 L. Yu, M. Shiung, D. Jondal, and C. H. McCollough, “Development and validation of a practical lower-dose-simulation tool for optimizing computed tomography scan protocols,” J. Comput. Assist. Tomogr. 36, 477– 487 (2012). 42 G. Yadava, S. Kulkarni, Z. R. Colon, J. Thibault, and J. Hsieh, “TU-A201B-03: Dose reduction and image quality benefits using model based iterative reconstruction (MBIR) technique for computed tomography,” Med. Phys. 37, 3372 (2010). 43 A. Winklehner, C. Karlo, G. Puippe, B. Schmidt, T. Flohr, R. Goetti, T. Pfammatter, T. Frauenfelder, and H. Alkadhi, “Raw data-based iterative reconstruction in body CTA: Evaluation of radiation dose saving potential,” Eur. Radiol. 21, 2521–2526 (2011). 44 Y. Zhang and R. Ning, “Investigation of image noise in cone-beam CT imaging due to photon counting statistics with the Feldkamp algorithm by computer simulations,” J. X-Ray Sci. Technol. 16, 143–158 (2008). 45 P. Massoumzadeh, S. Don, C. F. Hildebolt, K. T. Bae, and B. R. Whiting, “Validation of CT dose-reduction simulation,” Med. Phys. 36, 174–174 (2009).