A 3D MRI denoising algorithm based on Bayesian theory
 Fabio Baselice^{1}Email authorView ORCID ID profile,
 Giampaolo Ferraioli^{2} and
 Vito Pascazio^{1}
DOI: 10.1186/s129380170319x
© The Author(s) 2017
Received: 24 October 2016
Accepted: 30 January 2017
Published: 7 February 2017
Abstract
Background
Within this manuscript a noise filtering technique for magnetic resonance image stack is presented. Magnetic resonance images are usually affected by artifacts and noise due to several reasons. Several denoising approaches have been proposed in literature, with different tradeoff between computational complexity, regularization and noise reduction. Most of them is supervised, i.e. requires the set up of several parameters. A completely unsupervised approach could have a positive impact on the community.
Results
The method exploits Markov random fields in order to implement a 3D maximum a posteriori estimator of the image. Due to the local nature of the considered model, the algorithm is able do adapt the smoothing intensity to the local characteristics of the images by analyzing the 3D neighborhood of each voxel. The effect is a combination of details preservation and noise reduction. The algorithm has been compared to other widely adopted denoising methodologies in MRI. Both simulated and real datasets have been considered for validation. Real datasets have been acquired at 1.5 and 3 T. The methodology is able to provide interesting results both in terms of noise reduction and edge preservation without any supervision.
Conclusions
A novel method for regularizing 3D MR image stacks is presented. The approach exploits Markov random fields for locally adapt filter intensity. Compared to other widely adopted noise filters, the method has provided interesting results without requiring the tuning of any parameter by the user.
Keywords
3D MRI denoising Maximum a posteriori Markov random fields Statistical signal processingBackground
As most of clinical systems, Magnetic Resonance images are usually affected by artifacts and noise due to several reasons, such as the electronic noise generated by the thermal agitation of the charge carriers (thermal noise) and the imaging scanner technical limitations. Depending on the region of interest and on the application, noise could severely degrade the quality of acquired MR images and produce further errors in quantitative assessments from the data. Therefore, techniques for reducing the amount of noise affecting the acquisitions are usually implemented. In general, denoising methods aim at removing the undesirable noise while preserving the image details and local geometries [1].
One of the most used approaches for image denoising exploits the principle of averaging similar pixels of the image in order to lower the noise variance in the output image. This is mainly implemented in linear filters, such as Gaussian and boxcar. The main drawback is that the spatial average is equivalent to a lowpass filter in the frequency domain, thus image details and edges are degraded after the filter application. In other words, a blurring effect appears on the image [2].
To improve detail and edge preservation, nonlinear models have been proposed in literature [3, 4]. Anisotropic nonlinear diffusion belongs to this family of filters and has been developed for both regularizing the image and preserving edges [5]. The main drawback of this technique is the number of parameters that have to be tuned in order to reach effective performances.
Denoising can also be performed in transformed domains, where the separation between image and noise is expected to be easier. Of course, a proper transformation to the acquired data has to be applied. Among all, Wavelets are often adopted [6]. In case of MRI image, such approach introduces a bias in the filtered image. In order to reduce this disadvantage, the squared values of the image can be considered as the initial noisy one. However, such approach hardly preserve fine details of the images, especially in case of low signal to noise ratio (SNR). By considering as the transformed domain the spatial frequency one, Wiener filter has also to be taken into account [7, 8].
A different denoising approach consists in estimating the noise statistical parameters from the acquired images in order to reconstruct the noisefree image in a more effective way. To estimate noise characteristics, maximum likelihood estimation (MLE) can be implemented [9]. The main drawback of this approach is the assumption of constant signal in small regions, which produces poor details preservation.
Recently, a new filtering approach has been proposed, with very interesting detail preservation performances [2]. The methods assumes that there is redundancy across the image, and that similar patches can be found and jointly exploited for reducing noise. Such kind of methods, based on non local mean approaches [10], have proven to be very effective in images with high redundancy, but in some cases such as complicated structures or partial volume effect it fails, resulting in detail loss [11]. Another drawback of the technique is the required high computational burden.
In addition to previously reported ones, Markov random field (MRF) based methods provide interesting results [12, 13]. MRF exploits the spatial correlation information between a pixel and its neighborhood. This helps in reaching an estimation robust against noise and at the same time able to preserve fine structures and edges [14]. In other words, MRF allows a spatially adaptive noise regularization, able to reduce or increment signal smoothing according to the local statistical characteristics of the image. Often the implementation is iterative, making the approach computationally heavy.
Within this manuscript, a denoising technique based on MRF is exploited. In particular, following the approach of [14] developed for diffusion tensor MRI (DTMRI), a maximum a posteriori (MAP) estimator is proposed for regularizing 3D amplitude MRI acquisition stacks. The peculiarity of the approach consists in defining a 3D local Gaussian MRF (LGMRF) that effectively adapts the model to the local behavior of the unknown image. In particular, with respect to a classic GMRF, this model considers a hyperparameters map that describes the spatial correlation between each pixel and its neighborhood. Such characteristic allows tuning the filter intensity, i.e. regularizing smooth areas while preserving edges and small details in an unsupervised way. In particular, if fine details are found a weak smoothing is applied in order to preserve them, while in case of flat areas a stronger regularization is implemented. The smoothing effect is automatically tuned by the MRF model in order to find the optimal tradeoff between noise reduction and details preservation.
The manuscript is organized as follows: in the next Section the proposed methodology is presented, while in the following Section the framework for testing the performances of the approach is reported. Within the “Conclusion” section, the results are reported both in case of simulated and real datasets, together with a discussion about achievable performances compared to other widely adopted denoising methodologies. Finally, conclusions are drawn.
Methods
The LGMRF model of Eq. (7) measures the differences between the value of each pixels and its surroundings, weighting them via the hyperparameters.
From literature [20], the model is called “local” as the hyperparameters \(\theta _{k,q}\) used for tuning the MRF are locally defined. A simple Gaussian MRF adopts a scalar hyperparameter (one value for the whole image [21]), while the LGRMF defines for each voxel multiple hyperparameters values (one for each neighboring location). The values \(\theta _{k,q}\) allow to differently weight the neighboring pixels in the energy function. In particular, a low value of \(\theta _{k,q}\) forces strong regularization since it is an index of high spatial correlation in the considered region (i.e. between \(b_k\) and \(b_q\)). On the contrary, a high value of \(\theta _{k,q}\) means low spatial correlation between pixels k and q, i.e. an edge or a small detail is more probable, lowering the filtering intensity. Clearly, the hyperparameters are not known and have to be estimated from the available data.
Simulation and experiments
Specifications of filters used for comparison
Method  Experiment 

