See discussions, stats, and author profiles for this publication at: [Link]

net/publication/236738925

An MRI denoising method using image data

redundancy and local SNR estimation

Article in Magnetic Resonance Imaging · May 2013

DOI: 10.1016/[Link].2013.04.004 · Source: PubMed

CITATIONS READS

11 34

3 authors, including:

Hosein M. Golshan
University of Denver
7 PUBLICATIONS 20 CITATIONS

SEE PROFILE

All in-text references underlined in blue are linked to publications on ResearchGate, Available from: Hosein M. Golshan
letting you access and read them immediately. Retrieved on: 19 November 2016
 Magnetic Resonance Imaging 31 (2013) 1206–1217

Contents lists available at SciVerse ScienceDirect

Magnetic Resonance Imaging

journal homepage: [Link]

An MRI denoising method using image data redundancy and local SNR estimation
Hosein M. Golshan a, Reza P.R. Hasanzadeh a,⁎, Shahrokh C. Yousefzadeh b
a
DSP Research Lab, Department of Electrical Engineering, University of Guilan, Rasht, Iran
b
Road Trauma Research Center, Guilan University of Medical Science, Rasht, Iran

a r t i c l e i n f o a b s t r a c t

Article history: This paper presents an LMMSE-based method for the three-dimensional (3D) denoising of MR images
Received 15 August 2012 assuming a Rician noise model. Conventionally, the LMMSE method estimates the noise-less signal values
Revised 6 April 2013 using the observed MR data samples within local neighborhoods. This is not an efficient procedure to deal
Accepted 6 April 2013
with this issue while the 3D MR data intrinsically includes many similar samples that can be used to
improve the estimation results. To overcome this problem, we model MR data as random fields and
Keywords:
Denoising
establish a principled way which is capable of choosing the samples not only from a local neighborhood but
Image data redundancy also from a large portion of the given data. To follow the similar samples within the MR data, an effective
Linear minimum mean square error similarity measure based on the local statistical moments of images is presented. The parameters of the
Magnetic resonance imaging proposed filter are automatically chosen from the estimated local signal-to-noise ratio. To further enhance
Rician distribution the denoising performance, a recursive version of the introduced approach is also addressed. The proposed
filter is compared with related state-of-the-art filters using both synthetic and real MR datasets. The
experimental results demonstrate the superior performance of our proposal in removing the noise and
preserving the anatomical structures of MR images.
© 2013 Elsevier Inc. All rights reserved.

1. Introduction Preservation of the meaningful image structures is a necessary

condition for developing MRI denoising filters. In other words, filtering
Magnetic resonance imaging (MRI) as a powerful medical methods should not remove useful anatomical information in their
imaging technique has an ever increasing importance in current intent to clean the noise. Moreover, the signal-dependent nature of
diagnosis and clinical applications. The ability of MRI in providing noise in MRI which is usually modeled following a Rician distribution is
accurate 3D representations of the internal body tissues and its non- another problem to deal with. Therefore, specific methods should be
invasive technology makes it more preferable for both patients and developed in order to manage MR images properly.
physicians [1]. Nevertheless, there is a trade-off between the There is a wide range of methods for the denoising of MR images.
imaging time and quality. Obtaining higher resolution images Several techniques have been proposed based on the partial
requires longer scanning time and averaging over multiple acquisi- differential equations (PDE). Considering the original anisotropic
tions, which is usually limited by parameters such as the patient diffusion (AD) filter introduced by Perona and Malik [3], Gerig et al.
comfort or by physiological constraints [2]. Hence, MR images often [2] presented a PDE-based filtering method for MR images. However,
suffer from low signal-to-noise ratio (SNR), which disturbs the visual the incorrect assumption of the noise model was a major drawback
inspection and also the diagnosis. of this approach. The noise model was taken to be Gaussian instead
Post-acquisition denoising of MR images using digital image of Rician. Consequently, a bias is introduced in the filtered image
processing techniques has always been one of the standard problems which is especially problematic in low SNR images. In order to
discussed in the literature. Denoising has been extensively used in remedy such a drawback, Sijbers et al. [4] proposed an adaptive AD
MRI to improve the image quality for more accurate diagnosis. Also, filter using a Rician noise model. Finally, Samsonov and Johnson [5]
it is used as a preprocessing step in many image processing addressed a noise adaptive nonlinear diffusion method for the
procedures such as segmentation and registration; as a consequence, denoising of MR images. All the aforementioned approaches usually
the quantitative measurements and computer aided analysis of the tend to erase fine details of the image and generate unnatural
data will be more robust and efficient. structures due to their edge enhancement effect. Several methods
have been presented to overcome these problems that are referred
to [6,7]. Most recently, Krissian and Aja-Fernandez [8] proposed a
⁎ Corresponding author.
E-mail addresses: [Link]@[Link] (H.M. Golshan),
noise-driven anisotropic diffusion filtering of MR images based on
hasanzadehpak@[Link] ([Link]. Hasanzadeh), yousefzadeh@[Link] the linear minimum mean square error (LMMSE) estimator that has
(S.C. Yousefzadeh). achieved state-of-the-art results.

0730-725X/$ – see front matter © 2013 Elsevier Inc. All rights reserved.
[Link]
 H.M. Golshan et al. / Magnetic Resonance Imaging 31 (2013) 1206–1217 1207

