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 this paper focuses on a new denoising method firstly introduced by buades et al .   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 . these contributions lead to a fully - automated method and overcome the main limitation of the classical nl - means : the computational burden . 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 . many methods for image denoising have been suggested in the literature , and a complete review of them can be found in . 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 , nonlinear diffusion 
 [ 46 ] , mode filters , 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 ] . an evolution of this approach has been presented by tomasi and manduchi   who devised the bilateral filter which includes both a spatial and an intensity neighborhood . recently , the relationships between bilateral filtering and local mode filtering , local m - estimators , and nonlinear diffusion   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   and averaging methods for smoothing , has also been described in . finally , statistical averaging schemes enhanced via incorporating a variable spatial neighborhood scheme have been proposed in [ 1517 ] . 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 , we have chosen to combine a multiresolution approach with the nl - means filter , which has recently shown very promising results . recently introduced by buades et al . , 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 . in this section , we introduce the following notations : 

 u :    is the image , 
 where  represents the image grid , considered as cubic for the sake of simplicity and without loss of generality ( || = n);for the original voxelwise nl - means approach , 
 u(xi ) is the intensity observed at voxel xi , vi is the cubic search volume centered on voxel xi of size |vi| = ( 2 m + 1 ) , m   , ni is the cubic local neighborhood of xi of size |ni| = ( 2d + 1 ) , d   , u(ni ) = ( u ( ni ) ,  , u(|ni 
 |)(ni ) ) is the vector containing the intensities of ni 
  ( that we term  patch  in the following),nl(u)(xi ) is the restored value of voxel 
 xi , w(xi , xj ) 
  is the weight of voxel xj when restoring u(xi ) ( see figure 1(a ) ) ; 
 for the blockwise nl - means approach , 
 bi is the block 
 centered on xi of size |bi| = ( 2 + 1 ) ,    , u(bi ) is the vector 
 containing the intensities of the block bi , nl(u)(bi ) is the vector 
 containing the restored value of bi , w(bi , bj ) is the weight 
 of block bj   when restoring 
 the block u(bi ) ( see figure 1(b ) ) , the 
 blocks bik are centered on 
 voxels xik 
  which represent 
 a subset of the image voxels , equally regularly distributed over  ( see figure 2),n   represents the 
 distance between the centers of the blocks bik 
  ( see figure 2 ) . u :    is the image , 
 where  represents the image grid , considered as cubic for the sake of simplicity and without loss of generality ( || = n ) ; for the original voxelwise nl - means approach , 
 u(xi ) is the intensity observed at voxel xi , vi is the cubic search volume centered on voxel xi of size |vi| = ( 2 m + 1 ) , m   , ni is the cubic local neighborhood of xi of size |ni| = ( 2d + 1 ) , d   , u(ni ) = ( u ( ni ) ,  , u(|ni 
 |)(ni ) ) is the vector containing the intensities of ni 
  ( that we term  patch  in the following),nl(u)(xi ) is the restored value of voxel 
 xi , w(xi , xj ) 
  is the weight of voxel xj when restoring u(xi ) ( see figure 1(a ) ) ; 
  u(xi ) is the intensity observed at voxel xi ,   vi is the cubic search volume centered on voxel xi of size |vi| = ( 2 m + 1 ) , m   ,   ni is the cubic local neighborhood of xi of size |ni| = ( 2d + 1 ) , d   , u(ni ) = ( u ( ni ) ,  , u(|ni 
 |)(ni ) ) is the vector containing the intensities of ni 
  ( that we term  patch  in the following ) , nl(u)(xi ) is the restored value of voxel 
 xi ,   w(xi , xj ) 
  is the weight of voxel xj when restoring u(xi ) ( see figure 1(a ) ) ; for the blockwise nl - means approach , 
 bi is the block 
 centered on xi of size |bi| = ( 2 + 1 ) ,    , u(bi ) is the vector 
 containing the intensities of the block bi , nl(u)(bi ) is the vector 
 containing the restored value of bi , w(bi , bj ) is the weight 
 of block bj   when restoring 
 the block u(bi ) ( see figure 1(b ) ) , the 
 blocks bik are centered on 
 voxels xik 
  which represent 
 a subset of the image voxels , equally regularly distributed over  ( see figure 2),n   represents the 
 distance between the centers of the blocks bik 
  ( see figure 2 ) . bi is the block 
 centered on xi of size |bi| = ( 2 + 1 ) ,    ,   u(bi ) is the vector 
 containing the intensities of the block bi ,   nl(u)(bi ) is the vector 
 containing the restored value of bi ,   w(bi , bj ) is the weight 
 of block bj   when restoring 
 the block u(bi ) ( see figure 1(b ) ) ,   the 
 blocks bik are centered on 
 voxels xik 
  which represent 
 a subset of the image voxels , equally regularly distributed over  ( see figure 2 ) , n   represents the 
 distance between the centers of the blocks bik 
  ( see figure 2 ) . in the classical formulation of the nl means filter , the restored 
 intensity nl(u ) ( xi ) of the voxel xi , is a weighted average of the voxels intensities u(xi ) in the  search volume  vi of size ( 2 m + 1 ) : 
			
 ( 1)nl(u)(xi)=xjviw(xi , xj)u(xj ) , 
 where w(xi , xj ) is the weight 
 assigned to value u(xj ) to restore 
 voxel xi . more precisely , the weight evaluates the similarity between the intensity of the local neighborhoods ni and nj centered on voxels xi and xj , such that w(xi , xj )  [ 0 , 1 ]   and xjviw(xi , xj ) = 1 ( cf . , the computation of the weight is based on the 
 euclidean distance between patches u(nj ) and u(ni ) , defined as 
			
 ( 2)w(xi , xj)=1zieu(ni)u(nj)22/h2 , 
 where zi is a normalization constant ensuring that jw(xi , xj ) = 1 , and h acts as a 
 filtering parameter controlling the decay of the exponential   function . 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 h first , it has been shown that the optimal smoothing parameter h is proportional to the standard deviation of the noise  . second , if we want the filter independent of the neighborhood size , the optimal h must depend on |ni| 
  ( see , ( 2 ) ) . thus , the automatic tuning of the filtering parameter h amounts to determining the relationship h = f (  , |ni| ,  ) , 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 i as defined in [ 18 , 19 ] . for each voxel xi 
  of the volume , let us define 
			
 ( 3)i=67(u(xi)16xjpiu(xj ) ) , 
 pi 
  being the 6-neighborhood at voxel xi 
  and the constant 6/7 is used to ensure that e[i2]=^2 in the homogeneous areas . thus , the standard deviation of noise ^ is computed as 
			
 ( 4)^2=1|3|i3i2 . then , in order to make the filter independent of |ni| , we used the euclidean distance || 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 ) can be written as 

 ( 6)w(xi , xj)=1zieu(ni)u(nj)22/2^2|ni| , 
 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 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 . a partition of 
 the volume  into 
 overlapping blocks bik 
  of size ( 2 + 1 ) is performed , 
 such as  = kbik 
 , under the constraint that each block bik 
  intersects with 
 at least one other block of the partition . these blocks are centered on voxels xik 
  which 
 constitute a subset of . the voxels xik 
  are equally 
 distributed at positions ik = ( k1 n , k2 n , k3 n ) , ( k1 ,   k2 ,   k3 )   , where n represents the 
 distance between the centers of bik 
 . to ensure a global continuity in the denoised image , 
 the overlapping support of blocks is nonempty : 2  n.for each block bik 
 , an nl - means - like restoration is performed as 
 follows : 
			 ( 7)nl(u)(bik)=bjvikw(bik , bj)u(bj),with  w(bik , bj)=1zikeu(bik)u(bj)22/2^2|ni| , 
				 where zik 
  is a 
 normalization constant ensuring that jw(bik 
 , bj ) = 1 ( see figure 1 , 
 right ) . for a voxel xi included in several 
 blocks bik 
 , several estimations of the restored intensity nl(u)(xi ) the estimations given by different nl(u)(bik 
 ) for a voxel xi are stored in a 
 vector ai . the final restored intensity of voxel xi is then defined 
 as(8)nl(u)(xi)=1|ai|paiai(p ) , 
 where ai(p ) denotes the pth element of the 
 vector ai . a partition of 
 the volume  into 
 overlapping blocks bik 
  of size ( 2 + 1 ) is performed , 
 such as  = kbik 
 , under the constraint that each block bik 
  intersects with 
 at least one other block of the partition . these blocks are centered on voxels xik 
  which 
 constitute a subset of . the voxels xik 
  are equally 
 distributed at positions ik = ( k1 n , k2 n , k3 n ) , ( k1 ,   k2 ,   k3 )   , where n represents the 
 distance between the centers of bik 
 . to ensure a global continuity in the denoised image , the overlapping support of blocks is nonempty : 2  n. for each block bik 
 , an nl - means - like restoration is performed as 
 follows : 
			 ( 7)nl(u)(bik)=bjvikw(bik , bj)u(bj),with  w(bik , bj)=1zikeu(bik)u(bj)22/2^2|ni| , 
				 where zik 
  is a 
 normalization constant ensuring that jw(bik 
 , bj ) = 1 ( see figure 1 , 
 right ) .   for a voxel xi included in several 
 blocks bik 
 , several estimations of the restored intensity nl(u)(xi ) the estimations given by different nl(u)(bik 
 ) for a voxel xi are stored in a 
 vector ai . the final restored intensity of voxel xi is then defined 
 as(8)nl(u)(xi)=1|ai|paiai(p ) , 
 where ai(p ) denotes the pth element of the 
 vector ai . the main advantage of this approach is to significantly reduce the complexity of the algorithm . indeed , for a volume  of size n , the global complexity is ((2 + 1 ) ( 2 m + 1 ) ( ( n  n)/n ) ) . for instance , with n = 2 , the complexity is divided by a factor 8 . in [ 2123 ] , 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(bj ) to the current patch u(bi ) to compute nl(u)(bi ) can be viewed as a spatially adaptation of the patch dictionaries . as in [ 2123 ] , the preselection of blocks in vi 
  is based on the mean and the variance of u(bi ) and u(bj ) . the selection tests are given by 
			
 ( 9)w(bik , bj)={1zikeu(bik)u(bj)22/2^2|ni|    if  1<u(bik)u(bj)<11 ,    12<var(u(bik))var(u(bj))<112,0  otherwise , where 
 u(bik) and var(u(bik 
 ) ) represent , respectively , the mean and the variance of the intensity function for the block bik 
  centered on the voxel xik 
 . 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 1are 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 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 . recently , hybrid approaches coupling the nl - means filter and a wavelet decomposition have been proposed [ 2 , 24 , 25 ] . in , 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 , 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 , 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 . in order to improve the denoising result of the nl - means filter , we propose a multiresolution framework similar to   to implicitly adapt the filtering parameters ( h , |bi| ) 
  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 . 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 can not completely remove all frequencies of noise . as a consequence , we propose the following workflow ( see figure 3 ) . denoising of 
 the original image i using two sets 
 of filtering parameters : one adapted to the noise components removal ( i.e. , 
 large blocks and  = 1 ) and the other 
 adapted to the image features preservation ( i.e. , small blocks and    = 0.5 ) . the noise is efficiently removed and , conversely , in iu , the image features are preserved.decomposing io and iu into low- and 
 high - frequency subbands . the first level decomposition of the images is 
 performed with a 3d discrete wavelet transform ( dwt).mixing the 
 highest - frequency subbands of io and the lowest 
 frequency subbands of iu.reconstructing 
 the final image by an inverse 3d dwt from the combination of the selected high 
 and low frequencies . denoising of 
 the original image i using two sets 
 of filtering parameters : one adapted to the noise components removal ( i.e. , 
 large blocks and  = 1 ) and the other 
 adapted to the image features preservation ( i.e. , small blocks and    = 0.5 ) . the noise is efficiently removed and , conversely , in iu , the image features are preserved . the first level decomposition of the images is 
 performed with a 3d discrete wavelet transform ( dwt ) . mixing the 
 highest - frequency subbands of io and the lowest 
 frequency subbands of iu . reconstructing 
 the final image by an inverse 3d dwt from the combination of the selected high 
 and low frequencies . 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 . 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 io and io images . for each image , eight subbands are obtained : lll1 , llh1 , lhl1 , hll1 , lhh1 , hlh1 , hhl1 , and hhh1 . in the eight wavelet subbands obtained with io , the frequencies corresponding to noise are 
 efficiently removed from the high frequencies whereas the low frequencies associated 
 to the main features are spoiled.in the eight 
 wavelet subbands obtained with iu , the low frequencies associated to main features are 
 efficiently preserved whereas residual frequencies corresponding to noise are 
 present in high frequencies . in the eight wavelet subbands obtained with io , the frequencies corresponding to noise are 
 efficiently removed from the high frequencies whereas the low frequencies associated 
 to the main features are spoiled . in the eight 
 wavelet subbands obtained with iu , 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 io 
 ( i.e. , lhh1 , hlh1 , hhl1 , and hhh1 ) and the lowest frequencies of iu ( i.e. , lll1 , llh1 , lhl1 , and hll1 ) . then , the 4 lowest subbands of iu are combined with the 4 highest subbands of io . 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 iu and io were su = ( u , mw , u ) = ( 1 , 3 , 0 , 5 ) and so = ( o , mw , o ) = ( 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 ) . 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 . 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  ( 0 , (3/100 ) ) , where  is the value of the highest voxel intensity of the phantom ( 150 for t1-w and 250 for t2-w ) . in the following , nl - means : 
 standard voxelwise implementation with automatic tuning of the filtering 
 parameter h (  = 1 ) .optimized nl - means : voxelwise implementation with automatic tuning of the filtering 
 parameter h (  = 1 ) and voxels 
 selection presented in . 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 . 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 . 
		 for the nl - means and optimized nl - means filters , the parameters are those used in 
 : d = 1 ,  = 1 , m = 5 , 1 = 0.95 and 1 = 0.5.concerning the optimized 
 blockwise nl - means filter , the sets of parameters are those used in 
 : n = 2 ,  = 1 ,  = 1 , m = 5 , 1 = 0.95 and 1 = 0.5.finally , for 
 the optimized blockwise nl - means with wavelet mixing filter the 
 parameter are the following : n = 2 , su = ( u , mw , u ) = ( 1 , 3 , 0.5 ) , so = ( o , mw , o ) = ( 2 , 3 , 1 ) , 1 = 0.95 , and 1 = 0.5 . 
 for 8-bit encoded images , the psnr is defined as follows : 

 ( 10)psnr=20  log10255rmse , 
 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 ) . nl - means : 
 standard voxelwise implementation with automatic tuning of the filtering 
 parameter h (  = 1 ) . optimized nl - means : voxelwise implementation with automatic tuning of the filtering 
 parameter h (  = 1 ) and voxels 
 selection presented in . 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 . 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 . for the nl - means and optimized nl - means filters , the parameters are those used in 
 : d = 1 ,  = 1 , m = 5 , 1 = 0.95 and 1 = 0.5 . concerning the optimized 
 blockwise nl - means filter , the sets of parameters are those used in 
 : n = 2 ,  = 1 ,  = 1 , m = 5 , 1 = 0.95 and 1 = 0.5 . finally , for 
 the optimized blockwise nl - means with wavelet mixing filter the 
 parameter are the following : n = 2 , su = ( u , mw , u ) = ( 1 , 3 , 0.5 ) , so = ( o , mw , o ) = ( 2 , 3 , 1 ) , 1 = 0.95 , and 1 = 0.5 . 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 ) . 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 . 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 . 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 . in this section , we compare the proposed method with two of the most used approaches in mri domain : the nonlinear diffusion ( nld ) filter r   and the total variation ( tv ) minimization . 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 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 su = ( u , mw , u ) = ( 1 , 3 , 0.5 ) and so = ( o , mw , 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 . 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 the psnr value between the noisy image and the ground truth is called  no processing  and is used as a reference . figure 10 shows the denoising results obtained by the nld filter , the tv minimization , and our optimized blockwise nlm with wm . 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 . the t1-weighted mr images used for experiments were obtained with t1 sense 3d sequence on 3 t 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 . this paper presented a fully automated blockwise version of the nonlocal means filter with subbands wavelet mixing . 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   and the total variation minimization .
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 .