Computational and Mathematical Methods in Medicine
Volume 2015 (2015), Article ID 232389, 11 pages
http://dx.doi.org/10.1155/2015/232389

Research Article

A MRI Denoising Method Based on 3D Nonlocal Means and Multidimensional PCA

Liu Chang,1,2 Gao ChaoBang,1,2 and Yu Xi1,2

¹School of Computer Science, Chengdu University, Chengdu 610106, China
²Key Laboratory of Pattern Recognition and Intelligent Information Processing of Sichuan, Chengdu, China

Received 14 June 2015; Revised 12 August 2015; Accepted 24 August 2015

Academic Editor: Anne Humeau-Heurtier

Copyright © 2015 Liu Chang et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

Recently nonlocal means (NLM) and its variants have been applied in the various scientific fields extensively due to its simplicity and desirable property to conserve the neighborhood information. The two-stage MRI denoising algorithm proposed in this paper is based on 3D optimized blockwise version of NLM and multidimensional PCA (MPCA). The proposed algorithm takes full use of the block representation advantageous of NLM3D to restore the noisy slice from different neighboring slices and employs MPCA as a postprocessing step to remove noise further while preserving the structural information of 3D MRI. The experiments have demonstrated that the proposed method has achieved better visual results and evaluation criteria than 3D-ADF, NLM3D, and OMNLM_LAPCA.

1. Introduction

As a significant imaging technique, magnetic resonance imaging (MRI) provides very important information to research the tissues and organs in the human body with noninvasive style. However, MRI is affected by several artifacts and noise sources. One of them is the random fluctuation of the MRI signal which is mainly due to thermal noise. Such noise seriously degrades the acquisition of any quantitative measurements from the data. Consequently, the denoising techniques are required to improve the quality of MRI.

Generally speaking, MRI denoising techniques can be classified as either filtering, transform, or statistical approach [1]. Filtering methods remove noise with linear or nonlinear filters [2–5]. Transform methods employ some kinds of transformation to denoising MRI including wavelet transform [6] and curvelet transform [7]. Statistical methods estimate the noise with maximum likelihood [8], linear minimum mean square error (LMMSE) [9], Markov random process, and empirical Bayes [10]. In particular, nonlocal means (NLM) filter [11] has been used to denoise MRI image, achieving notable results [12–14]. NLM exploits the redundancy of the neighborhood pixel to remove the noise. The restored pixel is considered as the weighted average of the intensities of all pixels within the neighborhood area. Since MRI image has multichannel nature, NLM has been modified to denoise MRI data where the similarity measure can be considered to combine the relative information between different slices [15]. Nevertheless, the high computational burden has restricted its application for 3D MRI data. Therefore, [12] has proposed an optimized blockwise NLM filter for 3D MRI.

Due to its ability to perform decorrelation, PCA has also been used in image denoising. However, PCA requires that the number of images be bigger than the number of significant components of the image. The drawback has limited the application of PCA in the field of image denoising. The paper [16] has developed a two-stage approach to improve the quality of MRI data. After denoising with optimized multicomponent nonlocal mean (OMNLM), the local PCA is conducted over small local windows instead of the whole image to overcome the drawback. Nevertheless, PCA on overlapping windows will reduce the computational efficiency. Furthermore, the vectorization will make the structural information of image lost.

Actually, MRI is naturally a 3D image, which can be considered as tensor data on multidimensional space. From the aspect of super-resolution reconstruction, the noise image can be considered as the degraded version of the original image. Therefore, this paper proposed a multidimensional structure preserving MRI denoising algorithm. The algorithm consists of two stages. On the first stage, the 3D variant of the nonlocal means technique is employed to reduce the noise, which takes full advantage of the neighbor information between different 3D MRI slices and has the capability of exploiting the underlying structure in the multidimensional image. On the second stage, for the result image obtained from the first stage, multidimensional principal component analysis is performed to suppress the remaining noise, which avoids the vectorization to preserve the neighborhood information for MRI image and is helpful to improve the computation cost. According to the experiments on 3D MRI image, the proposed algorithm is superior to restore the original image from noise compared with other state-of-the-art methods.

2. Material and Methods

2.1. The 3D Variant of Nonlocal Mean

For image denoising problem, the noisy data is defined by the original noise free data with some noise :