Many transform domain methods have also been presented for approach. Other modified versions of the ML method have been
noise suppression in MRI. Wavelet-based filters have been applied addressed in Refs. [25,26]. Wong et al. [27] presented Quasi-Monte
frequently to the denoising of MR images [9–11]. Wood et al. [12] Carlo estimation based on regional statistics. Aja-Fernandez et al.
have used wavelet packets for noise suppression in low SNR MR [28,29] proposed the LMMSE estimator for the Rician noise model.
images. In comparison with the single-basis wavelet decomposition, Considering the noise-less signal as a realization of a random
wavelet packet provides a more compact signal representation. variable, it leads to a closed-form analytical solution for the results.
Anand and Sahambi [13] addressed a wavelet domain nonlinear This makes the LMMSE estimator computationally efficient for
filtering that has obtained appropriate results. Other transforms processing the large 3D MR images.
[14,15] have also been applied to denoise images. Recently, Manjon In this paper, we present an LMMSE-based approach to the 3D
et al. [16] successfully developed an oracle-based discrete cosine denoising of magnitude MR images. It should be noted that the
transform (ODCT) for the denoising of 3D MRI. LMMSE method originally estimates the noise-less signal value using
During the recent years, the nonlocal means (NLM) filter local statistics of the image contents. In other words, the true value
proposed by Baudes et al. [17] has received considerable attention. for each noisy pixel is estimated by a set of pixels selected from a
This method has been developed based on the redundancy of local neighborhood. This scheme leads to a sub-optimal filtering
patterns within the natural images. The NLM filter tries to estimate performance, since the 3D MR data usually includes many similar
the noise-less signal value by weighted averaging over the input patterns that can be used to improve the estimation results. In the
image. These weights are computed using patch-based similarities current work, we address this issue by developing a nonlocal
among the pixels and that one to be estimated. Many MRI denoising processing of the LMMSE method. The presented approach takes
methods have been proposed based on the NLM filter such as the advantage of the high degree of redundancy in the contents of MR
works given by Wiest-Daessele et al. [18], Coupe et al. [19,20] and images and provides a suitable similarity measure to find the similar
Manjon et al. [16,21,22]. patterns within the given MR data. Furthermore, all the parameters
Apart from the aforementioned classes of methods, several of the presented system are automatically chosen with respect to the
filtering techniques have been presented based on the statistical estimated local SNR values. Accordingly, the image structures and
estimation theorem. Sijbers et al. [23,24] suggested to estimate the the noise characteristics will be effective on the filtering process; as a
Rician distribution parameters using a maximum likelihood (ML) result of which, a better denoising performance is obtained.

Fig. 1. Noise distribution in MRI. Bottom part shows a 3D MR volume of the head region corrupted by 10% of the Rician noise. Top part depicts the histograms of two highlighted
sub-volumes. As shown, the Rician noise tends to be Gaussian distributed when SNR is high (top, right side) and Rayleigh distributed when SNR is zero (top, left side).
1208 H.M. Golshan et al. / Magnetic Resonance Imaging 31 (2013) 1206–1217

The rest of this paper is organized as follows. Section 2 briefly underlying Gaussian noise, A is the noise-less signal value and M is
explains the noise characteristics in MRI. Section 3 elaborates the the observed noisy value.
proposed method. Section 4 presents the experimental results with Unlike the additive Gaussian noise, the Ricain noise is signal-
both synthetic and real MR data. Finally, Section 5 contains the dependent and its mean value depends on the local intensity in the
conclusions and remarks. image. When SNR is high, the Rician distribution asymptotically
tends to be Gaussian. By contrast, it becomes Rayleigh distributed
2. Noise characteristics in MRI when SNR is zero (that is, only noise is present, A → 0) [13,23].
Fig. 1 shows these different manners of the Rician noise for a
In MRI, the acquired raw data is complex valued and is typical MR volume.
represented in frequency domain (k-space). This data is intrinsically Finally, it should be noted that the Rician noise is especially
corrupted by the additive Gaussian noise due to the presence of problematic in the low SNR ranges where it introduces a signal-
different noise sources during the image formation [30]. After the dependent bias. Hence, many MRI denoising techniques have been
inverse Fourier transform, the MR data is still complex valued and developed to deal with the squared MR images where the Rician
the distribution of noise remains Gaussian because of the linearity random variables obey a non-central Chi-square distribution. As a
and the orthogonaity of the Fourier transform [24]. Accordingly, the result, the remained bias is a constant term (2σ 2) and can therefore
raw noisy MR data Y can be modeled as the summation of a noise- be easily removed [23]:
less data X and a Gaussian distributed noise N:
   
2 2 2
Y ¼ X þ N: ð1Þ E M ¼ E A þ 2σ : ð3Þ

It is common to transform the complex data into magnitude and where E(∙) denotes the expectation operator.
phase data since they are more directly related to the physiological and
anatomical quantities of interest [24]. After this nonlinear operation, the 3. Theory
probability density function of the MR data changes to Rician [23,24]:
!  3.1. Problem definition
2 2
M A þM

