Computational and Mathematical Methods in Medicine

Volume 2015, Article ID 182659, 21 pages

http://dx.doi.org/10.1155/2015/182659

## The EM Method in a Probabilistic Wavelet-Based MRI Denoising

^{1}Laboratorio de Procesado de Imagen, Escuela Técnica Superior de Ingenieros de Telecomunicación, Campus Miguel Delibes s.n., 47011 Valladolid, Spain^{2}Departamento de Álgebra, Análisis Matemático, Geometría y Topología, Facultad de Ciencias, Campus Miguel Delibes s.n., 47011 Valladolid, Spain

Received 1 July 2014; Revised 23 September 2014; Accepted 30 September 2014

Academic Editor: Yi Gao

Copyright © 2015 Marcos Martin-Fernandez and Sergio Villullas. 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

Human body heat emission and others external causes can interfere in magnetic resonance image acquisition and produce noise. In this kind of images, the noise, when no signal is present, is Rayleigh distributed and its wavelet coefficients can be approximately modeled by a Gaussian distribution. Noiseless magnetic resonance images can be modeled by a Laplacian distribution in the wavelet domain. This paper proposes a new magnetic resonance image denoising method to solve this fact. This method performs shrinkage of wavelet coefficients based on the conditioned probability of being noise or detail. The parameters involved in this filtering approach are calculated by means of the expectation maximization (EM) method, which avoids the need to use an estimator of noise variance. The efficiency of the proposed filter is studied and compared with other important filtering techniques, such as Nowak’s, Donoho-Johnstone’s, Awate-Whitaker’s, and nonlocal means filters, in different 2D and 3D images.

#### 1. Introduction

Magnetic resonance imaging (MRI) is one of the most important imaging acquisition techniques [1], which allows studying the structural features of the internal body parts noninvasively. This procedure is based on the principle of nuclear magnetic resonance (NMR) [2], and the power of this technique over other noninvasive techniques, such as ultrasound, is the high quality of its images, despite the inconvenience of being a larger and expensive equipment. However, because the ultrasonic signal is not transmitted through bones, in cases such as imaging the brain which is surrounded by the skull, ultrasound modalities are not viable and resorting to magnetic resonance (MR) imaging is needed.

For a given acquisition time, in MRI there is a fundamental agreement between resolution and signal to noise ratio (SNR) [3]. MRI is affected by noise mainly produced by interference due to human body heat emission (Gaussian at frequency space and Rayleigh at envelope of its inverse Fourier transform [4]), which prevents correct identification of shapes and details. Moreover, there is a relationship between noise level and image resolution in MRI acquisition, that is, the larger the resolution of the acquired image, the lower the SNR [5]. The simplest method to reduce the noise level is to increase the acquisition time of the machine (thus increasing the number of images averaged, that is, increasing the machine number of experiments (NEX)), which would cause a large increase in spending and long waiting lists, but long acquisition times could be problematic for patients who are not able to remain in a resting state (due to stress, pain, and so on). To avoid these problems, a filter, which acts on the acquired image, can be applied. This filter must eliminate noise trying to preserve details. The main problem of removing noise by means of softening an image is the resulting loss of information on the edges and contours (image blur), which is typical in Gaussian filter convolution. Perona and Malik [6] proposed a new type of filtering, based on partial differential equations and the heat diffusion equation, which is the origin of a family of filters that allow homogenizing regions while maintaining or enhancing borders between them. Gerig et al. [7] propose the use of the so-called nonlinear anisotropic filter, which gives very good results in the context of MRI and is an important practical application of the ideas proposed by Perona and Malik. From the same diffusion equation, [7] proposes an alternative discretization with simpler formulation, whose stability is subject to certain restrictions of the parameters. Linearly optimal methods in the sense of minimum mean square error (Wiener filter) have also been adapted to the case of MRI [8, 9]. Another filter with good MR image denoising results is the so-called nonlocal means (NLM) filter [10, 11], which averages similar image pixels as a function of their intensity distance (some filters, like the bilateral filter [4], are based on the same proposition, but the advantage of the NLM over other methods is that the similarity measure used is more robust to noise due to region comparisons rather than pixel comparisons). Principles of nonparametric statistical methods are also the base of the iterated conditional entropy reduction (ICER) proposed by Awate and Whitaker [12], a Bayesian-inference algorithm based on Markov random field which estimates the uncorrupted-image statistics by optimizing an information-theoretic metric using the expectation-maximization algorithm. This filter method incorporates a Rician noise model, unlike NLM method, which is more general. He and Greenshields [13] designed another filter method that improves NLM filter by adding Rician noise information. Dabov et al. [14] also proposed a filter similar to NLM. This method creates 3D arrays formed by stacking together similar image 2D neighborhoods. The importance of grouping is to enable the use of a higher-dimensional filtering of each group, which exploits the potential similarity (correlation, affinity, etc.) between grouped fragments. More generally, Sivaramakrishnan and Weissman [15] designed a universal filter that does not need* a priori* noise information which is asymptotically optimal. Furthermore, Awate and Whitaker also proposed a patch-based method [16, 17] that tries to optimize the entropy of the noisy image to reduce noise.

