Efficient and Robust Non-Local Means Denoising Methods for Biomedical Images

Denoising is an important step to improve image quality and to increase the performance of image analysis. However, conventional partial differential equation based image denoising methods, especially total variation functional minimization techniques, do not work very well on biomedical images such as magnetic resonance images (MRI), ultrasound, and X-ray images. These images present small structures with signals barely detectable above the noise level which involve more complex noise and unclear edges. The recently developed non-local means (NLM) filtering method can treat these types of images better. The standard NLM filter uses the weighted averages of similar patches present in the image. Since the NLM filter is a non-local averaging method, it is very accurate in removing noise but has computational complexity. We develop efficient and optimized NLM based methods and their associate numerical algorithms. The new methods are still accurate enough and more efficient than the original NLM filter. Numerical results show that the new methods compare favorably to the conventional denoising methods.


Introduction
Image denoising is among the most fundamental problems in image processing with many applications. A large range of methods covering various fields of mathematics and engineering are available for denoising an image. The initial denoising models were derived from energy minimization using nonlinear partial differential equations (PDEs) [1]- [6]. The filtering based models using smoothing operators have also been used for quite a long time especially for engineering applications [7]- [11]. The most successful one among them was the non-local means method proposed by Buades, et al. in 2005 [11]. There have also been some efforts to extend some known PDEs and variational techniques in image processing to the non-local framework [12]- [14].
The filtering based models are mainly based on convolutions. The image is convoluted with a symmetric kernel in effort to smooth out the noise. To remove noise, the affected pixels are smoothed out using information from surrounding pixels. The choice of neighboring pixels play an important role in these models. The most basic smoothing filter is the Gaussian filter, where the image is locally convoluted with Gaussian kernel. It has been shown by Buades,et al. in [11] that the Gaussian convolution is optimal for harmonic functions, but it destroys edges and textures. The first non-local filtering method, where neighbors are chosen from the whole image as opposed to the adjacent pixels in local filters, was developed by Yaroslavsky [7]. The neighborhood filter has a big advantage of not blurring the edges much, compared to the local filters. However, the method will not give expected results when both pixel values are noisy. Moreover, it was also shown in [11] that neighborhood filters can create artificial shocks. In [11], the authors proposed the new state of art image denoising algorithm known as non-local means (NLM) algorithm. The method is non-local since each pixel is denoised using the weighted average of all the pixels in the entire image. Several numerical examples in [11] showed that the NLM method is very accurate in removing noise in the given image compared to other conventional filtering methods and even to some PDE based methods. Also, the method works very well on biomedical images which have complex noise and unclear edges. However, the main drawback of the non-local means denoising method is on computational complexity. Due to the non-locality, the method is very slow in achieving the results and hence quite impractical. Coupé et al. [15] reduced the computational time by developing an optimized NLM filter using a blockwise implementation and a pre-selection of the most relevant blocks.
In this paper, we develop new blockwise NLM filters based on the model developed by Coupé et al. [15]. Section 2 comprises of preliminary with more detailed information on the standard NLM model and the blockwise NLM method developed by Coupé et al. In Section 3, new methods and their numerical procedures are introduced and discussed. The numerical results in Section 4 show that the new methods outperform the conventional ones. Finally, we summarize our outcomes in Section 5.

Preliminary
In this section, we describe the two filtering models, the original NLM model [11] and Coupé et al.'s blockwise NLM [15].

Non-local means filtering method
The NLM filter [11] restores every pixel in the image by computing the weighted average of all the pixels in the image. The weight calculates the similarity between the intensity of the local neighborhoods. The method is given by the formula: where the weight w(i, j) depends on the similarity of the i and j pixels and satisfies the conditions 0 ≤ w(i, j) ≤ 1 and j w(i, j) = 1. For estimating the weight between the two pixels, the authors considered the similarities of a square neighborhood N k of a fixed size at each pixel k. The similarity is measured as a decreasing function of the Gaussian-weighted Euclidean distance [11], ||u(N i ) − u(N j )|| 2 2,a . The weight function is defined as where C(i) is the normalizing constant given by The similar neighborhoods have a larger weight which contributes more to the approximation NL[u](i). Therefore, the method produces more accurate denoising results. For computational purposes, the authors reduced the search window of each pixel to a window of size 21 × 21 pixels and the similarity square neighborhood N i to a 7 × 7 block. Therefore for the image of size N × N, the final complexity of the algorithm is in the order of O(49 · 441 · N 2 ).