AM
PðM ¼ jY jjA; σ Þ ¼ 2 exp − I0 uðMÞ: ð2Þ
σ 2σ 2 σ 2
The presented filtering method in this paper is based on the
LMMSE estimator assuming a Rician noise model. The closed-form
where I0(∙) is the 0th order modified Bessel function of the first kind, analytical solution and the low computational cost are the main
u(∙) is the Heaviside step function, σ is the standard deviation of the advantages that make the LMMSE method a suitable estimation

Fig. 2. Illustration of the local estimation drawbacks. Left side shows two adjacent regions with different underlying gray levels (i.e. 50 for dark area and 150 for bright area)
corrupted by 10% of the Rician noise. A zoomed sub-region is also shown that provides the content of the samples to estimate a noise-less signal value. The histogram of each
region is depicted in the right side. As shown, due to the different underlying gray levels, the selected samples for the estimation process come from different Rician distributions.
As a result of which, undesirable estimation will be obtained.
 H.M. Golshan et al. / Magnetic Resonance Imaging 31 (2013) 1206–1217 1209

alternative for MRI denoising. In terms of the Rician distributed Table 1

random variables, the LMMSE estimator is given by Refs. [28,29]: Summary of the proposed method.

for each iteration, do:

^ D E Get the output of the previous iteration as the input.
2 2 2
A ¼ M −2σ Estimate the variance of noise using Eq. (11).
D E  1 Compute the local means and variances using sample estimators.
4σ 2 M 2 −σ 2
0
Compute the local SNR values using Eq. (7).
 D E
2 2
þ max@1−    2 ; 0A M − M : ð4Þ Compute the inverted mean values using Eq. (6).
M4 − M2
Compute the control parameters μ1 and σ12 using Eq. (9).
for each voxel do
if (SNR b 3)
where M is the observed noisy value, Â is the estimated value of the
Find the voxels within a search volume that comply with
original noise-less signal A and σ is the standard deviation of the the condition of ‘SNR b 3’ of Eq. (8).
underlying Gaussian noise. Finally, 〈·〉 denotes the sample estimator else
of the expectation. Find the voxels within a search volume that comply with
Eq. (4) shows that the estimation results depend mainly on the the condition of ‘otherwise’ of Eq. (8).
endif
statistical expectations of different orders. In practice, the quanti-
Take the selected voxels as the content of samples and
tative values of these expectations must be estimated using a set of compute the LMMSE estimation (Eq. (4)) by them.
samples selected from the noisy MR data. Therefore, the partic- endfor
ipated samples play an important role on the accuracy and Take the square root of the result.
endfor
robustness of the results.
Conventionally, the LMMSE method estimates the true signal
value from a local neighborhood within which the signal is assumed
to be a realization of a random variable. In accordance with this
are two threshold values that control the amount of closeness of the
principle, all the requirements of the filtering method are locally
calculated statistics.
estimated using a neighborhood centered on the voxel under
Substantially, the nonlocal LMMSE method presented in [31]
processing (e.g. a local neighborhood of 3 × 3 × 3 voxels). This
improves the estimation results by choosing the samples in a more
approach provides a straightforward implementation for the LMMSE
consistent manner; however, it suffers from some limitations mainly
estimator; nonetheless, it is not an efficient solution. A major
due to the nature of Eq. (5):
drawback appears on the high frequency components of the image
like edges and small details where the different underlying gray 1. The mean condition in Eq. (5) shows different manners in the
levels in a local neighborhood lead to the different realizations of the low and high values of the image gray levels. In other words,
Rician noise; as a consequence of which, undesirable estimation will the intensity-sensitive property of the mean estimator causes
be obtained (see Fig. 2 for more details). Furthermore, the 3D MR more strict limitations on the selection of samples in low gray
images intrinsically include many similar patterns that can be used levels than high.
to achieve better estimation results. In contrast to this fact, the 2. Using hard threshold values for the control parameters μ1 and
conventional LMMSE estimator doesn’t take advantage of high σ12 doesn’t seem to be an efficient idea while the local statistical
degree of redundancy in the contents of images because it limits the moments of the image vary with respect to the noise level.
samples to a local neighborhood only.
In an attempt to remedy the above-discussed problems, a
3.2. Proposed method
nonlocal version of the LMMSE estimator was presented in [31].
This method uses the voxels located in the nonlocal neighborhoods
As noted earlier, the performance of the LMMSE method directly
of a voxel under processing to estimate its true signal value.
relates to the samples participated in the estimation process.
However, the efficiency of this filtering scheme is highly sensitive
Therefore, an efficient selection method of samples can guarantee
to the selection methods of samples. In Ref. [31], this goal was
the accuracy of the results. In the current work, we present a suitable
accomplished using a similarity measure based on the local means
approach to face this particular issue which takes better into account
and variances. It was shown that these two statistical moments are
the statistical properties of the image structures. This approach is
suitable estimators to discriminate different structures in the
developed based on the given similarity measure in Eq. (5) and
images [19]. In [31], the selection of samples was performed using
provides a principled way to overcome the intrinsic drawbacks of
the following similarity measure:
this equation.
As discussed in Section 3. 1, the mean condition in Eq. (5) limits
< similar samples if μ1 ≤ Eðηx Þ ≤ 1 and σ12 ≤ Varðηx Þ≤ 1
8
the selection of samples in low gray levels due to its intensity-
>
E ηy μ1 Var ηy σ12 ð5Þ
> sensitive property. To resolve this problem, we propose to evaluate
dissimilar samples else
:

