|Home | About | Journals | Submit | Contact Us | Français|
A critical issue in image restoration is the problem of noise removal while keeping the integrity of relevant image information. Denoising is a crucial step to increase image quality and to improve the performance of all the tasks needed for quantitative imaging analysis. The method proposed in this paper is based on a 3D optimized blockwise version of the Non Local (NL) means filter . The NL-means filter uses the redundancy of information in the image under study to remove the noise. The performance of the NL-means filter has been already demonstrated for 2D images, but reducing the computational burden is a critical aspect to extend the method to 3D images. To overcome this problem, we propose improvements to reduce the computational complexity. These different improvements allow to drastically divide the computational time while preserving the performances of the NL-means filter. A fully-automated and optimized version of the NL-means filter is then presented. Our contributions to the NL-means filter are: (a) an automatic tuning of the smoothing parameter, (b) a selection of the most relevant voxels, (c) a blockwise implementation and (d) a parallelized computation. Quantitative validation was carried out on synthetic datasets generated with BrainWeb . The results show that our optimized NL-means filter outperforms the classical implementation of the NL-means filter, as well as two other classical denoising methods (Anisotropic Diffusion  and Total Variation minimization process ) in terms of accuracy (measured by the Peak Signal to Noise Ratio) with low computation time. Finally, qualitative results on real data are presented.
Quantitative imaging involves image processing workflows (registration, segmentation, visualization, etc.) with increasing complexity and sensitivity to possible image artifacts. As a consequence, image processing procedures often require to remove image artifacts beforehand in order to make quantitative post-processing more robust and efficient. A critical issue concerns the problem of noise removal while keeping the integrity of relevant image information. This is particularly true for ultrasound images or magnetic resonance images (MRI) in presence of small structures with signals barely detectable above the noise level. In addition, a constant evolution of quantitative medical imaging is to process always larger cohorts of 3D data in order to find significant discriminants for a given pathology (e.g. see ). In this context, complex automatic image processing workflows are required  since human interpretation of images is no longer feasible. For effectiveness, these workflows have to be robust to a wide range of image qualities and parameter-free (or at least using auto-tuned parameters). This paper focuses on these critical aspects by introducing a new restoration scheme in the context of 3D medical imaging. The Non Local (NL-) means filter was originally introduced by Buades et al.  for 2D image denoising. The adaptation of this filter we propose for 3D images is based on (a) an automatic tuning of the smoothing parameter, (b) a selection of the most relevant voxels for the NL-means computation, (c) a blockwise implementation and (d) a parallelized computation. These different contributions allow to make the adapted filter fully-automated and above all to overcome the main limitation of the classical NL-means: the computational burden.
Section II gives a short overview of the literature on image restoration. Section III presents the proposed method with details about our contributions. Sections IV, V and VI show (a) the impact of our adaptations compared to the classical NL-means implementation and (b) a comparison with respect to other well established denoising methods on Gaussian and Rician noise. Both validation experiments are performed on a phantom data set in a quantitative way. Section VII shows results on real data such as 3 Tesla T1-weighted (T1-w) MRI and T2-weighted (T2-w) MRI of a patient with Multiple Sclerosis (MS) lesions. In section VIII we propose a discussion on the applicability and the further improvements of the NL-means filter in the context of 3D medical imaging.
Many methods have been proposed for edge-preserving image denoising. Some popular approaches include Bayesian approaches , PDE-based approaches , , , , robust and regression estimation , adaptive smoothing , wavelet-based methods –, bilateral filtering –, local mode filtering , hybrid approaches –.
Strong theoretical links exist between most of these techniques, as recently shown for local mode filtering , bilateral filtering, anisotropic diffusion and robust estimation ,  and adaptive smoothing , anisotropic diffusion and total variation minimization scheme .
More recently, some promising methods have been proposed for improved image denoising, based on statistical averaging schemes enhanced via incorporating a variable spatial neighborhood scheme –. Other approaches consist in modeling non-local pairwise interactions from training data  or a library of natural image patches , . The idea is to improve the traditional Markov random field (MRF) models by learning potential functions from examples and extended neighborhoods for computer vision applications –. Awate and Whitaker proposed another non-parametric patch-based method relying on comparisons between probability density functions .
Some of these techniques, generally developed for 2D images, have often been extended to 3D medical data, especially to MR images: anisotropic diffusion , , total variation , bilateral filtering and variants , wavelet-based filtering –, hybrid approaches , .
Most of the denoising methods restore the intensity value of each image voxel by averaging in some way the intensities of its (spatially) neighboring voxels. The basic and intuitive approach is to replace the value of the voxel by the average of the voxels in its neigbourhood (so-called box filtering ). In practice, this filter has been shown to be outperformed by the Gaussian filter, which consists in weighting each voxel in the neighborhood according to its distance to the voxel under study. Both filters can be iterated until the desired amount of smoothing is reached. Such data-independent approaches can be implemented very efficiently. Their major drawback is that they blur the structures of interest in the image (e.g. edges or small structures and textures).
This has naturally led to data-dependent approaches, which aim at eliminating (or reducing the influence of) the neighboring voxels dissimilar to the voxel under study. Simple order statistic operators can be used for this purpose, such as the median filter, leading to a simple generalization of the box filter. More sophisticated approaches, based on image derivatives have been successfully proposed for many applications, such as adaptive smoothing  and anisotropic diffusion . Neighborhood filters ,  and variants , , have been also proposed and consist in averaging input data over the image voxels that are spatially close to the voxel under study and with similar gray-level values.
All these techniques rely on the idea that the restored value of a voxel should only depend on the voxels in its spatial neighborhood that belong to the same population, that is the same image context. This has been termed by Michael Elad as the locally adaptive recovery paradigm . Another approach has been recently proposed, that has shown very promising results. It is based on the idea that any natural image has redundancy, and that any voxel of the image has similar voxels that are not necessarily located in a spatial neighborhood. First introduced by Buades et al. in , the NL-means filter is based on this redundancy property of periodic images, textured images or natural images to remove noise. In this approach, the weight involving voxels in the average, does not depend on their spatial proximity to the current voxel but is based on the intensity similarity of their neighborhoods with the neighborhood of the voxel under study, as in patched-based approaches. In other words, the NL-means filter can be viewed as an extreme case of neighborhood filters with infinite spatial kernel and where the similarity of the neighborhood intensities is substituted to the point-wise similarity of gray levels as in commonly-used bilateral filtering. This new non-local recovery paradigm allows to combine the two most important attributes of a denoising algorithm: edge preservation and noise removal.
In the following, we introduce the notations:
In the classical formulation of the NL-means filter, the restored intensity NL(u)(xi) of the voxel xi, is the weighted average of all the voxel intensities in the image u defined as:
where u(xj) is the intensity of voxel xj and w(xi, xj) is the weight assigned to u(xj) in the restoration of voxel xi. More precisely, the weight quantifies the similarity of the local neighborhoods Ni and Nj of the voxels xi and xj under the assumptions that w(xi, xj) [0, 1] and Σxj Ω3 w(xi, xj) = 1 (cf Fig. 1 left). The classical definition of the NL-means filter considers that each voxel can be linked to all the others, but for practical computational reasons the number of voxels taken into account in the weighted average can be limited to the so-called “search volume” Vi of size (2M+1)3, centered at the current voxel xi.
For each voxel xj in Vi, the Gaussian-weighted Euclidean distance defined in , is computed between u(Nj) and u(Ni). This distance is a classical L2 norm convolved with a Gaussian kernel of standard deviation a, and measures the distance between neighborhood intensities. Given this distance, w(xi, xj) is computed as follows:
where Zi is a normalization constant ensuring that Σj w(xi, xj) = 1, and h acts as a smoothing parameter controlling the decay of the exponential function. When h is very high, all the voxels xj in Vi will have the same weight w(xi, xj) with respect to the voxel xi. The restored value NL(u)(xi) will be then approximately the average of the intensity values of the voxels in Vi leading to strong smoothing of the image. When h is very low, the decay of the exponential function will be strong, thus only few voxels xj in Vi with u(Nj) very similar to u(Nj) will have a significant weight. The restored value NL(u)(xi) will tend to be the weighted average of some voxels with a similar neighborhood to current voxel xi leading to a weak smoothing of the image. In Section III-B.1, a trade-off has then to be found, and we propose a method to automatically estimate the optimal value of h.
In , Buades et al. show that, for 2D natural images, the NL-means filter outperforms state-of-the-art denoising methods such as the Rudin-Osher-Fatemi Total Variation minimization scheme , the Perona-Malik Anisotropic diffusion  or translation invariant wavelet thresholding . Nevertheless, the main drawback of the NL-means filter is the computational burden due to its complexity, especially for 3D images. Indeed, for each voxel of the volume, distances between the intensity neighborhoods u(Ni) and u(Nj) for all the voxels xj contained in Vi need to be computed. Let N3 denote the size of the 3D image, then the complexity of the filter is in the order of ((N(2M+1)(2d+1))3). For a 3D MRI of size 181 × 217 × 181 with the smallest possible value for d and M (d = 1 and M = 5), the computational time reaches up to 6 hours on 3 GHz CPU. This time is far beyond a reasonable duration expected for a denoising filter in a medical practice. For this reason, we propose several adaptations to reduce the computational burden which are detailed in Section III-B. We also show that these adaptations improve the quality of the denoising compared to the classical implementation.
According to , the smoothing parameter h depends on the standard deviation of the noise σ, and typically a good choice for 2D images is h ≈ 10σ. Equation 2 shows that h also needs to take into account |Ni|, if we want a filter independent of the neighborhood size. Indeed, the L2 norm increasing with |Ni|, h needs also to be increased to obtain an equivalent filter. The automatic tuning of the smoothing parameter h comes to determine the relationship h2 = f(σ2, |Ni|, β) where β is a constant. Let us show how we can estimate this relationship:
To deal with computational burden, Mahmoudi and Sapiro  recently proposed a method to preselect a subset of the most relevant voxels xj in Vi to avoid useless weight computations. In other words, the main idea is to select only the voxels xj in Vi that will have the highest weights w(xi, xj) in Equation 1 without having to compute all the Euclidean distances between u(Ni) and u(Nj). A priori neglecting the voxels which are expected to have small weights tends to speed up the filter and to improve the results (see Table II). In , this selection is based on the similarity of the mean value of u(Ni) and u(Nj), and on the similarity of the average over the neighborhoods Ni and Nj of the gradient orientation at pixel xi and xj. Intuitively, similar neighborhoods have the same mean and the same gradient orientation. The computation of the gradient orientation is very sensitive to noise and thus requires robust estimation techniques. This is too computationally expensive for medical applications. For this reason, in our implementation, the preselection of voxels in Vi is based on the mean and the variance of u(Ni) and u(Nj) which allows to decrease the computational burden. Figure 2 shows that the maps of local means and local variances are simple estimators allowing to discriminate different tissue classes and edges in images. In this way, the maps of local means and local variances are precomputed in order to avoid repetitive calculations for the same neighborhood. The selection tests can be expressed as follows:
where and Var(u(Ni)) represents respectively the mean and the variance of the local neighborhood Ni of voxel xi. As suggested in , with this kind of selection, the NL-means filter tends to better preserve the detailed regions while slightly spoiling the denoising of the flat regions. Indeed, in flat regions increasing the number of voxels tends to improve denoising because there are a large number of similar voxels. In more cluttered regions, increasing the number of voxels tends to remove the details during smoothing because there are very few similar voxels.
A blockwise implementation of the NL-means is developed as suggested in . This approach 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.
The main advantage of this approach is to significantly reduce the complexity of the algorithm. Indeed, for a volume Ω3 of size N3, the global complexity is . For instance, with n = 2, the complexity is divided by a factor 8. The voxels selection principle can also be applied to the blockwise implementation:
where and Var(u(Bik)) represent respectively the mean and the variance of the intensity function, for the block Bik centered on the voxel xik.
Another way to reduce the computational time is to distribute the operations on several processors via a cluster or a grid. The intrinsic nature of the NL-means filter makes it perfectly suited for parallelization and multithreading implementation. One of the main advantage of this filter, when compared to others method such as Anisotropic Diffusion or Total Variation minimization, is that the operations are performed without any iterative schemes. Thus, the parallelization of the NL-means filter is straightforward to perform and very efficient. We divide the volume into sub-volumes, each of them being treated separately by one processor. A server with 8 Xeon processors at 3 GHz and a Intel(R) Pentium(R) D CPU 3.40GHz were used in our experiments.
In order to evaluate the performances of the NL-means filter on 3D MR images, tests were performed on the BrainWeb database1 . Two images were simulated: T1-w MR image using SFLASH sequence (volume size = 181 × 217 × 181) and T2-w MR image with MS from SFLASH sequence (volume size = 181 × 217 × 181). As reported previously, it is a known fact that the MR images are corrupted by a Rician noise , , which can be well approximated by a white Gaussian noise in high intensity areas, typically in brain tissues . In order to verify if this approximation can be used for a NL-means based denoising, experiments are performed on phantom images with Gaussian and Rician noise.
A white Gaussian noise was added on the “ground truth”, and the notations of BrainWeb are used: a noise of 3% is equivalent to , where ν is the value of the brightest tissue in the image (150 for T1-w and 250 for T2-w). Several images were simulated to validate the performances of the denoising on various images (see Fig. 4):
T2-w images were used in order to show that our approach and its calibration are not specific to T1-w MRI sequences. Moreover, the tests on T2-w MRI with MS show how the NL-means filter could be useful in a pathological context due to its preservation of anatomic and pathologic structures.
The Rician noise was build from white Gaussian noise in the complex domain. Firstly, two images are computed:
where I0 is the “ground truth” and σ is the standard deviation of the added white Gaussian noise. Then, the noisy image is computed as:
The notation 3% for the Rician noise means that the Gaussian noise used in complex domain is equivalent to , where ν is the value of the brightest tissue in the image (150 for T1-w). According to the Peak Signal to Noise Ratio (PSNR) (see Eq. 13) between “ground truth” and noisy images, for a same level of noise, the Rician noise is stronger than the Gaussian noise (see Tab. I). Several images were simulated (see Fig. 5):
To show the efficiency of the NL-means filter on real data, tests were performed on image acquired with a high field MR system (Bruker 3 Tesla). The data used was a 256 × 256 × 256 T1-w image.
In a pathological context, the denoising step is crucial especially when the structures of interest have a small size: the integrity of pathological structures must be preserved by the denoising method. As said earlier, one objective of denoising is to include such processing in complex medical imaging workflows. This kind of workflows is widely used to process large cohort of subjects in many neurological diseases such as MS lesions. The data used for MS lesions qualitative validation was a T2-w MR image from an axial dual-echo, turbo spin-echo sequence (Philips 1.5 Tesla).
In the following, let us define:
In this section different aspects of NL-means filter implementation were investigated. First, the impact of the automatic tuning of the filtering parameter (Section V-A) and the influence of the size of the search volume and the neighborhood were studied (Section V-B). Then, the impact of voxels selection and blockwise implementation is analyzed via the comparison of the NL-means, Optimized NL-means, Blockwise NL-means and Optimized Blockwise NL-means filters (Sections V-C and V-D). Finally, we compare the proposed Optimized Blockwise NL-means filter with other well-established denoising methods: Anisotropic Diffusion filter  (implemented in VTK2) and Rudin-Osher-Fatemi Total Variation (TV) minimization process  (3D extension of the Megawave 2 implementation3) (Section V-G). The different variants of the NL-means filter can be freely tested at: http://www.irisa.fr/visages/benchmarks
In the following, several criteria are used to quantify the performances of each method: the PSNR obtained for different noise levels, histogram comparisons between the denoised images and the “ground truth”, and finally visual assessment. For the sake of clarity, the PSNR and the histograms are estimated only in a region of interest obtained by removing the background. For 8-bit encoded images, the PSNR is defined as follows:
where the RMSE is the root mean square error estimated between the ground truth and the denoised image.
In this section, the central parameters of interest are:
In each experiment (Sections V-A, V-B and V-C), we let one parameter vary while remaining the others constant, with default values: β = 1, M = 5, d = 1 and μ1 = 0.95, . Concerning the blockwise implementation the default parameters are n = 2 and α = 1.
Our experiments have shown that all the versions of the NL-means filter (NL-means, Optimized NL-means, Blockwise NL-means and Optimized Blockwise NL-means) tend to have a similar behavior with respect to the variation of the parameters. In that context, all the results are displayed with the proposed Optimized Blockwise NL-means filter, even if equivalent conclusions can be drawn with the NL-means, Optimized NL-means and Blockwise NL-means filters.
Validation was performed on T1-w and T2-w MRI, but the results concerning the study of the parameter influences are shown for T1-w MRI only. The results on T2-w MRI are shown in section V-F in order to underline that the parameters calibrated for T1-w MRI work fine on T2-w MRI.
Figure 6 shows the influence of the automatic determination of the smoothing parameter h2 = 2β2|Ni|. As described in III-B.1, h is a function of the global standard deviation of the noise in the volume esimated from pseudo-residuals (see 3 and Eq. 4). Here, β allows to adjust the automatic estimation of h in order to determine the optimal smoothing parameter 2β2|Ni| for each level of noise (see 6). For low levels of noise the best value of β is close to 0.5. For high levels of noise this value is 1. These results show that the estimation of the standard deviation of the noise is correctly performed by pseudo-residuals. These observations underline (a) how efficient the automatic estimation of the smoothing parameter h is, and (b) how the NL-means can be used without manual parameter tuning.
Figure 7 shows the influence of the size of the search volume and the local neighborhood. Increasing the number of voxels in the search volume Vi does not seem to affect the PSNR when M is greater than 5. Indeed, the theoretic definition of the NL-means filter states that the weighted average (see Eq. 1) computed for the restoration of voxel xi should be performed on all voxels xj Ω3. Practically, the limit M = 5 prevents useless computations. Moreover, increasing d degrades the denoising process. When d increases the NL-means filter drastically slows down. That is why we have not investigated the impact of d for d > 3.
The selection of the voxels in the search volume Vi is achieved by supposing that only the voxels whose the neighborhood is similar to the neighborhood of the voxel under study could be considered (see Eq. 1). To do so, as defined in III-B.2, the weight w(xi, xj) is calculated only for voxels such that: and . The influence of the limits μ1 and is studied in Figure 8. In a first experiment μ1 varies according to γ such as μ1 = 1 − γ while . In a second experiment varies according to γ following while μ1 = 0.95.
Figure 8 (left) shows that a restrictive selection based on the mean (low values of γ) increases the PSNR. In other words, the number of voxels taken into account in the weighted average is drastically reduced, as well as the computational time (also see Tab. II). The optimal limits were obtained for μ1 = 0.95 while . Concerning the variance (Figure 8, right), we observe that a too restrictive selection degrades the PSNR. In addition, a too permissive selection does not increase the PSNR while increasing uselessly the computational burden. A compromise was found by fixing . There is a clear dependency between the bounds for the mean and the variance. An optimal trade-off was determined experimentally.
Tab. II shows that the blockwise approach of the NL-means filter, with and without voxels selection (see Eq. 11), allows to drastically reduce the computational time. With a distance between the block centers n = 2, the blockwise approach divides this time by a factor 23 = 8 (see Tab. II). However, computational time reduction needs to be balanced with a slight decrease of the PSNR (see Fig. 9, left). For the optimized versions, the voxels/blocks selection in the search volume has several impacts. First, by reducing the average number of voxels/blocks used in the weighted averages, this decreases the computational time compared to the non-optimized versions (see Tab. II). Second, the selection of the most relevant voxels/blocks increases the quality of denoising for all the noise levels (see Fig. 9 (left) and Tab. II).
As described in section III-B.4, the multithreading in the case of the NL-means filter is particularly adapted due its non iterative nature. For the classical pixelwise NL-means implementation, the parallelization allows to divide the computational time by a factor close to the numbers of CPU. As eight processors were used for our experiments, the computational time with multithreading is about 8 times smaller (see Tab. II, and ). For the blockwise implementations the speedup is less and ). The difference of speedup between the classical NL-means and the blockwise NL-means filters have two origins.
In order to underline that the utilization of 8-CPUs is not required by our filter, the denoising have been also performed on a more common architecture: a DualCore Intel(R) Pentium(R) D CPU 3.40GHz. The results show that our filter takes less than 3 minutes to denoise a volume 181 × 217 × 181 voxels on this architecture.
To conclude, the different improvements included in the proposed Optimized Blockwise NL-means filter (i.e., blockwise approach and blocks selection) allow to speed up the denoising procedure, compared to NL-means filter, by a factor of 66 on 1 Xeon at 3GHz, 44 on 8× Xeon at 3GHz and 31 on a DualCore at 3.40GHz.
Figure 9 (right) presents the results obtained by the different NL-means filter versions on T2-w MRI with MS lesions. The optimal parameters (i.e. the default parameters described in section V), experimentally determined on T1-w MRI, and the automatic tuning of h were used on T2-w MRI. The Optimized NL-means and Optimized Blockwise NL-means filters outperform the NL-means and Blockwise NL-means filters also on T2-w MRI. The most important difference between the optimized and non-optimized versions are observed on T2-w MRI, which could be explained by the higher level of noise in the simulated T2-w MRI compared to T1- w MRI. Actually, the variance of noise varies with respect to the highest intensity tissues which is 150 in T1-w and 250 in T2-w. For 9% the variance of noise is 13.5 in T1-w images and is 22.5 in T2-w images because the highest tissue intensity is superior in T2-w images. These results suggest that the parameters experimentally tuned on T1-w images can be used for T2-w images. Figure 10 shows an example of denoising obtained by the optimized blockwise NL-means and the blockwise NL-means filters. The MS lesions are visually more preserved with the optimized version; this was confirmed by an experienced MRI reader.
As reported in Section II, the Anisotropic Diffusion filter (AD) was introduced to overcome the blurring effect of the Gaussian smoothing approach. First introduced by Perona and Malik , in this approach the image u is only convolved in the direction orthogonal to the gradient of the image which ensures the preservation of edges. The iterative denoising process of initial image u0 can be expressed as:
where u(x, t) is the image gradient at voxel x and iteration t, is the partial temporal deviation of u(x, t) and
where K is the diffusivity parameter. The AD filter method produces a good preservation of edges , . Nonetheless, the main disadvantage of Ad filter is to poorly denoise the constant regions (see Fig. 13).
The difficult task to preserve edges while correctly denoising constant areas has been addressed also by Rudin, Osher and Fatemi. They proposed to minimize the Total Variation (TV) norm subject to noise constraints , that is:
where u0 is the original noisy image, u is the restored image and σ the standard deviation of the noise. In this model, the TV minimization tends to smooth inside the image structures while keeping the integrity of boundaries. The TV minimization scheme can be expressed as an unconstrained problem:
where λ is a Lagrange multiplier which controls the balance between the TV norm and the fidelity term. Thus, λ acts as the filtering parameter. Indeed, for high values for λ the fidelity term is encouraged. For small values for λ the regularity term is desired. In practice, the TV minimization scheme tends to remove texture and small image structures as seen in Fig. 13 . To solve this problem, iterative total variation schemes have been recently developed , .
The main difficulty to achieve this comparison is related to the tuning of smoothing parameters in order to obtain the best results for AD filter and TV minimization scheme. In order not to penalize AD filter and TV minimization scheme, an exhaustive search for all parameters into a certain range. Then, the best results obtained with AD filter and TV minimization have been selected, whereas the fully-automatic results have been mentioned for the NL-means filters. The results of the NL-means filters are not “optimal” due to the non perfect estimation of the noise standard deviation. For AD filter, the parameter K varies from 0.05 to 1 with a step of 0.05 and the number of iterations varies from 1 to 10. For TV minimization, the parameter λ varies from 0.01 to 1 with step of 0.01 and the number of iterations varies from 1 to 10. The results obtained for 9% of Gaussian noise are presented Fig. 11, but this screening have been 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 necessary give the best visual output. Indeed, the best PSNR for AD filter is obtained for a visually under-smoothed image since this method tends to spoil the edges. To obtain a high PSNR, the denoised image needs to balance edge preserving and noise removing. For AD filter, this trade-off leads to inhomogeneities in flat areas in denoised image (see Fig. 13). For TV minimization, the best PSNR is obtained with a visually under-smoothed image since noise is present in denoised image (see Fig. 13).
As presented in Fig. 12 (left), our block optimized NL-means filter produces the best PSNR values whatever the noise level. On average, a gain of 1.85 dB is measured compared to TV and Anisotropic Diffusion methods. The PSNR value between the noisy image and the ground truth is called “No processing” and is used as reference.
To better understand how these differences in the PSNR values between the three methods can be explained, the histograms of the denoised images were compared to the histogram of the ground truth. Figure 12 (right) shows that the Optimized Blockwise NL-means filter is the only method able to retrieve a histogram similar to the ground truth. The NL-means-based restoration schemes clearly distinguish the three main peaks representing the white matter, the gray matter and the cerebrospinal fluid. The sharpness of the peaks shows how the Optimized Blockwise NL-means filter increases the contrast between denoised biological structures (see also Fig. 13). The distances between these histograms are estimated with the Bhattacharyya coefficient (BC) defined as:
where p and q are the two histograms be to compared and b is the bin index. A BC close to 1 means p and q are very similar. Each histogram of denoised images is compared to the “ground truth” one (see Tab. III). The distance between the histogram of the noisy image and the histogram of the “ground truth” is used as a reference. The BC distance shows that the restored histogram obtained with the Optimized Blockwise NL-means filter is the closest to the “ground truth”, as visually assessed in Figure 12 (right). Finally, Table III suggests that the NL-means-based approach could improve the registration of images, since the Mutual Information (MI) computed between the restored image and the “ground truth” is higher in comparison with AD filter and TV minimization. The MI is a similarity measure commonly used in image registration.
Figures 13 and and1414 show the restored images and the removed noise obtained with the three compared methods. As shown in the previous analysis, we observe that the homogeneity of white matter is higher when the image is denoised with the Optimized Blockwise NL-means filter. Moreover, focusing on the structure of the removed noise, it clearly appears that NL-means-based restoration schemes better preserves the high frequency components of the image corresponding to anatomical structures while removing efficiently the high frequencies due to noise. According to the “method noise” introduced in , the NL-means is a better denoising method since the removed noise is the most similar to a white Gaussian noise. Finally, the difference between the “ground truth” and the denoised image is presented in order to show which structures are removed during the denoising process. In Fig. 13, this difference shows that (a) the AD filter tends to spoil the edges especially on the skull, (b) the TV minimization slightly better preserves the edges but does not remove all the noise, and (c) the Optimized Blockwise NL-means filter visually better preserves the edges while efficiently removing the noise (especially for white matter).
In this section, the same experiments are performed on phantom data set corrupted by Rician noise in order to study the impact of the Gaussian assumption. Table IV shows the computation time and the denoising performance of the different compared NL-means filters. These results show that the optimized NL-means versions outperform the classical ones also for Rician noise. Figure 15 presents the comparison with AD filter and TV minimization in terms of PSNR values and histogram analysis. As for the AD filter and TV minimization, the NL-means-based denoising is able to correctly restore an image corrupted by Rician noise using a Gaussian approximation. When the histograms are compared, low values of intensity (< 20) are incorrectly restored for all the filters; the Gaussian approximation is not appropriate in that case. Nevertheless, it seems the underlying assumption is well suited to high values (> 60).
As for Gaussian noise, the NL-means-based restoration clearly emphasizes the three main peaks corresponding to the white matter, the gray matter and the cerebrospinal fluid. Figure 16 shows the visual results obtained when the methods are compared on phantom data with Rician noise. Compared to Fig. 13, the denoising of background is worse in the Rician case, but the cerebral structures are correctly restored with the NL-means filter especially the white matter (see Fig. 16). Finally, Figure 17 shows the PSNR results of the parameter screening for the AD filter and the TV minimization at 9% of Rician noise. All these results on Rician noise show that the PNSR values slightly decrease due to more pronounced noise compared to Gaussian case for a same level (see IV-A.2 for explanation), but the general performance of the filters is preserved.
The restoration results presented in Fig. 18 show good preservation of the cerebellum contours. Fully automatic segmentation and quantitative analysis of such structures are still a challenge, and improved restoration schemes could greatly improve these processings.
Figure 19 shows that the optimized blockwise NL-means filter preserves the lesions while removing the noise. The impact on further processing is not the scope of this paper and is not studied here. Nevertheless, visually the lesions appears more contrasted and as seen on the difference image the proposed NL-means approach does not include any structure of lesion in the estimated noise image. This was confirmed by an experienced neurologist.
This paper presents an optimized blockwise version of the Non Local (NL-) means filter, applied to 3D medical data. Validation was performed on the BrainWeb dataset  and showed that the proposed Optimized Blockwise NL-means filter outperforms the classical implementation of the NL-means filter and some state-of-the-art techniques, such as the Anisotropic Diffusion approach  and the Total Variation minimization process  on both Gaussian and Rician noise. These first results show that the image-redundancy assumption required for NL-means based restoration holds for 3D MRI. Compared to the classical NL-means filter, our implementation (with voxel preselection, multithreading and blockwise implementation) considerably decreases the required computational time (up to a factor of 60 on a Xeon at 3GHz) and increases the PSNR value of the denoised image. Nevertheless, the problem of the computational burden can still be investigated with other faster implementations such as the “plain multiscale” scheme also suggested in . Further works should be pursued for comparing NL-means based restoration with recent promising denoising methods, such as Total Variation in wavelet domain  or adaptive estimation method , . Moreover, the efficiency of the technique limiting the staircasing effect proposed in  needs to be studied for MRI.
We show on sample pathological cases (patients with MS lesions) that the filter preserves the major visual signature of the given pathology. However, the impact on specific pathologies needs to be further investigated.
Finally, the impact of the NL-means-based denoising on the performances of post-processing algorithms, like segmentation and registration schemes also should be studied. Nonetheless, the first results presented on the Mutual Information (MI) suggest that the proposed Optimized Blockwise NL-means filter could improve the image registration process. Indeed, the MI computed between the restored image and the “ground truth” is higher with the Optimized Blockwise NL-means filter than with the Anisotropic Diffusion approach and the Total Variation minimization process.