Blockwise non-local means method
The standard NLM method is very accurate in removing noise compared to other conventional denoising methods. However, the main drawback is on computational inefficiency due to its non-locality. To overcome the drawback, the authors of [15] suggested the use of blocks B i k of size (2α + 1) 2 with overlapping supports and performed NLM restoration on these blocks instead of pixels. The following is the proposed formula to restore images: where Z i k is a normalization constant satisfying B j ∈V i k w(B i k , B j ) = 1 andĥ is a smoothing parameter. Inspired by Mahmoudi and Sapiro [16] who proposed a method to pre-select the most similar neighborhoods only in order to avoid unnecessary weight computation, the weight w(B i k , B j ) was only considered when the ratios of mean and variance of u(B i k ) and u(B j ) are within sufficiently small ranges around 1. This blockwise NLM method significantly reduced the computation time for the NLM method.

New blockwise based non-local means algorithms
In this section, we develop more efficient denoising algorithms based on the blockwise NLM method developed by Coupé et al. [15]. The blockwise NLM method significantly reduced the computational time of the standard NLM method. The blockwise NLM filter is more practical in medical applications since the MRI and X-ray images are mostly 3D or stacks of 2D, which requires a large amount of computation to achieve desirable denoising results. Therefore, computational efficiency plays a crucial role in biomedical image denoising as long as it does not lose much accuracy. In [15], authors considered blocks overlapping with their neighboring blocks by one column or one row of pixels to ensure global continuity. As a result, they obtained multiple approximations of the restored intensity NL(u)(x i ) at the pixel x i for up to 4 values. The final restored intensity value at the pixel x i was then averaged We consider blocks that are still overlapping with neighboring blocks and perform a NLM-like restoration Eqs. (5)- (6). But instead of calculating the average of multiple intensity approximations, we replace the value with the most current value as we proceed to the next block. We also consider non-overlapping blocks to have just one intensity approximation at each pixel. In both cases, the denoised images are globally continuous although the new approaches may produce more checkerboard effects. Both methods do not require extra computations of averaging multiple intensities at overlapping blocks and therefore enhance the computational efficiency of the original blockwise NLM developed by Coupé et al. To further reduce the computational complexity, we also adopt the strategy in [15] that the Gaussian-weighted Euclidean distance || · || 2 2,a is replaced by the classical Euclidean distance || · || 2 2 . We summarize the steps of new algorithms as follows: 1. Use a standard Euclidean norm for calculating the weights. 3. Weights are calculated only when mean and variance of intensities of two blocks are similar: Var(u(B j )) < 1 for some appropriate parameters µ 1 and σ 1 . Here, u(B i k ) and Var(u(B i k )) represent the mean and the variance of u(B i k ), respectively.

Numerical results
In this section, we test our new algorithms, OBR-NLM and NOB-NLM, using two different images: Lenna and an MRI image for a human brain (Brain-MRI). The results are compared to standard NLM and Coupé et al.'s Blockwise (CB)-NLM model. The CB-NLM adopted other techniques such as Automatic Tuning of NLM filter and Parallel Computation [15], to further improve the computational efficiency of the standard NLM. For fair comparison with our new methods, we only apply the selective weight computation Eq. (7) to all blockwise NLM methods. All algorithms presented in this section are implemented in Matlab using an Intel Core i7 CPU 1.60 GHz Desktop Processor.
The peak-signal-to-noise ratio (PSNR) is defined as follows to measure the accuracy of the methods: PSNR = 10 log 10 Here, g is a noise free clean image and u is the restored image obtained by denoising a noisy version of g. For a pre-selection of weights as described in Eq. (7), the optimal values of  the parameters µ 1 and σ 2 1 were chosen experimentally in [15] and they were 0.95 and 0.5, respectively. We follow their suggestions but slightly modify the value of µ 1 to be 0.8 to obtain the best results. For Lenna, methods were tested on the image where 15% of Gaussian noise was synthetically added to the clean image. For Brain-MRI, we added 25% of Gaussian noise, more noise than Lenna, to mimic the case when medical images have considerable amount of noise and artifacts. For blockwise algorithms, the smoothing parameterĥ was set to 100 for Lenna and 500 for Brain-MRI. The smoothing parameter h for the standard NLM was chosen as 12 for Lenna and 15 for Brain-MRI.
As one can see in Table 1, the new algorithms OBR-NLM and NOB-NLM both outperform the standard NLM and CB-NLM in terms of computational efficiency. The visual verification in Figure 1 and Figure 2, and the PSNR values in Table 1 show that the new methods still do not lose much accuracy. One can notice that the PSNR values for Brain-MRI in Table 1 were greatly improved although the image was quite seriously contaminated by 25% of Gaussian  noise. It can also be seen in Table 1

Conclusions
In this article, the authors have studied efficient and robust filtering based methods to deal with biomedical images, which have complex noise and unclear edges. The computational efficiency of the methods for denoising biomedical images needs to be optimal to be applied to a larger size images, such as 3D scan of MRI or 2D stacks of X-ray images.