Abstract

A critical issue in image restoration is the problem of noise removal while keeping the integrity of relevant image information. The method proposed in this paper is a fully automatic 3D blockwise version of the nonlocal (NL) means filter with wavelet subbands mixing. The proposed wavelet subbands mixing is based on a multiresolution approach for improving the quality of image denoising filter. Quantitative validation was carried out on synthetic datasets generated with the BrainWeb simulator. The results show that our NL-means filter with wavelet subbands mixing outperforms the classical implementation of the NL-means filter in terms of denoising quality and computation time. Comparison with wellestablished methods, such as nonlinear diffusion filter and total variation minimization, shows that the proposed NL-means filter produces better denoising results. Finally, qualitative results on real data are presented.

1. Introduction

Image denoising can be considered as a component of processing or as a process itself. In the first case, the image denoising is used to improve the accuracy of various image processing algorithms such as registration or segmentation. Then, the quality of the artifact correction influences performance of the procedure. In the second case, the noise removal aims at improving the image quality for visual inspection. The preservation of relevant image information is important, especially in a medical context.

This paper focuses on a new denoising method firstly introduced by Buades et al. [1] for 2D image denoising: the nonlocal (NL) means filter. We propose, to improve this filter with an automatic tuning of the filtering parameter, a blockwise implementation and a mixing of wavelet su-bands based on the approach proposed in [2]. These contributions lead to a fully-automated method and overcome the main limitation of the classical NL-means: the computational burden.

Section 2 presents related works. Section 3 presents the proposed method with details about our contributions. Section 4 shows the impact of our adaptations compared to different implementations of the NL-means filter and proposes a comparison with well-established methods. The validation experiments are performed on a phantom dataset in a quantitative way. Finally, Section 5 shows results on real data.

Many methods for image denoising have been suggested in the literature, and a complete review of them can be found in [1]. Methods for image restoration aim at preserving the image details and local features while removing the undesirable noise. In many approaches, an initial image is progressively approximated by filtered versions which are smoother or simpler in some sense. Total variation (TV) minimization [3], nonlinear diffusion [4–6], mode filters [7], or regularization methods [3, 8] are among the methods of choice for noise removal. Most of these methods are based on a weighted average of the gray values of the pixels in a spatial neighborhood [9, 10]. One of the earliest examples of such filters has been proposed by Lee [11]. An evolution of this approach has been presented by Tomasi and Manduchi [9] who devised the bilateral filter which includes both a spatial and an intensity neighborhood.

Recently, the relationships between bilateral filtering and local mode filtering [7], local M-estimators [12], and nonlinear diffusion [13] have been established. In the context of statistical methods, the bridge between the Bayesian estimators applied on a Gibbs distribution, resulting with a penalty functional [14] and averaging methods for smoothing, has also been described in [10]. Finally, statistical averaging schemes enhanced via incorporating a variable spatial neighborhood scheme have been proposed in [15–17].

All these methods aim at removing noise while preserving relevant image information. The tradeoff between noise removal and image preservation is performed by tuning the filter parameters, which is not an easy task in practice. In this paper, we propose to overcome this problem with a 3D subbands wavelet mixing. As in [2], we have chosen to combine a multiresolution approach with the NL-means filter [1], which has recently shown very promising results.

Recently introduced by Buades et al. [1], the NL-means filter proposes a new approach for the denoising problem. Contrary to most denoising methods based on a local recovery paradigm, the NL-means filter is based on the idea that any periodic, textured, or natural image has redundancy, and that any voxel of the image has similar voxels that are not necessarily located in a spatial neighborhood. This new nonlocal recovery paradigm allows to improve the two most desired properties of a denoising algorithm: edge preservation and noise removal.

3. Methods