where E(ηx) and Var(ηx) represent respectively the mean and

Table 2
variance of the local neighborhood ηx and the parameters μ1and σ12 Method parameters: v is the radius of the search volume, f is the radius of the patches
(patch refers to a local neighborhood), h is the smoothing parameter used in the NLM-
based filters, σ is the standard deviation of the underlying Gaussian noise. Finally, μ1
and σ12 are the control parameters.

Methods Parameters

OBNLM [19] v = 5, f = 1, h = σ
WSM [20] v = 3, f1 = 1, f2 = 2, h = σ
ABONLM [22] v = 3, f1 = 1, f2 = 2
LMMSE [28] v = 5, f=1
Fig. 3. The presented decision-making functions for choosing the control parameters
RLMMSE [28] v = 5, f = 1, Iteration = 8
μ1 and σ12. The control parameters are automatically chosen with respect to the given
A1 [31] v = 5, f = 1,μ1 = 0.92, σ12 = 0.35
SNR value.
1210 H.M. Golshan et al. / Magnetic Resonance Imaging 31 (2013) 1206–1217

Fig. 4. Quality measures. Comparison of the quality measures with different image types and noise levels. Left shows the results of the NMSE measure; the lower values are the
best ones. Right shows the results of the SSIM index; the higher values are the best ones.
 H.M. Golshan et al. / Magnetic Resonance Imaging 31 (2013) 1206–1217 1211

the mean condition using the original and the inverted mean values, Eq. (6) only for the voxels that their corresponding SNR values are
where the inverted mean value is given by: lower than 3. Accordingly, Eq. (5) is modified as follows:

Inv ðE ðηx ÞÞ ¼ maxðM Þ−Eðηx Þ: ð6Þ 8 8 0 1

> ðSNRb3Þ @μ1 ≤ Eðηx Þ ≤ 1 or μ1 ≤ InvðEðηx ÞÞ 1A
>
> >
≤
> >
>
>
> > 
E ηy μ1 μ1
>
Inv E ηy
> >
where max(M) is the maximum value of the given MR data and
> >
>
>
> >
>
1
>
ð Þ
>
Var
Inv(E(ηx )) is the inverted mean value calculated in a local η
> <
and σ12 ≤  x ≤ 2
>
>
<
similar samples if > Var ηy σ1
neighborhood ηx of 3 × 3 × 3 voxels. >
>
>
>
E 1 1
>
ð Þ Var ðηx Þ
>
Eq. (6) intrinsically tends to improve the results in low gray
>
> > ηx 2
otherwise μ1 ≤   ≤ and σ1 ≤  ≤ 2
> >
>
>
> >
E
>
μ1 σ1
levels by increasing the number of similar samples, but it leads to
>
Var
>
>
>
>
>
: η y η y
>
useless calculations while dealing with high gray levels. Hence,
>
:
dissimilar samples else
we develop an SNR-adapted scheme in order to avoid additional ð8Þ
computations. SNR is theoretically defined as the ratio of the
noise-less signal A to the noise standard deviation σ (SNR = A/σ) Using the hard threshold values for the control parameters μ1 and
[9], however, it is not usable in practice due to the absence of A. σ12 is another drawback of the presented method in [31]. It was
Therefore, the SNR values are estimated here using the relation shown that both the control parameters depend on the noise level.
between noise and signal of the second order moment in a Ricain Through an exhaustive search on the synthetic noise-less MR images
distribution [23]: with different noise conditions, the pairs (μ1 = 0.9 and σ12 = 0.3)
and (μ1 = 0.95 and σ12 = 0.5) were empirically found to be the best
    options for low and high SNR values respectively. In practical cases
E M2 ¼ E A2 þ 2σ 2
g⇒
vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
!ffi
E M2
u
u where the SNR values are not clearly available, independently of the
A SNR ¼ max
t −2; 0 : ð7Þ
SNR ¼ σ2 noise level the pair (μ1 = 0.92 and σ12 = 0.35) has been suggested
σ as a general alternative. Such a method doesn’t provide any
perspective for an automatic scheme and proposes to use hard
where the expectation value is calculated locally using a neighbor- threshold values for the control parameters, which is a sub-optimal
hood of 3 × 3 × 3 voxels. way to deal with this issue.
Eq. (7) assigns a specific SNR value to each voxel of the MR Here, we present an SNR-adapted approach which is able to set
volume. In order to avoid useless computations, it is suggested to use the control parameters automatically with respect to the estimated

Original Image Noisy Image

Filtered Images
WSM ABONLM SNLMMSE RSNLMMSE

Difference Images
WSM ABONLM SNLMMSE RSNLMMSE

Fig. 5. Example denoising results for an axial slice of the synthetic T1w MR volume (Rician noise level of 10%).
1212 H.M. Golshan et al. / Magnetic Resonance Imaging 31 (2013) 1206–1217

