3D Wavelet Subbands Mixing for Image Denoising

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.


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 pro-poses 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.

RELATED WORKS
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][5][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][16][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.

METHODS
In this section, we introduce the following notations: where Ω 3 represents the image grid, considered as cubic for the sake of simplicity and without loss of generality (|Ω 3 | = N 3 ); (ii) for the original voxelwise NL-means approach, (1) (N i ), . . . , u (|Ni|) (N i )) T is the vector containing the intensities of N i (that we term "patch" in the following), (e) NL(u)(x i ) is the restored value of voxel x i , (f) w(x i , x j ) is the weight of voxel x j when restoring u(x i ) (see Figure 1(a));    Figure 1(b)), (e) the blocks B ik are centered on voxels x ik which represent a subset of the image voxels, equally regularly distributed over Ω 3 (see Figure 2), (f) n represents the distance between the centers of the blocks B ik (see Figure 2).

The nonlocal means filter
In the classical formulation of the NL means filter [1],the restored intensity NL(u)(x i ) of the voxel x i , is a weighted average of the voxels intensities u(x i ) in the "search volume" V i of size (2M + 1) 3 : where w(x i , x j ) is the weight assigned to value u(x j ) to restore voxel x i . More precisely, the weight evaluates the similarity between the intensity of the local neighborhoods N i and N j centered on voxels x i and x j , such that w(x i , x j ) ∈ [0, 1] and Figure 1, Left).  For each voxel x j in V i , the computation of the weight is based on the Euclidean distance between patches u(N j ) and u(N i ), defined as where Z i is a normalization constant ensuring that j w(x i , x j ) = 1, and h acts as a filtering parameter controlling the decay of the exponential function.

Automatic tuning of the filtering parameter h
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 auto-matic as possible. In this section, we propose an automatic tuning of the filtering parameter h.
First, it has been shown that the optimal smoothing parameter h is proportional to the standard deviation of the noise σ [1]. Second, if we want the filter independent of the neighborhood size, the optimal h must depend on |N i | (see, (2)). Thus, the automatic tuning of the filtering parameter h amounts to determining the relationship 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 i as defined in [18,19]. For each voxel x i of the volume Ω 3 , let us define Pierrick Coupé et al. P i being the 6-neighborhood at voxel x i and the constant √ 6/7 is used to ensure that E[ 2 i ] = σ 2 in the homogeneous areas. Thus, the standard deviation of noise σ is computed as Then, in order to make the filter independent of |N i |, we used the Euclidean distance · 2 2 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 h 2 = 2σ 2 [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.

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 Ω 3 into overlapping blocks under the constraint that each block B ik intersects with at least one other block of the partition. These blocks are centered on voxels x ik which constitute a subset of Ω 3 . The voxels x ik are equally distributed at positions where n represents the distance between the centers of B ik . To ensure a global continuity in the denoised image, the overlapping support of blocks is nonempty: 2α ≥ n.
(2) For each block B ik , an NL-means-like restoration is performed as follows: where Z ik is a normalization constant ensuring that j w(B ik , B j ) = 1 (see Figure 1, Right).
where A i (p) denotes the pth element of the vector A i .
The main advantage of this approach is to significantly reduce the complexity of the algorithm. Indeed, for a volume Ω 3 of size N 3 , the global complexity is O((2α + 1) 3 (2M + 1) 3 ((N − n)/n) 3 ). For instance, with n = 2, the complexity is divided by a factor 8.

Block selection
In [21][22][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 u(B j ) to the current patch u(B i ) to compute NL(u)(B i ) can be viewed as a spatially adaptation of the patch dictionaries. As in [21][22][23], the preselection of blocks in V i is based on the mean and the variance of u(B i ) and u(B j ). The selection tests are given by where u(B ik ) and Var(u(B ik )) represent, respectively, the mean and the variance of the intensity function for the block B ik centered on the voxel x ik . The new parameters 0 < μ 1 < 1 and 0 < σ 1 < 1 control the level of rejection related to tests. When μ 1 and σ 1 are close to 0, 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 μ 1 and σ 1 are close to 1, the selection is more severe and the number of patches taken into account decreases: the denoised image is less smoothed and the computation time Pierrick Coupé et al. decreases. This kind of selection tends to better enhance the contrast. In practice, μ 1 and σ 1 were chosen as in [21,22]: μ 1 = 0.95 and σ 1 = 0.5.

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.

Overall processing
In order to improve the denoising result of the NL-means filter, we propose a multiresolution framework similar to [2]   ferent 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 β = 1 efficiently remove all frequencies of noise but tend to spoil the main features of the image, whereas small blocks and low smoothing parameter (β = 0.5) 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). 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.

Selection of wavelet subbands
Once the original image I has been denoised using two sets of filtering parameters, a 3D DWT at the first level is performed on both I o and I u images. For each image, eight subbands are obtained: LLL 1 , LLH 1 , LHL 1 , HLL 1 , LHH 1 , HLH 1 , HHL 1 , and HHH 1 .
(i) In the eight wavelet subbands obtained with I o , 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 I u , 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 o (i.e., LHH 1 , HLH 1 , HHL 1 , and HHH 1 ) and the lowest frequencies of I u (i.e., LLL 1 , LLH 1 , LHL 1 , and HLL 1 ). Then, the 4 lowest subbands of I u are combined with the 4 highest subbands of I o . 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 α = 1, M = 5, μ 1 = 0.95, and σ 1 = 0.5. In our experiments, the two sets of parameters used to obtain I u and I o were S u = (α u , M W , β u ) = (1, 3, 0, 5) and (2, 3, 1). 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).

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 3% is equivalent to N (0, ν(3/100)), where ν is the value of the highest voxel intensity of the phantom (150 for T1-w and 250 for T2-w).

Comparison with different NL-means filters
In the following, let us define the following.
(ii) Optimized NL-means: voxelwise implementation with automatic tuning of the filtering parameter h (β = 1) 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 h (β = 1) 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 h (β = 1), 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]: 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 β ≈ 1 for high levels of noise and β ≈ 0.5 for low levels of noise. These results account for the error in the estimation of σ ( σ 2 = 3.42% at 3%, σ 2 = 7.93% at 9%, σ 2 = 12.72% at 15%, and σ 2 = 17.44% at 21%) (see Figure 4). The parameter β was fixed to 1 for all the experiments. 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 n = 2, the blockwise approach divides the timings by a factor superior to 5 (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 4 compared to the optimized NL-means filter.

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.

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 S u = (α u , M W , β u ) = (1, 3, 0.5) and S o = (α o , M W , β o ) = (2, 3, 1) are used for all noise levels. The automatic tuning of h adapts the smoothing to the noise level.
For NLD filter, the parameter K varied from 0.05 to 1 with a step of 0.05 and the number of iterations varied from 1 to 10. 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 9% 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.

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 2.15 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.
International Journal of Biomedical Imaging 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.

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.

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 wellestablished 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].