In this section, we introduce the following notations: (i) is the image, where represents the image grid, considered as cubic for the sake of simplicity and without loss of generality ;(ii)for the original voxelwise NL-means approach, (a) is the intensity observed at voxel , (b) is the cubic search volume centered on voxel of size , , (c) is the cubic local neighborhood of of size , , (d) is the vector containing the intensities of (that we term “patch” in the following),(e) is the restored value of voxel , (f) is the weight of voxel when restoring (see Figure 1(a)); (iii)for the blockwise NL-means approach, (a) is the block centered on of size , , (b) is the vector containing the intensities of the block , (c) is the vector containing the restored value of , (d) is the weight of block when restoring the block (see Figure 1(b)),(e) the blocks are centered on voxels which represent a subset of the image voxels, equally regularly distributed over (see Figure 2),(f) represents the distance between the centers of the blocks (see Figure 2).

(a)
(a)
(b)
(b)
Figure 1 
(a) Usual voxelwise NL-means filter: 2D illustration of the NL-means principle. The restored value of voxel (in red) is the weighted average of all intensities of voxels in the search volume , based on the similarity of their intensity neighborhoods and . In this example, we set and . (b) Blockwise NL-means filter: 2D illustration of the blockwise NL-means principle. The restored value of the block is the weighted average of all the blocks in the search volume . In this example, we set and .

Figure 2 
Blockwise NL-means filter. For each block centered on voxel , an NL-means-like restoration is performed from blocks . In this way, for a voxel included in several blocks, several estimations are obtained. The restored value of voxel is the average of the different estimations stored in vector . In this example, , , and .
3.1. The Nonlocal Means Filter

In the classical formulation of the NL means filter [1], the restored intensity of the voxel , is a weighted average of the voxels intensities in the “search volume” of size : where is the weight assigned to value to restore voxel . More precisely, the weight evaluates the similarity between the intensity of the local neighborhoods and centered on voxels and , such that and (cf., Figure 1, Left).

For each voxel in , the computation of the weight is based on the Euclidean distance between patches and , defined as where is a normalization constant ensuring that , and acts as a filtering parameter controlling the decay of the exponential function.

3.1.1. Automatic Tuning of the Filtering Parameter

As explained in the introduction, denoising is usually the first step of complex image processing procedures. The number and the dimensions of the data to process being continually increasing, each step of the procedures needs to be as automatic as possible. In this section, we propose an automatic tuning of the filtering parameter .

First, it has been shown that the optimal smoothing parameter is proportional to the standard deviation of the noise [1]. Second, if we want the filter independent of the neighborhood size, the optimal must depend on (see, (2)). Thus, the automatic tuning of the filtering parameter amounts to determining the relationship , where is a constant.

Firstly, the standard deviation of the noise needs to be estimated. In case of an additive white Gaussian noise, this estimation can be based on pseudoresiduals as defined in [18, 19]. For each voxel of the volume , let us define being the 6-neighborhood at voxel and the constant is used to ensure that in the homogeneous areas. Thus, the standard deviation of noise is computed as Then, in order to make the filter independent of , we used the Euclidean distance normalized by the number of elements: Based on the fact that, in the case of Gaussian noise and with normalized L2-norm, the optimal denoising is obtained for [20], (2) can be written as where only the adjusting constant needs to be manually tuned. If our estimation of the standard deviation of the noise is correct, should be close to 1. The optimal choice for will be discussed later.

3.1.2. Blockwise Implementation

The main problem of the NL-means filter is being its computational time, a blockwise approach can be used to decrease the algorithmic complexity. Indeed, instead of denoising the image at a voxel level, entire blocks are directly restored.

A blockwise implementation of the NL-means filter consists in (a) dividing the volume into blocks with overlapping supports, (b) performing NL-means-like restoration of these blocks, and c) restoring the voxels values based on the restored values of the blocks they belong to, as follows. (1)A partition of the volume into overlapping blocks of size is performed, such as , under the constraint that each block intersects with at least one other block of the partition. These blocks are centered on voxels which constitute a subset of . The voxels are equally distributed at positions , where represents the distance between the centers of . To ensure a global continuity in the denoised image, the overlapping support of blocks is nonempty: .(2)For each block , an NL-means-like restoration is performed as follows: where is a normalization constant ensuring that (see Figure 1, Right).(3) For a voxel included in several blocks , several estimations of the restored intensity are obtained in different . The estimations given by different for a voxel are stored in a vector . The final restored intensity of voxel is then defined as where denotes the element of the vector .