local SNR values. Consequently, the image structures and the noise volume is used as the noisy input of the current iteration. In this
characteristics are taken into account when choosing the samples. To way, the filter should reach a steady state when the estimated noise
present our method, the estimated local SNR values are divided into declines to zero:
three ranges. They are SNR b 5, 5 ≤ SNR ≤ 10 and SNR N 10 which
^
respectively provide the low, medium and high cases. In order to 2
D E
2 2
reduce the complexity and without loss of generality, in the medium A ¼ E M −2σ
D E  1
2 2 2
0
range of SNR, we define a linear relation between the control 4σ M −σ 
2
D E ∧
2
parameters and the estimated local SNR values. In the low and high þ max 1−    2 ; 0A M −E M
@ ⇒ A ¼ M: ð10Þ
M − M
4 2
ranges, it was proved that the reported values in Ref. [31] are
σ→0
approximately the best choices and therefore, they are used here
with no changes. Accordingly, the proposed decision-making In the case of the proposed method, it was empirically found that
functions are given by: three iterations are large enough to attain the steady state.
The performance of the recursive filtering method exactly relates
to the capability of properly estimate the noise level. From the

f
0:9 SNRb5
SNR þ 85 available approaches to estimate the noise in magnitude MR images,
μ1 ¼ 5≤SNR≤10; ð9Þ
100 e.g. Refs. [28] and [33–35], we used the one proposed in Ref. [28] due
0:95 SNRN10 to its consistency with the recursive filtering scheme. This method
estimates the noise using the local variances in a sub-volume inside
2
σ1 ¼ f 0:3
4SNR þ 10
0:5
100
SNRb5
5≤SNR≤10
SNRN10
the signal region where the Rician assumption asymptotically holds
once the image is filtered [28]. The noise estimator is given by:
 
2 2
σ ¼ mode σx : ð11Þ

Fig. 3 depicts the decision-making functions graphically. As where σx 2 is the unbiased local sample variance that may be
shown, they set the control parameters with respect to the calculated locally using a 3 × 3 × 3 neighborhood within a seg-
corresponding SNR values. mented region with signal and no background. Finally, the mode of
The LMMSE filtering method can be performed recursively if the all the local estimations is used as the noise variance.
variance of noise is dynamically estimated in each iteration [28,29]. The proposed selection method of samples originally tends to
Indeed, by the recursive implementation, the previous filtered evaluate all the voxels in the volume, but this manner is

Original Image Noisy Image

Filtered Images
WSM ABONLM SNLMMSE RSNLMMSE

Difference Images
WSM ABONLM SNLMMSE RSNLMMSE

Fig. 6. Example denoising results for an axial slice of the synthetic T2w MR volume (Rician noise level of 10%).
 H.M. Golshan et al. / Magnetic Resonance Imaging 31 (2013) 1206–1217 1213

computationally inefficient. Hence, we used a search volume of The second quality measure is the structural similarity (SSIM)
11 × 11 × 11 voxels for each voxel under processing. As shown in index [37], which is a criterion more consistent with the human
Refs. [19,32], this choice leads to a good compromise between visual system (HVS). SSIM reveals the structural and perceptual
computational cost and quality of the results. similarity between two images and hence, it is a commonly used
At the following sections we refer to our SNR-adapted Nonlocal quality measure in the context of MRI literature. It is given by:
LMMSE method as SNLMMSE and to its Recursive version as
  
RSNLMMSE. Finally, the routine of the proposed method is presented N 2μx μy þ c1 2σxy þ c2
1 X
 
∧
in Table 1. SSIM A; A ¼ ð13Þ
N x;y¼1 μ 2 þ μ 2 þ c
  :
x σ2 þ σ2 þ c
y 1 x y 2
4. Experiments and results
where μx and μy are the local mean value of images A and Â, σx and σy
4.1. Synthetic materials and quality measures are the respective standard deviations, σxy is the covariance value
and c1 and c2 are two constants. As suggested in Ref. [37], the SSIM is
To compare the filtered results with a reference dataset, we used locally estimated using a Gaussian kernel of 3 × 3 × 3 voxels. Finally,
the Brainweb [36] standard 3D MRI phantoms, T1-weighted (T1w), the mean value of all the local estimations is used as a quality metric.
T2w and PDw volumes of 181 × 217 × 181 voxels with 1 mm 3 voxel
resolution. We employed the 8 bit precision data where the original 4.2. Validation on synthetic data
values were in the range [0,255].
To evaluate the restoration performance of methods, two quality This section is devoted to compare the introduced approach with
measures are used. The first is the normalized mean squared error some recently proposed methods used for MRI denoising. For the
(NMSE) metric, which computes the distance between two images A sake of brevity, the competitive methods and their optimum
and  of the same size N as follows: parameters suggested by authors are tabulated in Table 2.
The experimental results are obtained using the Brainweb
N 
X ∧ 2

synthetic datasets. In order to evaluate the stability of methods

1
 N Ax − Ax with different noise conditions, a wide range of the Rician noise (5–
∧
NMSE A; A ¼ x¼1 : ð12Þ 20% of maximum intensity with 5% in step) were considered in the
μ μ∧
A A experiments. The quantitative results with different image types and
noise levels are depicted in Fig. 4. As shown, our proposed SNLMMSE
where μ A and μ denote the mean value of images A and and RSNLMMSE outperform the other compared approaches in most
 respectively. situations, apart from two cases with the T2w volume and the noise

Original Image Noisy Image

Filtered Images
WSM ABONLM SNLMMSE RSNLMMSE

Difference Images
WSM ABONLM SNLMMSE RSNLMMSE

Fig. 7. Example denoising results for an axial slice of the synthetic PDw MR volume (Rician noise level of 10%).
1214 H.M. Golshan et al. / Magnetic Resonance Imaging 31 (2013) 1206–1217

