Accepted Manuscript Title: Denoising Filters Evaluation for Magnetic Resonance Images Author: Danni Ai Jian Yang Jingfan Fan Weijian Cong Xuehu Wang PII: DOI: Reference:
S0030-4026(15)00708-1 http://dx.doi.org/doi:10.1016/j.ijleo.2015.07.155 IJLEO 55871
To appear in: Received date: Accepted date:
1-7-2014 26-7-2015
Please cite this article as: D. Ai, J. Yang, J. Fan, W. Cong, X. Wang, Denoising Filters Evaluation for Magnetic Resonance Images, Optik - International Journal for Light and Electron Optics (2015), http://dx.doi.org/10.1016/j.ijleo.2015.07.155 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
Denoising Filters Evaluation for Magnetic Resonance Images Danni Ai, Jian Yang*, Jingfan Fan, Weijian Cong, Xuehu Wang Beijing Engineering Research Center of Mixed Reality and Advanced Display, School of Optics and Electronics, Beijing Institute of Technology, Beijing 100081, China
ip t
*
[email protected]
Abstraction
Eleven denoising filters, proposed during last fifteen years, are introduced and compared for magnetic
cr
resonance images. Among them, the state-of-art denoising algorithms, NLM and BM3D, have attracted
much attention. Several expansions are proposed to improve the noise reduction based on these two algorithms. On the other hand, optimal dictionaries, sparse representations and appropriate shapes of
us
the transform’s support are also considered for the image denoising. Based on the estimated noise variance, the comparison of various filters is implemented by measuring the signal-noise-ratio (SNR), resolution and uniformity of a phantom image. The subjective judgment of denoising effectiveness is
an
executed for a clinical image. And the computational time is finally evaluated. Keywords: Denoising filters; Magnetic Resonance Image; Noise variance 1. Introduction
M
The noise is an inevitable problem derived in the formation procedure of magnetic resonance imaging (MRI). The observed MR image y is expressed as:
y=x+n
(1)
d
where x is an ideal image in the presence of noise n. Image denoising is to design an algorithm that can remove or attenuate the noise from the observation. It is an accepted preprocessing step in many
te
magnetic resonance (MR) image processing and analysis tasks, such as classification, segmentation or registration. It also intends to improve the quality of image for aiding visual inspection. Because of the uncertainty of the source and distribution of noise, MRI denoising is still a challenging problem in the
Ac ce p
field of medical image processing.
Choice of denoising algorithm is application dependent and depends on the type of noise present. Various types of denoising algorithms are proposed in the last decade or so for both natural image and MRI. The underlying characteristics of the given data are learned for image denoising. In this paper, we make a selection of the typical denoising methods, which effectively reduce the noise and enhance the image quality. They can be classified into two categories: the spatial methods [1-4] and the spatial-transformed methods [5-11]. One of the state-of-art denoising filter in the spatial space is non-local means (NLM) [1], which is expanded into Unbiased Non-local Means (UNLM) [2] and Foveated non-local means (FOVNLM) [3]. Bilateral filter is also a typical spatial method [4]. The most attractive spatial-transformed method is called block-matching and 3D filtering (BM3D) [5] that has several expansions to improve the noise reduction, such as BM3D_SHARP [6] and BM3D_SAPCA [7]. Additionally, optimal dictionaries, sparse representations [9,10,11] and appropriate shapes of the transform’s support [8] are also considered for the image denoising. The estimation of the noise level in an MRI is a crucial parameter for the denoising procedure. A current research, conducted by A. Foi [12], has contributed to the noise estimation. In this paper, we do the comparison based on the estimated noise level. The objective comparison of various filters is implemented by measuring the signal-noise-ratio (SNR), resolution and uniformity of a phantom image. The subjective judgment of
Page 1 of 12
denoising effectiveness is executed for a clinical image. And the computational time is finally evaluated. The organization of the paper is as follows: A brief introduction of eleven denoising filters is given in Section 2, which is followed by a short presentation of noise estimation. The comparison experimental results are discussed in Section 4, and finally, Section 5 concludes the work. Multiple Filters
ip t
2.
We give a brief description of eleven denosing filters, for which reliable faithful implementations are
2.1.
cr
available and will be compared in section 4. Spatial Methods
For spatial methods, the photometric similarity is directly compared for denoising in the image space
2.1.1.
us
rather than in the transformation space. Non-local means (NLM) filter
an
For every small patch in the MR image, one can find several similar patches in the same image. The NLM algorithm tries to make use of the redundancy of the similar patches [1]. Given a noisy MR image I, the filtered value at a point i within the reference patch Ni is computed as a weighted average
where the weights
M
of all the similar patches Nj in the image
(2)
depend on the similarity between the neighborhoods Ni and Nj of pixels i
d
and j. The similarity between two pixels i and j will depend on the similarity of the intensity grey-level
te
vectors I( ) and I( ). The pixels with a similar grey-level neighborhood to I( ) will have larger is then calculated as
Ac ce p
weights on the average. The similarity
where d
and
(3)
, the parameter h controls the decay of the
exponential function and k is the size of patches. 2.1.2.
Unbiased Non-local Means (UNLM) filter
In low signal-to-noise image, the square of MRI magnitude image is used for filtering. The noise bias is equal to
. Thus the UNLM will be defined as (4)
where
, μ is the mean value of the background of the squared magnitude image [2].
Page 2 of 12
2.1.3.
Foveated non-local means (FOVNLM) filter
Inspired by the human visual system, a patch foveation operator and measure patch similarity are introduced through a foveated distance [3]. By considering nonlocal means algorithm, an construction of the corresponding foveation operator is derived, that is called Foveated non-local means algorithm. 2.1.4.
Bilateral filter
ip t
Bilateral filter extends the concept of Gaussian smoothing by simultaneously considering geometric
closeness and photometric [4]. Pixels that are very different in intensity from the reference pixel are applied, one in the spatial domain, and the other one in the intensity domain. 2.2.
Spatial-Transformed Methods
cr
weighted less even though they may be in close proximity to the central pixel. Two Gaussian filters are
2.2.1.
us
After transforming image space into other domain, the image denosing is implemented in this domain. BM3D
an
By extending DCT and NLM, block-matching and 3D filtering (BM3D) is an efficient sliding window denoising method [5]. There are five steps in BM3D: (1) the similar image patches to a given patch are found and grouped into 3D block; (2) the 3D blocks are converted into a transformed subspace; (3) The transform spectrum coefficients are shrunk; (4) 3D transformation is inversed into original space; (5)
M
an adaptive aggregation procedure is applied to the filtered patches that are returned to the initial positions. BM3D_SHARP
d
2.2.2.
In order to simultaneously sharpen image details and attenuate noise, BM3D_SHARPENING
te
combines the BM3D with a transform-domain sharpening technique [6]. It is an expansion of BM3D. Following the shrinkage step of BM3D, the alpha-root is applied on the 2D or 3D transform spectra in
Ac ce p
order to sharpen image details. 2.2.3.
BM3D_SAPCA
As another expansion of BM3D, BM3D_SAPCA exploits nonlocal image modeling, principal component analysis (PCA), and local shape-adaptive (SA) anisotropic estimation [7]. Instead of using square image patches and a fixed 3D transform in BM3D, BM3D_SAPCA groups mutually similar adaptive-shape neighborhoods, and does the transform with adaptive learning basis. 2.2.4.
Pointwise Shape-Adaptive DCT (SA-DCT) filter
Pointwise shape-adaptive DCT (SA-DCT) makes use of the anisotropic LPA-ICI (local polynomial approximation - intersection of confidence intervals) technique to obtain an appropriate shape of the transform’s support which is not a square [8]. The SA-DCT is applied on this shape. A local estimate within the adaptive-shape region is reconstructed by the thresholded SA-DCT coefficients. Because of the overlapping of regions, the denoised image is obtained by averaging the local estimated with adaptive weights. 2.2.5.
K-SVD
K-SVD intends to find the best dictionary to represent the input data as sparse compositions as possible
Page 3 of 12
[9]. It is an iterative method that alternates between sparse coding of the examples based on the current dictionary and a process of updating the dictionary atoms to better fit the data. 2.2.6.
CSR
Clustering-based sparse representation (CSR) treats the local and nonlocal as peers, which simultaneously considers the local and nonlocal under a unified theoretic framework [10]. An iterative
ip t
shrinkage solution is developed for the double-header l1 optimization problem, which understands sparsity by unifying dictionary learning and structural clustering into a variational framework. LPG-PCA
cr
2.2.7.
To preserve the image local structures, principal component analysis with local pixel grouping (LPG-PCA) [11] models a block of a target pixel and its nearest neighbors as a vector variable. In a
us
local window centered on the target pixel, the similar blocks are selected by using block matching based local pixel grouping. As the training samples for PCA learning, all the similar blocks are
3.
an
unfolded as vectors to calculate the eigenvector matrix for target transformation. Noise Estimation
The noise variance is an important measure of MRI denoising. Additionally, it is also a crucial parameter in most of the denoising algorithms mentioned above. Usually, the noise in a single MRI can
M
be well modeled by a Rician distribution. A. Foi [12] proposed optimal forward and inverse variance-stabilizing transformations for the rice distribution. And a stable and fast iterative procedure for estimating the noise level is proposed from a single Rician-distributed image. This estimation
d
method is called as the automatic method. In this paper, we use two noise estimation methods to obtain the denoising results for the phantom data:
te
the manual adjustment of the noise values and the automatic method [12]. While for the clinical data,
Ac ce p
only the automatic method is utilized.
Fig. 1 Phantom image 4.
Experimental Results
We applied eleven denoising filters to both the phantom data and the clinical data. All the denoising filters can be tested with minor modification for fitting the input images. The original source codes are referred to the websites: NLM [13], UNLM [13], FOVNLM [14], Bilateral [15], BM3D [16], BM3DSHARP [16], BM3D_SAPCA [16], SA_DCT [17], KSVD [18], CSR [19], LPGPCA [19]. The experiments were carried out on a standard PC (Intel(R) Core (TM) i5-2400 CPU 3.10GHZ and2.00
Page 4 of 12
GB RAM) running Windows 7. 4.1.
Phantom data
4.1.1.
Evaluation methods
For the phantom data shown in Fig. 1, the signal-noise-ratio (SNR), resolution and uniformity are used
ip t
for evaluating the denoising effectiveness with different filters. The phantom data size is 256×256 with intensity [0, 255]. (1) SNR
SNR require manual selection of the region of interest (ROI), which is shown in Fig. 2. In the
us
is the mean value of selected
an
is the mean value of selected signal in the ROI;
background;
(5)
is the standard deviation of selected signal.
d
signal ROI (for , )
M
where
cr
experiments, SNR is calculated by:
Ac ce p
te
background ROI (for )
Fig. 2 ROI selection for the SNR
(2) Uniformity
Uniformity was calculated as the difference between the maximum and the minimum mean values of the ROI. It also requires manual selection of ROI shown in Fig. 3. Uniformity is calculated by the following equation,
where
and
(6)
are the maximum mean and the minimum mean values of the selected
ROI, respectively.
Page 5 of 12
ROI ,
)
Fig. 3 ROI selection for the uniformity
(3) Resolution
cr
ip t
(for
us
The spatial resolution quantifies the distiguishment of the different features for an image. We calculate the modulation depth with the signal’s peak and valley values on a bar-shaped phantom. In Fig. 4, the pixels across the bar from (88, 82) to (102, 68), shown in yellow line, are chosen to
an
calculate the resolution. The resolution here is defined as following:
(7)
Ac ce p
te
d
M
where pi are the local maximum values and vi are the local minimum values, i=1,2,3.
Fig. 4 Selected pixels for resolution calculation
4.1.2.
Comparison results
For the phantom data, two methods are used for determining the noise variance: automatic method and manual method.
(1) Automatic method
Based on the method proposed in [12], we automatically calculate the noise variance as 5.502 for the phantom data shown in Fig. 1. This parameter is used for all the filters. After measuring the denoised images, we obtain the comparison results for SNR, uniformity and resolution, which are shown in Fig. 5. It is clear in Fig. 5 that the uniformity is similar for each filter, while the SNR and resolution change a lot. Among all the filers, UNLM filter can obtain the highest SNR and the second higher resolution. BM3DSHARP produces the highest resolution, but the SNR is the worst. As shown in Fig. 6 that displays denoised phantom images, the image filtered by BM3DSHARP is getting
Page 6 of 12
worse than the original image. The parameter may not be the best one for BM3DSHARP, because it is not just a denoising filter, but also sharpen image details simultaneously [6]. Uniformity
SNR 25
1
20
0.8
15
ip t
0.6
10
0.4
5
0.2
0
cr
0
Type of filters
us
Left Range for Resolution & Uniformity Right Range for SNR (red)
Resolution 1.2
Ac ce p
te
d
M
an
Fig. 5 Comparison of resolution, uniformity and SNR with parameter 5.502
Page 7 of 12
ip t UNLM
an
us
cr
NLM
Original
Bilateral
BM3D
BM3DSHARP
BM3DSAPCA
SADCT
KSVD
CSR
LPGPCA
Ac ce p
te
d
M
FOVNLM
Fig. 6 Denoised phantom images with different filters and the original image (2) Manual method In order to find the best parameter (noise variance) for the phantom image, we manually change the value of parameter in the range [1, 40]. According to Fig. 5, the uniformity changes slightly,
Page 8 of 12
which will not be as an evaluation method for the manual measurement. The best parameters shown in table 1 are obtained if both SNR and resolution of the denoised images are more than these of the original image. It is clear from Table 2 that the best parameter with the automatic method is inside or near the best parameter range with manual method for most filters, except for BM3DSHARP. It further explains the reason for the worse denoising result using BM3DSHARP with automatic parameter. Figure 7
ip t
shows the denoising result using BM3DSHARP with manual parameter 21, which is much better than that in Fig. 6. Table 2. The best parameters of each filters BM3D_SAPCA
BM3DSHARP
FOVNLM
Best Parameters
2-4
3
21
2-4
Filters
UNLM
Bilateral
KSVD
CSR
Best Parameters
1-9
3-4
2-3
2
SA_DCT
NLM
cr
BM3D
2-5
2-6
LPGPCA 3-4
d
M
an
us
Filters
te
Fig. 7 The denoising result using BM3DSHARP with manual parameter 21 Figure 8 and 9 show the SNR and resolution with the best range of parameters, respectively.
Ac ce p
“High” and “Low” assign the SNR and resolution range of each filter. “Threshold” denotes these of the original image. Considering both Fig. 8 and Fig. 9, UNLM is superior to other filters. Here, UNLM is proposed for MR image and treats the noise with Rician distribution. Also, NLM and BM3DSHARP are better filters for MRI denoising. SNR
23.5
23
Range
22.5
22
21.5
High
21
Low
20.5
Threshold
20
Type of filters
Fig. 8 SNR with the best parameters
Page 9 of 12
Resolution 1.05
Range
1 0.95 0.9
High Low
ip t
0.85
Threshold
cr
0.8
us
Type of filters
4.2.
an
Fig. 9 Resolution with the best parameters
Clinical Data
We compare eleven denoising filters by using a sT2 MR image of 512×512 pixels acquired on 1.5T from OsiriX dataset [20]. The noise variance of the clinical image is 3.8, which is automatically
M
calculated for the image with method [12]. All the filters utilize this parameter, except for the BM3DSHARP using value 12 as the parameter. The denoised images as well as the original image are shown in Fig. 11. According to the opinions of the clinically trained people, the BM3D_SAPCA is the
d
best one among all the filters, which can preserve the image structures while reducing noise. Additionally, BM3D_SHARP and UNLM also present competitive denoising solution compared with
te
the rest filters. The computational time by using all the filters is illustrated in Fig. 10. Obviously, BM3D_SHARP is the fastest denoising filter.
Ac ce p
900
823
800
747.33
700
Time (s)
600 500 400 300 200 100
0
154.66
115.79
43.39
183.89
151.19 6.26
4.03
0.95
2.39
Type of filters
Fig. 10 The comparison of the computational time
Page 10 of 12
ip t
NLM
UNLM
an
us
cr
Original
Bilateral
BM3D
Ac ce p
te
d
M
FOVNLM
BM3DSHARP
BM3D_SAPCA
SADCT
KSVD
CSR
LPGPCA
Fig. 11 The original clinical image and the denoised results with eleven filters 5.
Conclusion
This paper introduced eleven typical image filters for MR image denosing. To find the noise variance, we used manual and automatic methods for the phantom image, and proved effective denoising results
Page 11 of 12
are achieved on removing the noise as well as keeping the useful information with both methods. Based on this conclusion, the automatic method is directly utilized for the clinical image. The comparative study of various denoising filters for MRI shows that the UNLM is the best with objective judgment for the phantom image, while BM3D_SAPCA is superior to other filters with subjective judgment for the clinical image.
ip t
Conflict of Interests
The authors declare that they have no competing interests.
cr
Acknowledgments
This work was supported by the National Basic Research Program of China (2010CB732505, 2013CB328806), the Key Projects in the National Science & Technology Pillar Program
us
(2013BAI01B01) and the National Hi-Tech Research and Development Program (2013AA013703). References
Ac ce p
te
d
M
an
1. Buades, A., Coll, B., Morel, J.M., 2005.A review of image denoising algorithms, with a new one. Multiscale Modeling and Simulation 4 (2), 490–530, 2005. 2. J. Manj´on, J. Carbonell-Caballero, J. Lull, G. Garc´ıa-Mart´ı, L. Mart´ı-Bonmat´ı, and M. Robles, MRI denoising using non-local means, Med. Image Anal., 12, pp. 514–523, 2008 3. A. Foi and G. Boracchi, “Foveated self-similarity in nonlocal image filtering”, Proc. SPIE Electronic Imaging 2012, Human Vision and Electronic Imaging, 8291-32, Burlingame (CA), USA, January 2012. 4. C. Tomasi and R. Manduchi, “Bilateral filtering for gray and color images,” in Proceedings of the IEEE International Conference on Computer Vision, pp. 839–846, 1998. 5. K. Dabov, A. Foi, V. Katkovnik, and K. Egiazarian, "Image denoising by sparse 3D transform-domain collaborative filtering," IEEE Transactions on Image Processing, 16(8), 2007. 6. K. Dabov, A. Foi, V. Katkovnik, and K. Egiazarian, "Joint image sharpening and denoising by 3D transform-domain collaborative filtering," Proc. 2007 Int. TICSP Workshop Spectral Meth.Multirate Signal Process, SMMSP 2007. 7. K. Dabov, A. Foi, V. Katkovnik, and K. Egiazarian, "BM3D Image Denoising with Shape - Adaptive Principal Component Analysis", Proc. Workshop on Signal Processing with Adaptive Sparse Structured Representations (SPARS'09), Saint-Malo, France, April 2009. 8. Foi, A., V. Katkovnik, and K. Egiazarian, “Pointwise Shape-Adaptive DCT as an overcomplete denoising tool”, Proc. Int. TICSP Workshop Spectral Meth. Multirate Signal Process., SMMSP 2005, Riga, 2005. 9. M. Aharon, M. Elad, and A.M. Bruckstein, “The K-SVD: An Algorithm for Designing of Overcomplete Dictionaries for Sparse Representation”, IEEE Trans. On Signal Processing, 54(11), pp. 4311-4322, 2006. 10. W. Dong, X. Li, L. Zhang, G. Shi, “Sparsity-based Image Denoising via Dictionary Learning and Structural Clustering.” IEEE Conference on Computer Vision and Pattern Recognition (CVPR), pp: 457-464, 2011. 11. L.Zhang, W.Dong, D.Zhang, G.Shi, “Two-stage image denoising by principal component analysis with local pixel grouping.” Pattern Recognition vol. 43, pp: 1531-1549, 2010. 12. A. Foi, “Noise estimation and removal in MR imaging: the variance stabilization approach,” in Proceedings of the IEEE International Symposium on Biomedical Imaging: From Nano to Macro, Chicago, IL, USA, 2011. 13. http://personales.upv.es/jmanjon/denoising/nlm2d.htm 14. http://www.cs.tut.fi/~foi/FoveatedNL/#ref_software 15. http://www.mathworks.com/matlabcentral/fileexchange/authors/21859 16. http://www.cs.tut.fi/~foi/GCF-BM3D/ 17. http://www.cs.tut.fi/~foi/SA-DCT/ 18. http://www.cs.technion.ac.il/~elad/software/ 19. http://www4.comp.polyu.edu.hk/~cslzhang/ 20. http://www.osirix-viewer.com/datasets/
Page 12 of 12