The classical NLM technology has believed that the intensity of the point can be restored from the weighted average of all the point intensities from the noisy image based on the redundant representations of image [13]:

where is the weight assigned to point in the restoration of point and is the search area centered at the current point . According to (2), this method has the capability to reconstruct the voxel from all similar voxels in the restricted neighbor volume . Consequently, redundant information from the same MRI image and different slices can be used to reconstruct the current voxel efficiently.

The key issue of NLM is the computation of , which represents the similarity of neighborhood points. Generally speaking, within the search area , the weight is related to the distance , with and being neighborhoods around and as follows [13]:

where is a normalization constant with and is a filtering parameter.

However, the basic NLM has a great influence for computational efficiency, especially for 3D MR images. Consequently, [12] has implemented the 3D blockwise version of NLM (NLM3D), which divides the volume into overlapping blocks and treats each block as a point to perform NLM-like restoration. For NLM3D, the and in (3) are 3D patches around points and , and is the similarity of the neighbor volumes between and . Based on the constrution of 3D neighbor blocks, the NLM3D will restore noisy data from neighbor blocks within intraslices and interslices simultaneously which is helpful to preserve the structural information of different slices for MRI. Therefore, compared with classical NLM, the method has not only reduced the complexity of NLM significantly but also achieved superior denoising performance. The NLM3D has already applied ito super-revolution reconstruction of MRI [17] with different preselection step and computation formulation of .

2.2. Multidimensional Principal Component Analysis

Although PCA has been applied in image denoising widely, most denoising algorithms based on PCA assume that data lie on vector space and usually process the vectorization operation to make image into a vector. The vectorization destroys the structural information about the neighborhood.

Instead of data in vector space, any multidimensional data can be considered as tensor data in multidimensional space. Each tensor data will be treated by tensor decomposition [18]. It is helpful to preserve the structural information and enhance the computational efficiency. At present, a large number of tensor algorithms have been presented and have a wide application in computer vision, pattern recognition, and machine learning [19].

Based on PCA on vector space, [20] has developed multidimensional principal component analysis (MPCA) for tensor data and has achieved outstanding performance. For MPCA, a high-dimensional image dataset can be expressed as a tensor dataset , where is a dimensional tensor and is the number of samples in the dataset.

In tensor algebra, any tensor data can be expressed based on Tucker decomposition model as follows [20]:

where is an orthogonal matrix. So we can get

The target of MPCA is to compute orthogonal projective matrices to maximize the total scatter tensor of the projected low-dimensional feature as follows:

where is the mean of tensor data and .

Due to the difficulty in the computation of orthogonal projective matrices simultaneously, these orthogonal projective matrices can be solved iteratively. Generally speaking, it is assumed that the projective matrices are known; then we can solve the following optimized problem to obtain :

where and is the mean of tensor data . is the mode- unfolding matrix of tensor . The paper [19] has proved the vector-based and 2DPCA can be considered as the special cases of MPCA.

2.3. Proposed Method

To the best of our knowledge, this is the first attempt to introduce MPCA to MRI image denoising. It should be noticed that image denoising based on MPCA is different from its application in machine learning.

For MRI denoising, 3D MRI is a 3rd-order tensor , and are the height and width of each MRI slice, respectively, and is the number of slices, so 3D MRI can be considered as an image set . The principal components can be computed by MPCA:

where , . Generally speaking, the first principal components conserve most information of images. It is desirable to abandon smaller principal components to remove noise. There are various approaches to determine the value of . The paper specifies a constant to represent the number of the largest principal components corresponding to the largest eigenvalues. After that, the restored image is expressed:

where , , and is the restored image. Consequently, the proposed algorithm can be summarized.

(1)NLM 3D filtering is applied to denoise and obtain the initial 3D images.(2)Then the initial 3D images are processed by MPCA to remove noise furtherly.

3. Experiments

Several experiments were conducted to compare the proposed methods with related state-of-the-art methods.

In order to illustrate the performance of the proposed method, several experiments were conducted to compare the proposed methods with related state-of-the-art methods, including the 3D version of anisotropic diffusion filtering (3D-ADF) [3, 5], NLM3D [12], and OMNLM_LAPCA [16] on synthetic data and real clinical data. All experiments are performed on MATLAB R2015a.