The main advantage of this approach is to significantly reduce the complexity of the algorithm. Indeed, for a volume of size , the global complexity is . For instance, with , the complexity is divided by a factor .

3.1.3. Block Selection

In [21–23], the authors have shown that neglecting the voxels/blocks with small weights (i.e., the most dissimilar patches to the current one) speeds up the filter and significantly improves the denoising results. Indeed, the selection of the most similar patches to the current patch to compute can be viewed as a spatially adaptation of the patch dictionaries. As in [21–23], the preselection of blocks in is based on the mean and the variance of and . The selection tests are given by where and represent, respectively, the mean and the variance of the intensity function for the block centered on the voxel . The new parameters and control the level of rejection related to tests. When and are close to , there is almost no selection and the number of patches taken into account increases: thus the denoised image becomes smoother. The filter is equivalent to the classical NL-means and the computation time increases. When and are close to , the selection is more severe and the number of patches taken into account decreases: the denoised image is less smoothed and the computation time decreases. This kind of selection tends to better enhance the contrast. In practice, and were chosen as in [21, 22]: and .

3.2. Wavelet Subbands Mixing
3.2.1. Hybrid Approaches

Recently, hybrid approaches coupling the NL-means filter and a wavelet decomposition have been proposed [2, 24, 25]. In [24], a wavelet-based denoising of blocks is performed before the computation of the nonlocal means. The NL-means filter is performed with denoised version of blocks in order to improve the denoising result. In [25], the NL-means filter is applied directly on wavelet coefficients in transform domain. This approach allows a direct denoising of compressed images (such as JPEG2000) and a reduction of computational time since smaller images are processed. In [2], a multiresolution framework is proposed to adaptively combine the result of denoising algorithms at different space-frequency resolutions. This idea relies on the fact that a set of filtering parameters is not optimal over all the space-frequency resolutions. Thus, by combining to the transform domain the results obtained with different sets of filtering parameters, the denoising is expected to be improved.

3.2.2. Overall Processing

In order to improve the denoising result of the NL-means filter, we propose a multiresolution framework similar to [2] to implicitly adapt the filtering parameters (, ) over the different space-frequency resolutions of the image. This adaptation is based on the fact that the size of the patches impacts the denoising properties of the NL-means filter. Indeed, the weight given to a block depends on its similarity with the block under consideration, but the similarity between the blocks depends on their sizes. Thus, given the size of the blocks, removal or preservation of image components can be favored.

In the transform domain, the main features of the image correspond to low-frequency information while finer details and noise are associated to high frequencies. Nonetheless, noise is not a pure high-frequency component in most images. Noise is spanned over a certain range of frequencies in the image with mainly middle and high components [2].

In NL-means-based restoration, large blocks and setting efficiently remove all frequencies of noise but tend to spoil the main features of the image, whereas small blocks and low smoothing parameter tend to better preserve the image components but cannot completely remove all frequencies of noise. As a consequence, we propose the following workflow (see Figure 3). (i)Denoising of the original image using two sets of filtering parameters: one adapted to the noise components removal (i.e., large blocks and ) and the other adapted to the image features preservation (i.e., small blocks and ). This yields two images and . In , the noise is efficiently removed and, conversely, in , the image features are preserved.(ii)Decomposing and into low- and high-frequency subbands. The first level decomposition of the images is performed with a 3D discrete wavelet transform (DWT).(iii)Mixing the highest-frequency subbands of and the lowest frequency subbands of .(iv)Reconstructing the final image by an inverse 3D DWT from the combination of the selected high and low frequencies.


Figure 3 
Workflow. First, the noisy image I is denoised with two sets of filtering parameters and . Then, and are decomposed into low- and high-frequency subbands by 3D DWT. The four lowest frequency subbands of (i.e., , , , and ) are mixed with the four highest-frequency subbands of (i.e., , , , and ). Finally, the result image is obtained by inverse 3D DWT of the selected subbands.