Original Images Noisy Images

Filtered Images Difference Images

OBNLM

WSM

ABONLM

SNLMMSE

RSNLMMSE

Fig. 8. Qualitative comparison of methods under consideration on a sample pathological case. The first row shows three continuous slices of the MR volume including MS lesions.
In order to facilitate the visual analysis, the other rows only show magnification of the white square region. The filtered results for noise level of 15% are provided in each case. As
can be observed, the proposed SNLMMSE and RSNLMMSE show better performance in preserving the visual signature of the given pathology and removing the noise.
 H.M. Golshan et al. / Magnetic Resonance Imaging 31 (2013) 1206–1217 1215

Table 3
Quantitative comparison of methods under consideration on a brain T2w MR volume with MS lesions.

5% 10% 15%

NMSE SSIM BC NMSE SSIM BC NMSE SSIM BC

Noisy data 0.00850 0.5738 0.6949 0.0335 0.3584 0.6629 0.0720 0.2653 0.6382
OBNLM 0.00093 0.9612 0.8653 0.0026 0.8888 0.7888 0.0047 0.7967 0.7420
WSM 0.00085 0.9634 0.8526 0.0024 0.8929 0.7830 0.0047 0.7973 0.7369
ABONLM 0.00086 0.9611 0.8509 0.0023 0.8985 0.8006 0.0041 0.8153 0.7613
LMMSE 0.00200 0.8936 0.7974 0.0023 0.7822 0.7452 0.0081 0.6760 0.7103
RLMMSE 0.00180 0.9304 0.8799 0.0046 0.8907 0.8303 0.0052 0.8552 0.7903
A1 0.00100 0.9629 0.9001 0.0035 0.9091 0.8078 0.0041 0.8087 0.7545
SNLMMSE 0.00089 0.9652 0.8938 0.0022 0.9189 0.8484 0.0036 0.8638 0.8198
RSNLMMSE 0.00088 0.9676 0.9109 0.0021 0.9320 0.9028 0.0032 0.8959 0.8847

The quality measures have been computed in a region of interest (ROI) including MS lesions. The results are shown for different levels of the Rician noise. Best value of each
column is highlighted.

levels of 5% and 10% where the OBNLM and WSM slightly give better between tissues and showing fewer oscillations over homogeneous
results for the NMSE metric. In terms of the SSIM index, however, the areas are the main advantages of the proposed method.
proposed method obtains the best results in all the cases.
In order to facilitate the reproducibility of the results, the MATLAB 4.3. Validation on pathological cases
(Mathworks Inc.) source code of the presented experiments will be
available on the following Web page: [Link] The anatomical structures and small features are the most
Figs. 5–7 have been provided for a visual comparison of the valuable parts of the medical images to the diagnosis. Hence, a
methods under consideration. Also, the absolute difference between suitable denoising method should retain theses image contents
the original and denoised images is shown in these figures to indicate while removing the noise. To evaluate the methods in this situation,
the error produced by the filtering process (a difference image with we used a sample MR data with MS (multi sclerosis) lesions. MS is a
lower gray levels is more desirable). As shown, all the compared neurological disease that appears in the brain MR images in the form
methods have removed the noise signature successfully, however, of small plates. From a pathological point of view, the number of
more preserved fine structures like vessels, better defined contrast these lesions is an important parameter to follow the progress of

Original Image SNLMMSE RSNLMMSE

Fig. 9. Example results of the proposed method applied to real cervical T1w MR images. A typical sagittal slice is shown in each case. The middle row shows magnification of the
white square region. Bottom row shows the corresponding residual images.
1216 H.M. Golshan et al. / Magnetic Resonance Imaging 31 (2013) 1206–1217

disease. Moreover, their contrast on the image is another alternative The quality measures obtained by the methods under consider-
that helps the physicians for a more accurate diagnosis. ation with different levels of the Rician noise are tabulated in Table 3.
In the experiments, we used a T2w MRI phantom including MS As shown, the proposed RSNLMMSE achieves the best results in most
lesions from Brainweb [36]. This data was a 3D volume of situations and the SNLMMSE shows very competitive performance.
181 × 217 × 181 voxels and 1 mm 3 voxel resolution and its original The BC values given by the proposed method show much
values were in the range [0,255] (8 bit precision data). improvement over the other filters. This ensures the contrast
Different levels of the Rician noise were evaluated while between different tissues is visually better preserved on the filtered
comparing the methods. Fig. 8 shows the filtered results with the images using the introduced approach.
Rician noise level of 15%. As can be observed, showing a more
smoothed background region with fewer trace of artifacts make the
MS lesions visually more contrasted with the proposed method. The 4.4. Validation on real MR data
NLM-based methods (OBNLM, WSM and ABONLM) caused blurred
results and failed to retain some of the small MS lesions, probably To verify the consistency of the proposed method with real MR
because of the excessive smoothing over the homogeneous areas. images, the experiments were carried out on two datasets. The first
In order to compare the methods quantitatively, we used a sub- one was a T1w cervical 3D MR volume acquired on a 1 Tesla intera
volume of the original MR data including the MS lesions. It was taken Philips scanner (TR = 400 ms, TE = 11 ms, flip angle = 90°, voxel
to be the region of interest (ROI) in our assessments. In the res olut i on = 0 .5 957 × 0. 59 57 × 3 mm 3 , vol um e s i z e =
experiments, besides the NMSE and SSIM, we also used the 512 × 512 × 9 voxels). The filtering results for this dataset can be
Bhattacharrya coefficient (BC) as an alternative quality measure observed in Fig. 9. It should be noted that a strong filtering method
because of its relevance to the image contrast [13]. The BC measures should extract as much noise as possible from the image while
the closeness of two discrete probability distributions p and q over a keeping the structural features unaltered. Here, the residual images
same domain X (e.g. X∈ [0,255] for the 8 bit precision data) as follows: (i.e. the difference between the noisy and filtered images) are used
as a visual criterion to follow the traces of useful anatomical
pðxÞqðxÞ: information removed by the filtering process. As shown in Fig. 9,
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
BC ðp; q Þ ¼ ∑ ð14Þ
x∈X both the proposed SNLMMSE and RSNLMMSE have removed the
noise successfully. Besides a good definition of the contrast between
The range of the coefficient is BC∈ [0,1] and the close value to 1 tissues, the filtered images also show a relevant performance in
shows that the data samples corresponding to the distributions p and preserving the detailed structures. Additionally, there is no trace of
q are very similar. the significant anatomical information on the residual images.