There are some free parameters that need to be set to obtain optimal performance. For 3D anisotropic diffusion filtering, the integration constant is the maximum value, the number of iterations is 4, and the gradient modulus threshold is 70. For NLM3D and OMNLM_LAPCA, the radius of the search area is 5 and the radius of similarity area is 2. For the proposed method, the number of preserved largest principal components is 140 (see below).

Three kinds of quality measurement are used to evaluate the denoising performance. The first one is the signal-to-noise ratio (SNR), the second one is the peak signal-to-noise ratio (PSNR), and the last one is the structural similarity index (SSIM) [21].

The SNR is computed as follows:

where is the original image, is the mean of image , and is the denoised image.

The PSNR is based on the root mean square error (RMSE) between the denoised image and original image:

The SSIM is defined as follows:

where , , is the dynamic range, , and ; and are the mean of images and , respectively; and are the standard noise variance of images and , respectively; is the covariance of and .

3.1. Synthetic Data

In this part, the 3D T1-weighted MRI image in the well-known BrainWeb [22] dataset is used to evaluate the performance of the proposed method. The size of the dataset is with 1 mm³ voxel resolution. To simulate Rician noise, zero mean Gaussian noise with 3–5% standard deviation is added to the real and imaginary parts of the 3D MRI images, as shown in Figure 1.

(a)

(b)

(c)

(d)

Figure 1: Example images of the BrainWeb database. (a) Noisy-free T1-w image, (b) noisy image corrupted with a Rician noise at 3%, (c) noisy image corrupted with a Rician noise at 4%, and (d) noisy image corrupted with a Rician noise at 5%.

The denoising performance of different methods with different noise levels is compared based on SNR, PSNR, and SSIM, as shown in Table 1. It is obvious that the proposed method is superior to the other three methods under the three evaluation measurements. The denoising images and corresponding residuals are shown in Figure 2. It is consistent with the measurement that the proposed method has the best visual effect.

Table 1: The comparison of denoising methods with different noise levels on T1 weighted MRI Images.

Denoising method
3% Rician noise
4% Rician noise
5% Rician noise

SNR
PSNR
SSIM
SNR
PSNR
SSIM
SNR
PSNR
SSIM

3D-ADF
13.7207
25.4512
0.7642
12.9584
25.4144
0.7624
12.1337
25.3515
0.7604

OMNLM_LAPCA
18.5248
35.0833
0.9495
16.0899
33.2808
0.9234
14.1684
31.7119
0.8927

NLM3D
25.9919
37.5705
0.976
24.0109
35.9133
0.965
20.8779
34.5518
0.9472

Proposed
27.0346
38.5828
0.9864
25.0351
36.924
0.9755
22.9218
35.5598
0.9555

Figure 2: The noise-free image, noisy image, denoised images, residuals, and detail denoised images of different methods with 5% Rician noise for different slices.

For the proposed method, it is required to determine the optimal number of the largest principal components. To study its influence on the denoising performance, the SNR, PSNR, and SSIM with different principal components with different noise levels are shown in Figures 3–5. It can be seen that the optimal number is between 140 and 160 when the better measurements are obtained. So it is unreasonable to set a constant to choose the principal components. However, it is still an open problem in machine learning.

Figure 3: Plots of SNR versus the number of principal components for different noise levels.

Figure 4: Plots of PSNR versus the number of principal components for different noise levels.

Figure 5: Plots of SSIM versus the number of principal components for different noise levels.

Based on the idea proposed in [16], the paper has presented a two-step MRI denoising algorithm, which employed NLM3D to restore 3D MRI followed by MPCA. To validate the proposed two-step method in the paper, we also have evaluated the order of MPCA and NLM3D. Figure 6 has shown the denoised result with different order of MPCA and NLM3D. It can be seen that if we apply MPCA and then apply NLM3D, the detail information is lost and the denoised image is blurred. In contrast, the proposed steps in the paper have the capability to preserve the detail. The possible reason is that both of noise and detail are high frequency signals; the threshold technique of MPCA will remove the detail information while removing the noise if we conduct MPCA on the noisy image directly. Consequently, NLM3D is unable to restore the detail from the blurred neighboring block.

(a) MPCA + NLM3D

(b) NLM3D + MPCA