A large interesting characteristic of the wavelet transform is its capacity to preserve detail at different scales, because of its ability to model the information locally present in the image due to the multiresolution decomposition [18]. In the literature, other authors have performed MR image filtering using wavelet techniques. Donoho and Johnstone [19] proved that a simple thresholding algorithm with an appropriated base may be a (nonlinear) filter which is almost optimal. Nowak [20] and Pizurica et al. [21] propose to perform the filtering using the discrete wavelet transform. In particular, the significance of Nowak’s work is that it uses the fact that MR magnitude image obeys a Rician distribution and its square image noise obeys a noncentral Chi-square distribution. In addition, other researches, like Sijbers et al. [22], use this fact. A diffusion method such as Perona-Malik’s, but adapted to the Rician distribution case, has been proposed in [23]. In [20], the wavelet transform is performed on the square of the amplitude image where the noise and the bias of approximation coefficients are reduced. Another interesting work is that of Anand and Sahambi [4]. In this case, the square of the amplitude image is also used, correcting the bias and applying a bilateral filter (Gaussian filtering in the spatial and amplitude domains) over approximation coefficients (it is also based on the Rician image distribution). On the other hand, Yang and Fei [24] combined the 1D wavelet transform with the Radon transform to denoise Rician noise in MR images. Finally, we can also mention the method proposed by Wirestam et al. [25], where a technique of shrinkage coefficients of the filter is based on a Wiener filter. The novelty of this method is that the filtering is performed in the complex Fourier image, where noise data are complex Gaussian. This method has the problem that complex data are not always available in MRI acquisition (usually, MR machines provide only data in the image plane after envelope calculation in DICOM format. Complex raw data in the -space are normally stored in a proprietary format which is not open and is brand-dependent. The complex data in the image plane are not even stored in the machine).

Based on wavelets state of art and given locality property of wavelet transform, the alternative that we propose in this paper is to perform filtering in the domain of the transform coefficients, adapting to MRI the method proposed in [26], originally designed for mammographic images. We also propose to estimate the parameters of the model by means of the expectation maximization (EM) method, proposed by Dempster et al. [27] and Moon [28]. Moreover, this fact makes the filter independent of noise variance estimators, in contrast to other filter methods.

This paper is structured as follows. In Section 2, some different MR image denoising methods are detailed. First, in Section 2.1, the methods proposed by Donoho and Johnstone [19] (Section 2.1.1) and Nowak [20] (Section 2.1.2) are described. Section 2.1.3 shows shift variance of wavelets coefficients and how to correct it. This section also contains Section 2.2 (with Sections 2.2.1 and 2.2.2), where Awate-Whitaker’s algorithm [16, 17] and nonlocal means filter [10, 11] are presented. Finally, Section 2.3 describes the used technique to estimate noise variance. In Section 3, the new wavelet denoising method, based on [26], is presented. Section 4 presents the algorithms corresponding to the methods showed in Sections 2 and 3 and some practical experiments. First, Section 4.1 shows a step-by-step explanation of the different denoising algorithms presented before. Second, in Section 4.2, the images used to test the efficiency of the proposed method are detailed. The measurements used to compare the different filter methods are described in Section 4.3. Fourth, in Section 4.4, the numerical experiments and some remarks, obtained after performing these experiments, are presented. Finally, Section 5 contains a brief summary of the obtained results and some possible future research lines. Appendix A contains the details about the optimization method proposed to estimate the parameters of the filter presented in Section 3 and Appendix B contains some useful auxiliary functions.