Anisotropic diffusion  \(n_{Iter}=4,~ \Delta _t=3/44,~k=70,~ c= 1/[1 + ( \ \nabla _I k)^2]\) 
3D bilateral  \(\sigma _{x,y}=5,~ \sigma _z=5,~ \sigma _R=15,~ sam_S=5,~ sam_R=15\) 
LMMSE  Window size = \(7 \times 7\) 
BM3D  \(\sigma _{noise}\) has been provided 
All the considered algorithms have been implemented in Mathworks™Matlab^{®} environment by using a Dell™Optiplex 990 workstation with Linux Debian as operative system.
First, a simulated case study is implemented in order to quantitatively evaluate noise removal effectiveness. In particular, mean square error (MSE) and Structural Similarity Index (SSIM) [28] are considered as quality indexes.
The simulated case study exploits Matlab^{®} 3D MRI head phantom, which is composed of 27 slices of \(128 \times 128\) voxels. Rice distributed noise has been considered with different scale parameter \(\sigma\) in order to evaluate performances under different levels of noise. In particular, the adopted values were between 0.5 and 8, corresponding to an SNR between 9 and 33 dB.
3 T real dataset: imaging protocol details
MRI scanner  Philips achieva 
Field intensity  3.0T 
Sequence  Spin Echo 
FOV  230 \(\times\) 230 \(\times\) 135 mm 
Voxel size  0.45 \(\times\) 0.45 \(\times\) 4.5 mm 
Stack resolution  512 \(\times\) 512 \(\times\) 30 pixels 
The second real dataset also refers to an MRI head axial acquisition, but an 1.5 T scanner was considered. This dataset is public available on the dicom.nema.org website. The acquisition parameters are reported in Table 3.
1.5 T real axial dataset: imaging protocol details
Field intensity  1.5 T 
Sequence  Fast Spin Echo 
Modality  T2 weighted 
FOV  220 \(\times\) 220 \(\times\) 128 mm 
Voxel size  0.43 \(\times\) 0.43 \(\times\) 2.0 mm 
Stack resolution  512 \(\times\) 512 \(\times\) 64 pixels 
1.5 T real sagittal dataset: imaging protocol details
Field intensity  1.5 T 
Sequence  Fast Spin Echo 
Modality  T2 weighted 
FOV  281 \(\times\) 281 \(\times\) 45 mm 
Voxel size  0.55 \(\times\) 0.55 \(\times\) 5.0 mm 
Stack resolution  512 \(\times\) 512 \(\times\) 9 pixels 
Results and discussion
By looking at the filtering results in case of the simulated case study, reported in Fig. 3, qualitative information on the algorithm filtering effectiveness can be inferred. The anisotropic diffusion filter produces the smoothest results, with very low performances in terms of edges and details retrieval. Also the LMMSE filtered image is blurred, with constant areas and very few small anatomical structures visible. On the other hand, 3D bilateral and BM3D approaches are capable of better details preservation, but the noise reduction effectiveness is limited. Proposed MAP approach can be placed in the middle. The preservation accuracy of edges and small elements is globally similar to BM3D filter, but the noise reduction effectiveness is higher.
Both SSIM and MSE graphs computed for different noise levels and reported in Fig. 5 confirm such findings. In particular, the MSE curves show a higher effectiveness of proposed MAP approach in reducing noise, with a small advantage of BM3D only in case of very low noise (\(\sigma \le 1.5\), SNR \(\ge\) 23 dB). We recall that the proposed MAP methodology has been developed under the assumption of low SNR, so an MSE deterioration in case of low \(\sigma\) (not noisy) is expected.
Moving to the SSIM index, which is mostly related to details preservation, MAP and BM3D approaches are characterized by similar performances, with an advantage of the former in case of strong noise (\(\sigma \ge 5\), SNR \(\le\) 13 dB) and of the latter in case of weak noise (\(\sigma \le 3\), SNR \(\ge\) 17 dB), as expected again.
By looking at the residual error maps of Fig. 4, which refer to the \(\sigma =4\) case, different behaviors of the considered approaches can be found. In particular, the anisotropic diffusion and the LMMSE filters show a residual error mostly concentrated on edges, indicating poor detail preservation capability and confirming the low values of SIM reported in Fig. 5. The 3D bilateral and the MAP filters, show small errors across the whole image (smooth areas and details). Visually, the errors of the proposed MAP approach are lower than the 3D bilateral filter, as expected from the MSE values. Moreover, they are mainly concentrated on the outer edges: this result was expected, because of the lack of a useful neighborhood on the external edges. The BM3D residual error lays in the middle of the two categories, as the error intensity is not low but at the same time it is not concentrated close to the edges.
Moving to the 3 T real dataset, let us focus on Fig. 6 and on the enlargements over a region of interest (ROI) placed in the top right corner, reported in Fig. 7. As in the simulated case study, the LMMSE filter oversmooths the images, loosing many details, while the 3D bilateral approach is not very effective in reducing noise. Concerning the anisotropic diffusion filter, interesting results both in terms of noise reduction and edges preservation are achieved. This is in contrast with the poor results of the simulated dataset previously reported. Such behavior can be explained considering the high supervision required by the approach: the same configuration parameters (the default ones) have been adopted for both real and simulated test cases, providing results with very different accuracy in the two case studies. This aspect greatly influences the applicability of the Anisotropic Diffusion filter. The BM3D algorithm result is very clean, and noise reduction is considerable, although its behavior is opposite with respect to anisotropic diffusion: the good results of the simulated dataset are not confirmed. As in the previous case, the effect is probably due to the several configuration parameters required by the algorithm. We underline that the default parameters adopted by the authors in [10] were considered in this manuscript.
Moving to the proposed MAP approach, again the filtered image is a good compromise between noise reduction and details preservation, providing good reconstructions without requiring any parameter setting. In particular, looking at the enlargement reported in Fig. 7, it is clear the strong smoothing of LMMSE and BM3D approaches, and the good details preservation of anisotropic diffusion, 3D bilateral and MAP approaches, with the latter appearing a good tradeoff between smoothness and sharpness.
The real 1.5 T case axial study mainly confirms previous results. With respect to 3 T case, a higher noise is corrupting the images. The LMMSE in this case fails in detecting the smooth areas, producing a very poor regularization effectiveness. That said, again the performances of 3D bilateral and proposed MAP approaches are similar, with slight higher noise removal effectiveness for the latter. anisotropic diffusion solution is very smooth, while BM3D approach produces a result that is in between LMMSE and 3D bilateral in terms of details preservation and noise smoothing. In the enlargements reported in Figs. 9 and 10 it is evident that MAP is also effective in handling smooth transition, with the boundary between white matter and gray matter that is well recognized. Moreover, a lot of small details are correctly retrieved in the filtered image.
Moving to the second 1.5 T dataset, the sagittal acquisition of a portion of the vertebral column, results over two consecutive slices are reported. LMMSE and BM3D algorithms produce an oversmoothed solution, with some fine details lost across the image. On the other side, The anisotropic diffusion approach preserves a lot of details at the cost of not very effective noise removal. Among all, the 3D bilateral approach produces the most similar results with respect to the MAP algorithm. Moving to the ROI of Fig. 13, it allows to better appreciate the differences among algorithms. In particular, it is evident that only the proposed MAP filter is able to correctly retrieve the spinal cord texture.
Conclusions
Within this manuscript, a denoising methodology for MR images has been presented. The algorithm implements a maximum a posteriori (MAP) estimator and models the 3D acquired data via a local Gaussian Markov random field (LGMRF). The peculiarity of LGMRF consists in its ability to adapt itself to the local behavior of the imaged slices, allowing the preservation of small anatomical structures and edges filtering without any supervision. Simulated and real data results show that, compared to other widely adopted MRI denoising algorithm such as linear minimum mean square error (LMMSE), blockmatching and 3D (BM3D), 3D bilateral and nonlinear anisotropic diffusion filters, a good noise removal effectiveness is achieved by MAP together with interesting performances in terms of edges and details preservation. In particular, compared to 3D bilateral and BM3D filters it showed similar edge preservation effectiveness, but with higher regularization capability in smooth areas without manual set of any filter parameter.
The achieved results together with its unsupervised nature potentially make the proposed approach an interesting and promising instrument for denoising clinical images.
Issues regarding some implementation aspects of the proposed methodology are still open, and will be addressed in the future. In particular, further studies on the optimal size of the neighborhood system and on the improvement of the computational efficiency of the filter.
Abbreviations
 BM3D:

blockmatching and 3D filtering
 DT:

diffusion tensor
 LGMRF:

local Gaussian MRF
 LMMSE:

linear minimum mean squared error
 MAP:

maximum a posteriori
 MLE:

maximum likelihood estimator
 MR:

magnetic resonance
 MRF:

Markov random field
 MRI:

magnetic resonance imaging
 MSE:

minimum square error
 ROI:

region of interest
 SNR:

signal to noise ratio
 SSIM:

Structural Similarity Index
Declarations
Authors’ contributions
All authors participated in developing the method, testing on considered dataset and writing the manuscript. All authors read and approved the final manuscript.
Acknowledgements
Authors would like to thank Rocchina Caivano and Aldo Cammarota, Department of Radiology, IRCSS CROB, Rionero in Vulture, Italy, for helping in the clinical evaluations of the obtained results.
Competing interests
The authors declare that they have no competing interests.
Ethics approval and consent to participate
Written informed consent was obtained from the patients for publication of this manuscript and any accompanying images.
Data availability
The simulated dataset is available within Mathworks™Matlab^{®} software. The 1.5 T dataset is available on the dicom.nema.org website. Please contact author for requests about the 3 T dataset.
Funding
No source of funding supported this research.
Open AccessThis article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
Authors’ Affiliations
References
 Mohan J, Krishnaveni V, Guo Y. A survey on the magnetic resonance image denoising methods. Biomed Signal Process Control. 2014;9:56–69.View ArticleGoogle Scholar
 Buades A, Coll B, Morel JM. A review of image denoising algorithms, with a new one. IAM J Multiscale Model Simul SIAM Interdiscip J. 2005;4(2):490–530.MathSciNetView ArticleMATHGoogle Scholar
 Balafar MA. Review of noise reducing algorithms for brain MRI images. Int J Tech Phys Probl Eng. 2012;4(4):54–9.Google Scholar
 Vaishali S, Rao KK, Rao GVS. A review on noise reduction methods for brain MRI images. In: International conference on signal processing and communication engineering systems (SPACES), 2015, p. 363–5; 2015.
 Gallea R, Ardizzone E, Pirrone R, Gambino O. Noise filtering using edgedriven adaptive anisotropic diffusion. In: 21st IEEE international symposium on computerbased medical systems, 2008. CBMS ’08. p. 29–34; 2008.
 Anand CS, Sahambi JS. MRI denoising using bilateral filter in redundant wavelet domain. In: TENCON 2008–2008 IEEE region 10 conference, p. 1–6; 2008.
 MartinFernandez M, AlberolaLopez C, RuizAlzola J, Westin CF. Sequential anisotropic Wiener filtering applied to 3D MRI data. Magn Reson Imaging. 2007;25(2):278–92.View ArticleGoogle Scholar
 Baselice F, Ferraioli G, Johnsy AC, Pascazio V, Schirinzi G. Speckle reduction based on wiener filter in ultrasound images. In: 2015 Annual international conference of the IEEE engineering in medicine and biology society (EMBC) (2015).
 Sijbers J, den Dekker AJ, Scheunders P, Van Dyck D. Maximumlikelihood estimation of rician distribution parameters. IEEE Trans Med Imaging. 1998;17(3):357–61.View ArticleGoogle Scholar
 Dabov K, Foi A, Katkovnik V, Egiazarian K. Image denoising by sparse 3d transformdomain collaborative filtering. IEEE Trans Image Process. 2007;16(8):2080–95.MathSciNetView ArticleGoogle Scholar
 Iftikhar MA, Rathore S, Jalil A, Hussain M. A novel extension to nonlocal means algorithm: application to brain MRI denoising. In: 2013 16th international multi topic conference (INMIC), p. 195–200; 2013.
 Descombes X, Kruggel F, von Cramon DY. fMRI signal restoration using a spatiotemporal markov random field preserving transitions. NeuroImage. 1998;8(4):340–9.View ArticleGoogle Scholar
 Baselice F, Ferraioli G, Pascazio V. A Bayesian approach for relaxation times estimation in MRI. Magn Reson Imaging. 2016;34(3):312–25.View ArticleGoogle Scholar
 MartinFernandez M, Westin CF, AlberolaLopez C. 3D Bayesian regularization of diffusion tensor MRI using multivariate Gaussian Markov random fields. In: Seventh international conference on medical image computing and computerassisted intervention (MICCAI’04), vol 7, no Pt 1, p. 351–9; 2004.
 Baselice F, Ferraioli G, Grassia A, Pascazio V. Optimal configuration for relaxation times estimation in complex spin echo imaging. Sensors. 2014;14(2):2182.View ArticleGoogle Scholar
 Baselice F, Ferraioli G, Pascazio V. Relaxation time estimation from complex magnetic resonance images. Sensors. 2010;10(4):3611–25.View ArticleGoogle Scholar
 Baselice F, Caivano R, Cammarota A, Ferraioli G, Pascazio V. T1 and T2 estimation in complex domain: first results on clinical data. Concepts Magn Reson A. 2014;43(5):166–76.View ArticleGoogle Scholar
 Sijbers J, den Dekker AJ, Van Audekerke J, Verhoye M, Van Dyck D. Estimation of the noise in magnitude MR images. Magn Reson Imaging. 1998;16(1):87–90.View ArticleGoogle Scholar
 Ambrosanio M, Baselice F, Ferraioli G, Lenti F, Pascazio V. Intra voxel analysis in magnetic resonance imaging. Magn Reson Imaging. 2017;37:70–80. doi:10.1016/j.mri.2016.11.009.View ArticleGoogle Scholar
 Baselice F, Ferraioli G, Pascazio V, Schirinzi G. Contextual informationbased multichannel synthetic aperture radar interferometry: addressing DEM reconstruction using contextual information. IEEE Signal Process Mag. 2014;31(4):59–68.View ArticleGoogle Scholar
 Saquib S, Bouman C, Sauer K. ML parameter estimation for Markov random fields, with applications to Bayesian tomography. IEEE Trans Image Process. 1998;7:1029–44.View ArticleGoogle Scholar
 Baselice F, Ferraioli G, Shabou A. Field map reconstruction in magnetic resonance imaging using Bayesian estimation. Sensors. 2010;10(1):266–79.View ArticleGoogle Scholar
 Li SZ. Markov random field modeling in image analysis. Secaucus: Springer; 2001.View ArticleMATHGoogle Scholar
 Geman S, Geman D. Stochastic relaxation, Gibbs distributions, and the Bayesian restoration of images. IEEE Trans Pattern Anal Mach Intell. 1984;6:721–41.View ArticleMATHGoogle Scholar
 Gerig G, Kubler O, Kikinis R, Jolesz FA. Nonlinear anisotropic filtering of MRI data. IEEE Trans Med Imaging. 1992;11(2):221–32.View ArticleGoogle Scholar
 Paris S, Durand F. A fast approximation of the bilateral filter using a signal processing approach. In: 9th European conference on computer vision. 2006.
 AjaFernandez S, Niethammer M, Kubicki M, Shenton ME, Westin CF. Restoration of DWI data using a rician LMMSE estimator. IEEE Trans Med Imaging. 2008;27(10):1389–403.View ArticleGoogle Scholar
 Wang Z, Bovik AC, Sheikh HR, Simoncelli EP. Image quality assessment: from error visibility to structural similarity. IEEE Trans Image Process. 2004;13(4):600–12.View ArticleGoogle Scholar