Figure 6: The denoised results with different order of MPCA and NLM3D (5% Rician noise). (a) First apply MPCA and then apply NLM3D but in (b) first apply NLM3D and then apply MPCA.

All denoising methods were performed in MATLAB 2015 on a Windows 7 computer equipped with an Intel Core i7-5600U, 2.6 GHz CPU and 8 GB RAM. To denoise typical 3D dataset with the size of , the corresponding computational time is listed in Table 2.

Table 2: The comparison of computational time.

Denoising methods
The computational time(s)

3D-ADF
17.6905

OMNLM_LAPCA
32.8070

NLM3D
110.7295

Proposed
111.0259

It cannot be denied that, compared with other algorithms, the proposed method will spend more time to denoise 3D MRI since it makes use of the 3D structure information from neighboring slices. It is also believed that the implementation of the proposed methods using MATLAB/C MEX techniques and parallel computations on graphic processing units may significantly further accelerate the filtering.

3.2. Validation on Real Clinical Data

To evaluate the proposed method on real clinical data, the experiments are conducted on real T1-w MRI data. The data were acquired on a GE MR750 3.0T scanner. The anatomical images were scanned using a T1-weighted axial sequence parallel to the anterior-commissure-posterior-commissure line. Each anatomical scan has 156 axial slices (spatial resolution = 1 mm × 1 mm × 1 mm, field of view = 256 mm × 256 mm, time repetition (TR) = 8.124 ms). The noisy-free image and noisy image are shown in Figure 7.

Figure 7: The example images of real clinical data and noisy data with 5% Rician noise.

For real clinical data, the denoised results are shown in Figure 8 and Table 3. It can be seen that the proposed method also has achieved the best visual result. It may be that NLM3D with block representation restores noisy data from different neighboring slices. Moreover, in contrast with the principal components of PCA in vector space, MPCA seeks the principal components in tensor space, which has the capability to preserve the structure information of neighboring voxels and slices. At the same time, most image information focuses on the first principal components, so the threshold technique is helpful to remove noise further.

Table 3: The comparison of different methods with 5% Rician noise.

Denoising method
SNR
PSNR
SSIM

3D-ADF
9.7543
24.8338
0.6717

OMNLM_LAPCA
9.7943
30.42
0.8415

NLM3D
18.7476
32.6647
0.9034

Proposed
19.7629
33.6679
0.9236

(a)

(b)

(c)

Figure 8: The denoised images, residuals, and detail denoised images of different methods with 5% Rician noise for different slices.

4. Discussion

The paper has proposed a structure preserving MRI denoising algorithm. The method has integrated NLM3D and MPCA to restore noisy image from 3D neighborhood and has achieved a better result compared with some famous MRI denoising methods, such as 3D-ADF, OMNLM_LAPCA, and NLM3D. However, the confusing question of the proposed method is how to determine the optimal number of principal components, which will affect the denoising effect. So our next work will research the selection problem of principal components. We will consider the cumulative energy or the scoring of principal components in the future.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgments

This work is supported by the National Natural Science Funds of China (Grant no. 61502059) and the Young Scientist Project of Chengdu University (no. 2013XJZ21). The authors are glad to to have learned a lot from the reviewers and the editor who made many excellent comments to improve the paper presentation.

References

J. Mohan, V. Krishnaveni, and Y. Guo, “A survey on the magnetic resonance image denoising methods,” Biomedical Signal Processing and Control, vol. 9, no. 1, pp. 56–69, 2014.

E. R. McVeigh, R. M. Henkelman, and M. J. Bronskill, “Noise and filtration in magnetic resonance imaging,” Medical Physics, vol. 12, no. 5, pp. 586–591, 1985.

P. Perona and J. Malik, “Scale-space and edge detection using anisotropic diffusion,” IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 12, no. 7, pp. 629–639, 1990.