#### 2. State of Art

This section reviews some filters used in the practical experiments of Section 4. The new filter proposed in this paper is compared with other two wavelet filters, Donoho-Johnstone’s hard thresholding filter [19] and Nowak’s [20] filter, the Awate and Whitaker’s [16, 17] patch-based (Gaussian and Rician) filters and the nonlocal means [10, 11] filter. Moreover, wavelet-based filtering is shift-variant, as Section 2.1.3 shows. This subsection also presents how to avoid this problem. Finally, an estimator of noise variance, required in both wavelet-based filters, is given in Section 2.3.

##### 2.1. Wavelet-Domain Filters

An image/volume can be interpreted as a 2-dimensional/3-dimensional function with compact support. The values of this function, represented in a matrix/3D array , are a good approximation to scale coefficients, , in discrete wavelet transform. The fast wavelet transform algorithm let us calculate the scale, , and detail, , coefficients in the following levels , that is, if , and can be calculated in function of . Noise interferences modify the details of the MR image/volume; as noise grows more levels are affected. The wavelet coefficients are filtered to denoise the image/volume. The wavelet coefficients (the index corresponds to level and orientation: horizontal, vertical, diagonal, and so on, and index corresponds to scale and position) can be determined by means of where is the discrete wavelet function at level and orientation and scale and position, and index represents pixel/voxel position. Similarly the scaling coefficient (the index corresponds to level and index corresponds to scale and position) can be determined by means of where is the scaling function for level and scale and position. Given a sequence of wavelet coefficients and scaling coefficients for the image/volume (where represent the number of scales and positions for each and is the number scales and positions for each ) the filters can be defined in the wavelet domain.

The following two filters, proposed by Donoho and Johnstone and Nowak, are used to analyze the efficiency of the new filter presented in next section.

###### 2.1.1. Donoho-Johnstone’s Filter

The classic hard thresholding filter described by Donoho and Johnstone [19] is given bywhere is the module operator and with standard deviation of the noise in the image/volume .

###### 2.1.2. Nowak’s Filter

Another filter is proposed by Nowak [20] bywhere and is the variance of the wavelet coefficient . An estimation for that, , is proposed in [20] where is the standard deviation of the noise in the image/volume .

###### 2.1.3. Shift-Invariant Filtering

As [20] shows, wavelet coefficients filtering (based on discrete wavelet transform) is shift-variant, which can produce the appearance of artifacts, because the wavelet coefficient values depend on the alignment between the data and the wavelet basis functions. Shift-invariant (translation-invariant, undecimated) methods [29–33] can provide better performance, avoiding the appearance of artifacts (see an example in Figure 1). Shift-invariant wavelet transforms provide a higher degree of regularity [29, 33, 34] than standard wavelet analysis approaches, so shift-invariant estimation algorithms usually outperform standard methods. A shift-invariant filtering can be obtained by applying the filter for every possible shift of the image, unshifting each filtered result, and averaging all the results obtained. As performing all shifts of the image would be computationally too expensive, to reduce computational burden of filtering, an approximately shift-invariant scheme is proposed; that is, all shifts are replaced by a small range of shifts. More specifically, the image is shifted in both horizontal and vertical directions (left, right, up, and down), at most pixels/voxels in each direction in steps of one pixel/voxel. While this does not guarantee shift-invariance, it does reduce the dependence of the filter output on the alignment between the data and the wavelet basis functions. In our experiments, we used , that is, movements of in each direction giving rise to a total of shiftings for the 2D images (the values enlarge computational time excessively with minor improvements), and , that is, movements of in each direction giving rise to a total of shiftings for the 3D volumes (the same remark as in 2D also applies in 3D for ). The final image/volume is constructed after unshifting each image/volume and averaging the resulting images/volumes.