The EM Method in a Probabilistic Wavelet-Based MRI Denoising

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.


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 2 Computational and Mathematical Methods in Medicine 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 expectationmaximization 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 branddependent. 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.

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.

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, 0 , in discrete wavelet transform. The fast wavelet transform algorithm let us calculate the scale, , and detail, , coefficients in the following levels 0 ≤ ≤ , that is, if 0 = , +1 and +1 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 = ( [ ]) =1 and scaling coefficients = ( [ ]) =1 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.

Donoho-Johnstone's Filter.
The classic hard thresholding filter described by Donoho and Johnstone [19] is given by where | ⋅ | is the module operator and := noise √2log ( ) with noise standard deviation of the noise in the image/ volume . [20] by

Nowak's Filter. Another filter is proposed by Nowak
where and An estimation for that,( ) 2 [ ], is proposed in [20] where noise is the standard deviation of the noise in the image/volume .

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. Shiftinvariant (translation-invariant, undecimated) methods [29][30][31][32][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 shiftinvariant 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 = 2, that is, movements of {−2, −1, 0, 1, 2} in each direction giving rise to a total of 25 shiftings for the 2D images (the values > 2 enlarge computational time excessively with minor improvements), and = 1, that is, movements of {−1, 0, 1} in each direction giving rise to a total of 27 shiftings for the 3D volumes (the same remark as in 2D also applies in 3D for > 1). The final image/volume is constructed after unshifting each image/volume and averaging the 25/27 resulting images/volumes.

Awate-Whitaker's Patch-Based Filter and Nonlocal Means
Filter. For quantitative evaluation, these two wavelet methods (and the new proposed filter) are compared (in 2D case) with a patch-based method proposed by Awate and Whitaker [16,17] and the nonlocal means (NLM) filter [10,11]. This comparison will be very significative as this method was recently shown to be competitive with respect to wavelet methods.

Awate-Whitaker's Filter.
In this approach, given an image , a random vector ( ) = [ ( ), ( )], where represents the pixel position on the image , is generated, where ( ) is the intensity of at the pixel and ( ) is the intensity of at the pixels in a neighborhood ( ( ) is a vector) of pixel (we usẽ( ),̃( ), and̃( ) for degraded image random variables). The target of the method is to minimize the entropy of the conditional PDF, ℎ(̃|̃), for which a descending gradient method, given bŷ where is a parameter, is used. In this paper, we used two versions of this filter method using Gaussian and Rician models for the PDF.

Nonlocal Means Filter.
For a given image , the NLM filtered image at pixel position ( NLM is the NLM filter operator) is given by the weighted average of all the pixels in a searched area Ω of pixel position in the image , where 0 ≤ ( , ) ≤ 1, ∑ ∈Ω ( , ) = 1. The weights ( , ) are based on the similarity between the neighbourhoods of pixels ( ) and ( ) and are defined as where and are the neighborhoods of the pixel positions and , respectively, is a Gaussian weighted squared Euclidean distance, and ℎ is the exponential decay control smoothing parameter. Region Ω can be the whole image, but, because of computational reasons, Ω uses to be a smaller region in the local neighborhood.

Estimation of 2
noise . Donoho and Johnstone's [19] and Nowak's [20] methods are highly dependent on the noise estimate. We assume that noise distribution, in complex domain, is zero mean Gaussian. Then, noise estimate means noise standard deviation/variance ( noise / 2 noise ) estimate. In MR images, noise is usually unknown a priori and it must be estimated from the data. A good estimator of 2 noise of is given in [35,36] where 2 local is the local variance of , defined as where the operator LV is defined in Appendix B.

A New Probabilistic Wavelet Filter
The wavelet filter proposed in this paper, which will be named hereafter as Villullas-Martin's filter, is defined by the shrinkage where are the wavelet coefficients at level/orientation of an image/volume and with parameter ( [ ] represents the posterior probability of being detail given [ ] with the prior probability of a coefficient being noise). This filter is based on the filter proposed by Gorgel et al. [26]. The novelty of our filter is the proposed distribution models for details and noise coefficients in the wavelet domain to MRI and the use of the EM method for estimating the parameters. This avoids the problem setting any free parameter such as the noise variance which is usually problematic. Noiseless MR images/volumes have a distribution in the wavelet domain with a pronounced maximum at the origin (due to smooth regions) and long tails produced by edges and different structures contained in the image/volume (as shown in Figure 2), so this distribution can be approximately modeled by a Laplacian function, given by with parameters and . This distribution was first proposed to model wavelet coefficients distribution of mammographic images by Gorgel et al. [26]. MR magnitude image/volume distribution is Rician [5,37]. In high signal to noise ratio (high intensity, bright) regions, Rician distribution tends to a Gaussian distribution and in low signal to noise ratio (low intensity, dark) regions, Rician distribution tends to a Rayleigh distribution [4]; that is, MR noise (signal free) can be modeled by a Rayleigh distribution. In the wavelet domain, the distribution of this noise has a maximum at the origin with short tails (see Figure 3). In this case, we approximate this distribution by a Gaussian distribution, given by The definition of as a conditioned probability lets us modify the wavelet coefficients taking into account noise intensity. Besides, the larger image/volume modification by wavelet coefficients shrinkage, the bigger module of the shrunk coefficients; that is, a change in low module wavelet coefficient does not change significantly the image/volume. So, the importance of is focused on the tails of the distributions. Figure 4 shows the mixture model of detail and noise, with the two approximated distribution models presented before, superimposed on real image/volume wavelet histograms at different levels and orientations.
Parameters estimation in the original paper [26] consider both distributions independently choosing the parameter depending on the ratio noise/detail. In this paper we consider the joint distribution resulting from considering both distributions allowing us to find the parameters that maximize the similarity between the real and the theoretical distributions. For this purpose, parameters , , noise , and are to be calculated by the EM method as Appendix A shows in detail. Here a brief summary is given as follows; given the (independent) known data X = { } =1 (vector X represents each set of wavelet coefficients for each scale and orientation of an image/volume. We also drop the upper index for the sake of simplicity along this the likelihood function to be maximized in this method is where Θ = [ , , noise , ] and ( | Θ) = Gauss ( | noise ) + (1 − ) Laplace ( | , ), and its expected value is Maximizing this expression, we obtain It is an implicit system, so estimators ,̂,̂2 noise ,̂, andâ re calculated by fixed-point iteration with initial conditionŝ Ini := 1 2 , Ini := median (X) , where 2 local := LV(X) is defined in (B.2) (whereX is the vector X as a matrix [of wavelet coefficients at scale and orientation ]), and

Experiments
In this section, some noisy MR 2D and 3D volumes are filtered to compare the algorithm proposed in this paper with other wavelet filters as well as patch-based and NLM filters found in the reviewed literature and described in Section 2. (ii) Compute the ( = 2)-scale DWT of the squared magnitude image/volume 2 .
(iii) Remove the bias from the scaling coefficients by subtracting = 2 +1 noise from each (see [20]). (iv) Filter the wavelet coefficients through the Nowak's filter, N , defined in Section 2.1.
(v) Compute the inverse DWT of the filtered wavelet and unbiased scaling coefficients to obtain an estimate of squared denoised image/volume.
(iv) The image/volume +1 is given by the intensitŷ+ 1 at pixel . If + 1 is less than the maximum number of iterations, go to Step (ii) with = + 1. Otherwise, the filtered image/volume is +1 .

Nonlocal Means Algorithm
(i) Given the parameters values, calculate the weights function for each pixel of the image/volume .
(ii) For each pixel, compute the NLM filtered image/ volume with the weight function.

MR Data Sets.
The experiments were conducted on fifth MRI data sets. The first and second data sets consist of simulated MR volumes and images obtained from the Brainweb database [38]. The third and fourth data sets were collected from Centro de Diagnóstico Valladolid (CDV) QDIAGNOSTICA in Valladolid (Spain). The last data set was obtained from a database of the Laboratory of Image Processing (LPI) of the University of Valladolid. Details about these data sets are described in following subsections.

Simulated MR Images.
Simulated MR images/volumes are a useful data set which allows a first evaluation of different analysis methods. In the 3D case, the noiseless image volume consists of a 3-dimensional volume of resolution 180×216×4 extracted from a volume generated in the Brainweb database of resolution 181 × 217 × 181, T1-weighted, 1 mm slice thickness, 0% of noise, and RF = 0%. Image intensity values vary in {0, 1, . . . , 255}. In the 2D case, the noiseless image consists of a 2-dimensional axial section of the 3-dimensional volume of resolution 180 × 216. Rician noise was generated by the equation̂: where is the noiseless MR image/volume and / ∈ {1, 2} are (0, 2 ) independent identically distributed Gaussian random variables, with ∈ {5, 6, . . . , 20}.

Real MR Images.
Real data set consists of two data sets. Real data images are 2D axial section of a brain, acquired using a General Electric Signa 1. . In this case, as no ground truth is available, the evaluation will be performed based on experts' rankings.

Quality Measures.
After filtering these images/volumes ( is the noiseless image/volume, is the noisy image/volume, and is the filtered noisy image/volume), the results obtained are compared using the following efficiency measurements (see Appendix B for details).

Averaged Error Local Variance (AELV).
AELV is an objective quality measure [39] that quantifies the deviation of estimated values from the true value. Specifically, the AELV of with respect to is measured as where is the total number of pixels/voxels of and the operator LV is defined in Appendix B.

Normalized Averaged Error Local Variance (NAELV).
Variations in AELV can be difficult to understand. To test at what rate reduces its value after filtering, AELV normalization is used as

Structural Similarity (SSIM).
Although AELV is a useful measure of similarity, it is not suitable to obtain a comparison similar to that performed by the human eye [40,41]. Most common alternative is SSIM which is consistent with the visual perception. The SSIM index is estimated as where the operators GM( ), GV( ), and GCV( , ) are defined in Appendix B, 1 := 6.5025 and 2 := 58.5225.

Averaged Local Signal to Noise Ratio (ALSNR).
Another simple method of checking the noise level in an image is the averaged local SNR, measured as where is the total number of pixels/voxels of and the operator LV is defined in Appendix B.

Contrast. The measure normally used to calculate images contrast is given by
where max and min are the maximum and minimum value in a specific ROI of the image/volume . This quality measurement is not effective enough because some filters (e.g., some wavelet filters) generate a bias that modifies the values of this quality measurement but does not affect the contrast of the image/volume. Besides, there are many tools, used by radiologists, which let them modify the window/level of the image/volume. We can observe the variation of contrast taking a section of the corresponding images/volumes to compare the different changes of intensity. Figure 5 shows that wavelet filters preserve edge changes of intensity very well due to its property of locality.  Figure 6 shows the comparison of the different averaged quality measurements as function of the parameter (each parameter has associated 20 simulated images experiments and the corresponding quality measurement value is the average of the 20 quality measurement values of the corresponding images. We use this 20 data sets to reduce the variability of the quality measurements). For low values of the parameter (almost noiseless images), Awate-Whitaker's filter (Rician over Gaussian in SSIM measurement and Gaussian over Rician in AELV measurement) have a good noise removal and it beats the other methods in this case. Nowak's filter gets results close to Villullas-Martin's filter (in fact, Nowak's filter improves Villullas-Martin filter in SSIM for tiny values of ) but as parameter value increases, difference between Villullas-Martin's filter and Nowak's filter grows. Villullas-Martin's filter is clearly the best method for medium and high values of (hard noisy images. This value corresponds to low NEX in image acquisition). NAELV graph shows us that Awate-Whitaker's filter has constant improvement proportion in contrast to wavelet filters, which amplify its improvement proportion with . Then, Awate-Whitaker's filter is a good denoising method for very low values of but as increases, the filter's strength decreases. The NLM filter has a behavior near to Nowak's filter, better than Nowak's filter in high values of and worst in low and medium values of . Figure 7 shows an example ( = 15) where visual differences between the filtering methods can be seen. In this case, Villullas-Martin's filtering and Nowak's filtering are similar to each other and better than Donoho-Johnstone's filtering which has lower noise removal. NLM filter is the worst choice because the noise has not been properly eliminated inside brain structure.

Experiment 2: Filtering Real Data.
In the second experiment, we have filtered the first real data set described in Section 4.2.2. Noisy images set is the images belonging to NEX ∈ {4, 8, 16} and noiseless image is generated with NEX = 64. Figure 8 shows the comparisons of the different quality measurements as function of NEX. It can be seen that Villullas-Martin's filter improves Nowak's and Donoho-Johnstone's filters, with a slight increase as NEX increase. Awate-Whitaker's filters cannot improve wavelet filters results (moreover, as AELV graph shows, Awate-Whitaker's filter ruins low NEX images). As Experiment 1, NLM filter works as good as Nowak's filter. Figure 9 shows the example of NEX = 8. In this example we can see that Villullas-Martin's filter provide a better filtering, where more noise is removed than Nowak's, Donoho-Johnstone's, and NLM filters. Awate-Whitaker's filter with Gaussian model obtains a bit worse denoising than Nowak's, Donoho-Johnstone's, and NLM filters and Awate-Whitaker's filter with Rice model does not remove enough noise inside brain structures and generates artifacts.

Experiment 3: 3D Wavelet
Filtering. This experiment shows the strength of the wavelet filters in the 3-dimensional case. In this case, the 16 3D simulated data volumes described in Section 4.2.1 for ∈ {5, 6, . . . , 20} are filtered. Figure 10 shows the comparison of the different averaged (20 volumes for each level whose quality measurements are averaged for each ) quality measurements as function of . This experiment shows the same behavior as Experiment 1 in wavelet filters. All of these methods improve the filtering in 3D case preserving the relationship between them.

Experiment 4: Filtered Real Image against Higher
NEX. The forth experiment evaluates the filtering power of Villullas-Martin's filter. This filter is compared with NEX rising. For this comparison, the second real data set described in Section 4.2.2 is used. The "same" 32 images with NEX = 1 let us control noise level, which can be reduced by averaging images in the complex domain, where higher values of imply lower noise level. The selected image to be filtered belongs to = 2. Table 1 shows that the first value of averaged images ( ) which improves the quality measurements of filtered image with = 2 is = 4; that is, Villullas-Martin's filter lets us obtain images with quality as good as with double acquisition time. Figure 11 shows the images involved in this experiment. It can be seen that filtered image with = 2 has less noise at smooth regions inside the structure than the image with no filtering with = 3.  much noise as other filters (apparently), which obtains similar filtering results in these cases, the result image given by our method was considered by the experts visually nicer and with higher contrast in the boundaries of the image, which eases radiologist work in identifying the different brain structures; that is, Villullas-Martin's filter preserves structures better than the other wavelet filters after denoising. The experiment results can be summarized as follows. In 84.2% of cases, Villullas-Martin's filter was chosen as the best filter method (rank = 1), whereas Nowak's filter was chosen as the worst filter in 88.3% of cases (rank = 3). Donoho-Johnstone's filter was chosen as the medium filter (rank = 2) in 75% of cases. Table 4 shows all percentages. Figure 12 shows some example images in this experiment.

Conclusion
A new wavelet-domain filtering has been proposed in this paper. Assuming that wavelet coefficients of a noiseless MR image/volume can be modeled by a Laplacian distribution (as we know in brain, the filter can be adapted to other human parts easily with minor or no changes), and the wavelet   coefficients distribution of Rayleigh noise can be approximated by a Gaussian distribution, a probabilistic method has been proposed; that is, wavelet coefficients are shrunk depending on its conditioned probability of being noise or detail (posterior probability). To calculate the parameters involved in the expression of VM , EM method has been used. This fact makes a filter independent of noise estimators,   As future research lines, we plan to generalize the distribution models of the wavelet coefficients both for the details and for the noise. For the details a good candidate will be the generalized Gaussian distribution which seems very promising. In addition to that, we plan to check whether the proposed method can be also applied to other body parts as well as to other MR modern modalities such as fast, ultrafast, and low-field MRI for which the noise level is known to be very high.

A. E.M. Parameter Estimation
We estimate the filter parameters by the EM method. This appendix shows the detailed derivation of the estimation equations. . (A.14) The estimator of ,̂cannot be determined using this technique as the expression ∑ =1 | ⋅ − | is piecewise-linear and its minimum is always in a vertex. So in this case a direct method is proposed: Finally the estimation of the shape parameter is given by where is the total number of pixels/voxels and represents the different pixel/voxel positions on the image/volume .