F. Zhang and L. Ma, “MRI denoising using the anisotropic coupled diffusion equations,” in Proceedings of the 3rd International Conference on BioMedical Engineering and Informatics (BMEI '10), pp. 397–401, October 2010.

A. A. Samsonov and C. R. Johnson, “Noise-adaptive nonlinear diffusion filtering of MR images with spatially varying noise levels,” Magnetic Resonance in Medicine, vol. 52, no. 4, pp. 798–806, 2004.

I. Delakis, O. Hammad, and R. I. Kitney, “Wavelet-based de-noising algorithm for images acquired with parallel magnetic resonance imaging (MRI),” Physics in Medicine and Biology, vol. 52, no. 13, pp. 3741–3751, 2007.

V. G. Ashamol, G. Sreelekha, and P. S. Sathidevi, “Diffusion-based image denoising combining curvelet and wavelet,” in Proceedings of the 15th International Conference on Systems, Signals and Image Processing (IWSSIP '08), pp. 169–172, June 2008.

J. Rajan, J. Veraart, J. Van Audekerke, M. Verhoye, and J. Sijbers, “Nonlocal maximum likelihood estimation method for denoising multiple-coil magnetic resonance images,” Magnetic Resonance Imaging, vol. 30, no. 10, pp. 1512–1518, 2012.

S. Aja-Fernández, C. Alberola-López, and C.-F. Westin, “Noise and signal estimation in magnitude MRI and rician distributed images: a LMMSE approach,” IEEE Transactions on Image Processing, vol. 17, no. 8, pp. 1383–1398, 2008.

S. P. Awate and R. T. Whitaker, “Feature-preserving MRI denoising: a nonparametric empirical bayes approach,” IEEE Transactions on Medical Imaging, vol. 26, no. 9, pp. 1242–1255, 2007.

A. Buades, B. Coll, and J. M. Morel, “A review of image denoising algorithms, with a new one,” Multiscale Modeling & Simulation, vol. 4, no. 2, pp. 490–530, 2005.

P. Coupe, P. Yger, S. Prima, P. Hellier, C. Kervrann, and C. Barillot, “An optimized blockwise nonlocal means denoising filter for 3-D magnetic resonance images,” IEEE Transactions on Medical Imaging, vol. 27, no. 4, pp. 425–441, 2008.

J. V. Manjón, J. Carbonell-Caballero, J. J. Lull, G. García-Martí, L. Martí-Bonmatí, and M. Robles, “MRI denoising using non-local means,” Medical Image Analysis, vol. 12, no. 4, pp. 514–523, 2008.

J. Rajan, A. J. Den Dekker, and J. Sijbers, “A new non-local maximum likelihood estimation method for Rician noise reduction in magnetic resonance images using the Kolmogorov-Smirnov test,” Signal Processing, vol. 103, pp. 16–23, 2014.

P. Coupé, J. V. Manjón, M. Robles, and D. L. Collins, “Adaptive multiresolution non-local means filter for three-dimensional magnetic resonance image denoising,” IET Image Processing, vol. 6, no. 5, pp. 558–568, 2012.

J. V. Manjn, N. A. Thacker, J. J. Lull, G. Garcia-Martí, L. Martí-Bonmatí, and M. Robles, “Multicomponent MR image denoising,” International Journal of Biomedical Imaging, vol. 2009, Article ID 756897, 18 pages, 2009.

J. V. Manjn, P. Coup, A. Buades, V. Fonov, D. Louis Collins, and M. Robles, “Non-local MRI upsampling,” Medical Image Analysis, vol. 14, no. 6, pp. 784–792, 2010.

T. G. Kolda and B. W. Bader, “Tensor decompositions and applications,” SIAM Review, vol. 51, no. 3, pp. 455–500, 2009.

C. Liu, T. Yan, W. D. Zhao et al., “Incremental tensor principal component analysis for handwritten digit recognition,” Mathematical Problems in Engineering, vol. 2014, Article ID 819758, 10 pages, 2014.

H. Lu, K. N. Plataniotis, and A. N. Venetsanopoulos, “MPCA: multilinear principal component analysis of tensor objects,” IEEE Transactions on Neural Networks, vol. 19, no. 1, pp. 18–39, 2008.

Z. Wang, A. C. Bovik, H. R. Sheikh, and E. P. Simoncelli, “Image quality assessment: from error visibility to structural similarity,” IEEE Transactions on Image Processing, vol. 13, no. 4, pp. 600–612, 2004.

C. A. Cocosco, V. Kollokian, R. K.-S. Kwan, and A. C. Evans, “Brainweb: online interface to a 3D MRI simulated brain database,” NeuroImage, vol. 5, no. 4, p. S425, 1997.