Original Image SNLMMSE RSNLMMSE

Fig. 10. Example results of the proposed method applied to the real brain T2w MR images. From top to bottom: typical axial slices of the original and filtered volumes,
magnification of the white square region, and the corresponding residual images.
 H.M. Golshan et al. / Magnetic Resonance Imaging 31 (2013) 1206–1217 1217

The second MR dataset was also obtained using the above- [7] Basu S, Fletcher T, Whitaker R. Rician noise removal in diffusion tensor MRI. Proc
MICCAI 2006:117–25.
described scanner. It was a T2w brain 3D MR volume (TR = [8] Krissian K, Aja-Fernandez S. Noise-driven anisotropic diffusion filtering of MRI.
2095 ms, TE = 110 ms, flip angle = 90°, voxel resolution = IEEE Trans Image Process 2009;18(10):2265–74.
0.8984 × 0.8984 × 5 mm 3, volume size = 256 × 256 × 18voxels). [9] Nowak RD. Wavelet-based Rician noise removal for magnetic resonance
imaging. IEEE Trans Image Process 1999;8(10):1408–19.
As shown in Fig. 10, the proposed method attains good results in the [10] Pizurica A, Philips W, Lemahieu I, Acheroy M. A versatile wavelet domain noise
terms of noise removal and contrast. According to the residual filtration technique for medical imaging. IEEE Trans Med Imaging 2003;22:
images, it can be appreciated that no anatomical information has 323–31.
[11] Bao P, Zhang L. Noise reduction for magnetic resonance images via adaptive
been lost due to the denoising process. multiscale products thresholding. IEEE Trans Med Imaging 2003;22(9):1089–99.
[12] Wood JC, Johnson KM. Wavelet packet denoising of magnetic resonance images:
5. Conclusion and discussion importance of Rician noise at low SNR. Magn Reson Med 1999;41:631–5.
[13] Anand CS, Sahambi JS. Wavelet domain non-linear filtering for MRI denoising.
Magn Reson Imaging 2010;28:842–61.
We have presented a nonlocal version of the LMMSE estimator [14] Muresan DD, Parks TW. Adaptive principal components and image denoising.
for Rician noise removal in the 3D MRI. The proposed method takes IEEE Int Conf Image Process 2003;1:101–4.
[15] Guleryuz OG. Weighted averaging for denoising with overcomplete dictionaries.
advantage of the fact of image data redundancy as an intrinsic
IEEE Trans Image Process 2007;16:3020–34.
property of MR images. We modeled the MR data as random fields [16] Manjon JV, Coupe P, Baudes A, Louis Collins D, Robels M. New methods for MRI
and developed a similarity measure to choose the appropriate subset denoising based on sparseness and self-similarity. Med Imaging Anal
2012;16(1):18–27.
of samples within a large portion of the given data. It is in contrast to
[17] Buades A, Coll B, Morel JM. A review of image denoising algorithms, with a new
the conventional LMMSE method which limits the selection of one Multiscale Model. Simulation 2005;4(2):490–530.
samples to a local neighborhood only. [18] Wiest-Daessle N, Prima S, Coupe P, Morrissey SP, Barillot C. Rician noise removal
All the required parameters of the presented method are by non-local means filtering for low signal-to-noise ratio MRI: applications to
DT-MRI. Proc MICCAI 2008:171–9.
automatically derived with respect to the estimated local SNR values. [19] Coupe P, Yger P, Prima S, Hellier P, Kervrann C, Barillot C. An optimized blockwise
In this way, the local structures of image and the characteristics of nonlocal means denoising filter for 3-D magnetic resonance images. IEEE Trans
noise participate in the selection method of samples. As a conse- Med Imaging 2008;27(4):425–41.
[20] P. Coupe, P. Hellier, S. Prima, C. Kervrann, and C. Barillot, “3D wavelet subbands
quence, the estimation results improve noticeably. mixing for image denoising,” Int. J. Biomed. Imaging., [Link]
Various validations have been performed with both synthetic and 2008/590183 (Article ID: 590183).
real MR images. Quantitative analysis based on the different [21] Manjon JV, Carbonell-Caballero J, Lull JJ, Garcia-Marti G, Marti-Bonmati L, Robles
M. MRI denoising using non-local means. Med Imaging Anal 2008;12:514–23.
measures, NMSE, SSIM and BC demonstrates that the presented [22] Manjon JV, Coupe P, Marti-Bonmati L, Collins DL, Robles M. Adaptive nonlocal
method outperforms the other compared filters in most situations. In means denoising of MR images with spatially varying noise levels. J Magn Reson
particular, preservation of the fine structural details while removing Imaging 2010;31(1):192–203.
[23] Sijbers J, den Dekker AJ, Scheunders P, Van Dyck D. Maximum-likelihood
the noise makes the proposed method a confident denoising
estimation of Rician distribution parameters. IEEE Trans Image Process
alternative for MR images. To evaluate this issue, we focused on 1998;17(3):357–61.
pathological cases (T2w MR images with MS lesions). Due to the [24] Sijbers J, den Dekker AJ. Maximum likelihood estimation of signal amplitude and
noise variance from MR data. Magn Reson Med 2004;51(3):586–94.
experimental results, the proposed method showed improvements
[25] He L, Greenshields IR. A nonlocal maximum likelihood estimation method for
over the previous methods both in preserving the given pathology Rician noise reduction in MR images. IEEE Trans Med Imaging 2009;28(2):
and removing the noise. 165–72.
Finally, a future research topic includes the design of new similarity [26] Rajan J, Jeurissen B, Verhoye M, Audekerke JV, Sijbers J. Maximum likelihood
estimation-based denoising of magnetic resonance images using restricted local
measures. Another extension to this work would be to develop a more neighborhoods. Phys Med Biol 2011;56:5221–34.
efficient method for the selection of filtering parameters. [27] Wong A, Mishra AK. Quasi-Monte Carlo estimation approach for denoising
MRI data based on regional statistics. IEEE Trans Biomed Eng 2011;58(4):
1076–83.
Acknowledgments [28] Aja-Fernandez S, Alberola-Lopez C, Westin C-F. Noise and signal estimation in
magnitude MRI and Rician distributed images: A LMMSE approach. IEEE Trans
The authors are grateful to Dr. Alizadeh for his useful comments Image Process Aug. 2008;17(8):1383–98.
[29] Aja-Fernandez S, Niethammer M, Kubicki M, Shenton ME, Westin C-F.
on the diagnostic details of the denoised MR images. Also, the authors Restoration of DWI data using a Rician LMMSE estimator. IEEE Trans Med
would like to thank the R&D center of Poursina Hospital, Rasht, Iran Imaging 2008;27(10):1389–403.
for providing access to the real MR data used in this research. [30] Edelstein WA, Glover GH, Hardy CJ, Redington RW. The intrinsic SNR in NMR
imaging. Magn Reson Med 1986;3(4):604–18.
[31] Golshan HM, Hasanzadeh RPR. “A non-local Rician noise reduction approach for
References 3-D magnitude magnetic resonance images,” in Machine Vision and Image
Processing (MVIP), 2011 7th Iranian, pp. 1–5, 2011.
[1] Wright GA. “Magnetic resonance imaging,” IEEE signal processing magazine, 1. [32] Golshan HM, Hasanzadeh RPR. “A modified Rician LMMSE estimator for the
p. 56–66. restoration of magnitude MR images”, Int. J. Light Electron Opt. (Optik), http://
[2] Gerig G, Kubler O, Kikinis R, Jolesz FA. Nonlinear anisotropic filtering of MRI data. [Link]/10.1016/[Link].2012.07.001.
IEEE Trans Med Imaging 1992;11(2):221–32. [33] Sijbers J, Poot D, den Dekker AJ, Pintjenst W. Automatic estimation of the noise
[3] Perona P, Malik J. Scale-space and edge detection using anisotropic diffusion. variance from the histogram of a magnetic resonance image. Phys Med Biol Feb.
IEEE Trans Pattern Anal Mach Intell 1990;12(7):629–39. 2007;52:1335–48.
[4] Sijbers J, den Dekker AJ, Van der Linden A, Verhoye M, Van Dyck D. Adaptive [34] Coupe P, Manjon JV, Gedamu E, Arnold D, Robles M, Collins DL. Robust Rician
anisotropic noise filtering for magnitude MR data. Mag Reson Imaging noise estimation for MR images. Med Imaging Anal 2010;14:483–93.
1999;17(10):1533–9. [35] Rajan J, Poot D, Juntu J, Sijbers J. “Noise measurement from magnitude MRI using
[5] Samsonov AA, Johnson CR. Noise-adaptive nonlinear diffusion filtering of MR local estimates of variance and skewness,. Phys Med Biol 2010;55:441–9.
Images with spatially varying noise levels. Magn Reson Med 2004;52:798–806. [36] [Link]
[6] Lysaker M, Lundervold A, Xue-Cheng T. Noise removal using fourth-order partial [37] Wang Z, Bovik AC, Sheikh HR, Simonceli EP. Image quality assessment: from
differential equation with applications to medical magnetic resonance images in error visibility to structural similarity. IEEE Trans Image Process 2004;13(4):
space and time. IEEE Trans Image Process 2003;12:1579–90. 600–12.