In this paper, we propose an implementation of this approach using our optimized blockwise NL-means filter and the 3D DWT Daubechies-8 basis. The latter is implemented in Qccpack (http://qccpack.sourceforge.net) in the form of dyadic subband pyramids. This DWT is widely used in image compression due to its robustness and efficiency.

3.2.3. Selection of Wavelet Subbands

Once the original image has been denoised using two sets of filtering parameters, a 3D DWT at the first level is performed on both and images. For each image, eight subbands are obtained: , , , , , , , and . (i)In the eight wavelet subbands obtained with , the frequencies corresponding to noise are efficiently removed from the high frequencies whereas the low frequencies associated to the main features are spoiled.(ii)In the eight wavelet subbands obtained with , the low frequencies associated to main features are efficiently preserved whereas residual frequencies corresponding to noise are present in high frequencies.

Thus, we select the highest frequencies of (i.e., , , , and ) and the lowest frequencies of (i.e., , , , and ). Then, the 4 lowest subbands of are combined with the 4 highest subbands of . Finally, an inverse 3D DWT is performed on these 8 selected subbands to obtain the final denoised image (see Figure 3).

In [21, 22], the optimal parameters for 3D MRI have been estimated as , , , and . In our experiments, the two sets of parameters used to obtain and were and . Compared to [21, 22], the size of “search volume” was reduced to decrease the computational time. Several sets of parameters have been tested, the mentioned numerical values are satisfying to balance the denoising performance (high PSNR values) and computational burden. Finally, to decrease the computational time, this workflow is parallelized and each version is computed on different CPUs or cores (see Figure 3).

4. Validation on a Phantom Data Set

4.1. Materials

In order to evaluate the performance of the different variants of the NL-means filter on 3D MR images, tests were performed on the BrainWeb database [26]. Several images were simulated to validate the performance of the denoising on various images: (a) T1-w phantom MRI for 4 levels of noise 3%, 9%, 15%, and 21% and (b) T2-w phantom MRI with multiple sclerosis (MS) lesions for 4 levels of noise 3%, 9%, 15%, and 21%. A white Gaussian noise was added, and the notations of BrainWeb are used: a noise of is equivalent to , where is the value of the highest voxel intensity of the phantom (150 for T1-w and 250 for T2-w).

4.2. Comparison with Different NL-Means Filters

In the following, let us define the following. (i)NL-means: standard voxelwise implementation with automatic tuning of the filtering parameter () [1].(ii)Optimized NL-means: voxelwise implementation with automatic tuning of the filtering parameter () and voxels selection presented in [21].(iii)Optimized blockwise NL-means: (This filter can be freely tested at http://www.irisa.fr/visages/bench marks) blockwise implementation with automatic tuning of the filtering parameter () and blocks selection presented in [22].(iv)Optimized blockwise NL-means with wavelet mixing: proposed filter based on a blockwise implementation, an automatic tuning of the filtering parameter (), a block selection, and a wavelet subbands mixing. The selected filtering parameters for the different implementations were as follows. (i)For the NL-means and optimized NL-means filters, the parameters are those used in [21]: ,,,, and .(ii)Concerning the optimized blockwise NL-means filter, the sets of parameters are those used in [22]: ,,,,, and .(iii)Finally, for the optimized blockwise NL-means with wavelet mixing filter the parameter are the following: , , ,, and . For 8-bit encoded images, the PSNR is defined as follows: where RMSE denotes the root mean square error estimated between the ground truth and the denoised image. For the sake of clarity, the PSNR values are estimated only in the region of interest (cerebral tissues) obtained by removing the background (i.e., the label 0 of the discrete model in Brainweb).

Firstly, we have experimentally verified that the optimal denoising is obtained for for high levels of noise and for low levels of noise. These results account for the error in the estimation of ( at , at , at , and at ) (see Figure 4). The parameter was fixed to 1 for all the experiments.


Figure 4 
Influence of the filtering parameter on the PSNR according to and for several levels of noise. These results are obtained with the optimized blockwise NL-means filter on the T1-w phantom MRI and account for the error in the estimation of .
4.2.1. Quantitative Results

Table 1 shows that the blockwise approach of the NL-means filter, with and without voxels selection (see (9)), allows to drastically reduce the computational time. With a distance between the block centers corresponding to , the blockwise approach divides the timings by a factor superior to (see Table 1). However, the computational time reduction is balanced with a slight decrease of the PSNR (see Figure 5) compared to the optimized NL-means filter presented in [21]. Our optimized blockwise NL-means with wavelet mixing allows to compensate this slight decrease of the PSNR and to divide the computational by a factor compared to the optimized NL-means filter.

Table 1 
Comparison of different implementations of NL-means in terms of computational time and denoising quality. The computational time was obtained with multithreading on a DualCore Intel(R) Pentium(R) D CPU 3.40 GHz. These results were obtained on a T1-w phantom image of voxels with of noise.
(a)
(a)
(b)
(b)
Figure 5 
Comparison of the different NL-means filters on T1-w phantom MRI and T2-w phantom MRI with MS.
4.2.2. Visual Assessment

Visually, the proposed method combines the most important attributes of a denoising algorithm: edge preservation and noise removal. Figure 6 shows that our filter removes noise while keeping the integrity of MS lesions (i.e., no structure appears in the removed noise). Figure 7 focuses on the differences between the optimized blockwise NLM and the optimized blockwise NLM with WM filters. The denoising result obtained with the optimized blockwise NLM with WM filter visually preserves the edges better than the optimized blockwise NLM filter. This is also confirmed by visual inspection of the comparison with the “ground truth”. The images of difference between the phantom and the denoised image (see bottom of Figure 7) show that less structures have been removed with the optimized blockwise NLM with WM filter. Thus, the multiresolution approach allows to better preserve the edges and to enhance the contrast between tissues.

(a)
(a)
(b)
(b)
(c)
(c)
(d)
(d)
Figure 6 
Fully automatic restoration obtained with the optimized blockwise NL-means with wavelet mixing filter in 3 minutes on a DualCore Intel(R) Pentium(R) D CPU 3.40 GHz. The image is a T2-w phantom MRI with MS of voxels and 9% of noise.

Figure 7 
Top. Phantom and phantom noisy with 9%. Middle. The denoising result obtained with the optimized blockwise NLM with WM filter and the optimized blockwise NLM filter. Bottom. The image of difference between the phantom and the denoising result (i.e., - ). The contrast of the zooms have been artificially increased. Visually, less structures have been removed with the optimized blockwise NLM with WM filter.

4.3. Comparison with Other Methods

In this section, we compare the proposed method with two of the most used approaches in MRI domain: the nonlinear diffusion (NLD) filter r [4] and the total variation (TV) minimization [3]. The main difficulty to achieve this comparison is related to the tuning of smoothing parameters in order to obtain the best results for NLD filter and TV minimization scheme. After quantifying the parameter space, we exhaustively tested all possible parameters within a certain range. This allows us to obtain the best possible results for the NLD filter and the TV minimization.

For the optimized blockwise NLM with WM, the same set of parameters and are used for all noise levels. The automatic tuning of adapts the smoothing to the noise level.

For NLD filter, the parameter K varied from to with a step of and the number of iterations varied from to . For TV minimization, the parameter varied from 0.01 to 1 with a step of 0.01 and the number of iterations varied from 1 to 10. The results obtained for a of Gaussian noise are presented in Figure 8, but this screening was performed for the four levels of noise. It is important to underline that the results giving the best PSNR are used, but these results do not necessarily give the best visual output. Actually, the best PSNR value for the NLD filter and TV minimization are obtained for a visually under-smoothed image since these methods tend to spoil the edges. This is explained by the fact that the optimal PSNR is obtained when a good tradeoff is reached between edge preserving and noise removing.

(a)
(a)
(b)
(b)
Figure 8 
Result for the NLD filter and the TV minimization on phantom images with Gaussian noise at . For the NLD filter, K varied from to with a step of and the number of iterations varied from to . For the TV minimization, varied from 0.01 to 1 with a step of 0.01 and the number of iterations varied from 1 to 10.

4.3.1. Quantitative Results

As presented in Figure 9, our block-optimized NL-means with wavelet mixing filter produced the best PSNR values whatever the noise level is. On average, a gain of  dB is achieved compared to TV minimization and AD filter. The PSNR value between the noisy image and the ground truth is called “No processing” and is used as a reference.


Figure 9 
Comparison between nonlinear diffusion, total variation, and optimized blockwise NL-means with wavelet mixing denoising. The PSNR experiments show that the optimized blockwise NL-means with wavelet mixing filter significantly outperforms the well-established total variation minimization process and the nonlinear diffusion approach.

4.3.2. Visual Assessment

Figure 10 shows the denoising results obtained by the NLD filter, the TV minimization, and our optimized blockwise NLM with WM. Visually, the NL-means-based approach produced the best denoising. The removed noise (see middle of Figure 10) shows that the proposed method removes significantly less structures than NLD filter or TV minimization. Finally, the comparison with the “ground truth” underlines that the NL-means restoration gives a result very close to the “ground truth” and better preserves the anatomical structure compared to NLD filter and TV minimization.


Figure 10 
Comparison between nonlinear diffusion, total variation, and our optimized blockwise NL-means with wavelet mixing denoising on synthetic T1-w images. Top. Zooms on T1-w BrainWeb images. Left. The “ground truth”. Right. The noisy images with of Gaussian noise. Middle. The results of restoration obtained with the different methods and the images of the removed noise (i.e., the difference (centered on 128) between the noisy image and the denoised image. Bottom. The difference (centered on 128) between the denoised image and the ground truth. Left. Nonlinear diffusion denoising. Left. Nonlinear diffusion denoising. Middle. Total variation minimization process. Right. Optimized Blockwise NL-means with WM filter. The NL-means-based restoration better preserves the anatomical structure in the image while efficiently removing the noise as it can be seen in the image of removed noise.

5. Experiments on Clinical Data

The T1-weighted MR images used for experiments were obtained with T1 sense 3D sequence on 3T Philips Gyroscan scanner. The restoration results, presented in Figure 11, show good preservation of the cerebellum. Fully automatic segmentation and quantitative analysis of such structures are still a challenge, and improving restoration schemes could greatly improve these processings.

(a)
(a)
(b)
(b)
(c)
(c)
Figure 11 
Fully automatic restoration obtained with the optimized blockwise NL-means with wavelet mixing filter on a 3 Tesla T1-w MRI data of voxels in less than 4 minutes on a DualCore Intel(R) Pentium(R) D CPU 3.40 GHz.

6. Discussion and Conclusion

This paper presented a fully automated blockwise version of the nonlocal means filter with subbands wavelet mixing. Experiments were carried out on the BrainWeb dataset [26] and real dataset. The results on phantom shows that the proposed optimized blockwise NL-means with subbands wavelet mixing filter outperforms the classical implementation of the NL-means filter and the optimized implementation presented in [21, 22], in terms of PSNR values and computational time. Compared to the classical NL-means filter, our implementation (with block selection, blockwise implementation, and wavelet subbands mixing) considerably decreases the required computational time (up to a factor of 20) and significantly increases the PSNR of the denoised image. The comparison of the filtering process with and without wavelet mixing shows that the subbands mixing better preserves edges and better enhances the contrast between the tissues. This multiresolution approach allows to adapt the smoothing parameters along the frequencies by combining several denoised images. The comparison with well-established methods such as NLD filter and TV minimization shows that the NL-means-based restoration produces better results. Finally, the impact of the proposed multiresolution approach based on wavelet subbands mixing should be investigated further, for instance, when combined to the nonlinear diffusion filter [4] and the total variation minimization [3].