PDF (6 MB) - Magnetic Resonance Imaging

Magnetic Resonance Imaging 32 (2014) 702–720
Contents lists available at ScienceDirect
Magnetic Resonance Imaging
journal homepage: www.mrijournal.com
Generalized total variation-based MRI Rician denoising model with
spatially adaptive regularization parameters
Ryan Wen Liu a, b, Lin Shi b, c,⁎, Wenhua Huang d,⁎⁎, Jing Xu d, Simon Chun Ho Yu a, Defeng Wang a, b, e, f
a
Department of Imaging and Interventional Radiology, The Chinese University of Hong Kong, Shatin, New Territories, Hong Kong SAR, P.R. China
Research Center for Medical Image Computing, The Chinese University of Hong Kong Shatin, New Territories, Hong Kong SAR, P.R. China
Department of Medicine and Therapeutics, The Chinese University of Hong Kong, Shatin, New Territories, Hong Kong SAR, P.R. China
d
Institute of Clinical Anatomy, Southern Medical University, Guangzhou, Guangdong, P.R. China
e
CUHK Shenzhen Research Institute, Shenzhen, Guangdong, P.R. China
f
Department of Biomedical Engineering and Shun Hing Institute of Advanced Engineering, The Chinese University of Hong Kong, Shatin, New Territories, Hong Kong SAR, P.R. China
b
c
a r t i c l e
i n f o
Article history:
Received 28 June 2013
Revised 13 January 2014
Accepted 7 March 2014
Keywords:
Total variation
Magnetic resonance imaging (MRI)
Diffusion tensor MRI (DT-MRI)
Image denoising
Hyper-Laplacian prior
Rician distribution
a b s t r a c t
Magnetic resonance imaging (MRI) is an outstanding medical imaging modality but the quality often
suffers from noise pollution during image acquisition and transmission. The purpose of this study is to
enhance image quality using feature-preserving denoising method. In current literature, most existing MRI
denoising methods did not simultaneously take the global image prior and local image features into
account. The denoising method proposed in this paper is implemented based on an assumption of spatially
varying Rician noise map. A two-step wavelet-domain estimation method is developed to extract the noise
map. Following a Bayesian modeling approach, a generalized total variation-based MRI denoising model is
proposed based on global hyper-Laplacian prior and Rician noise assumption. The proposed model has the
properties of backward diffusion in local normal directions and forward diffusion in local tangent
directions. To further improve the denoising performance, a local variance estimator-based method is
introduced to calculate the spatially adaptive regularization parameters related to local image features and
spatially varying noise map. The main benefit of the proposed method is that it takes full advantage of the
global MR image prior and local image features. Numerous experiments have been conducted on both
synthetic and real MR data sets to compare our proposed model with some state-of-the-art denoising
methods. The experimental results have demonstrated the superior performance of our proposed model in
terms of quantitative and qualitative image quality evaluations.
© 2014 Elsevier Inc. All rights reserved.
1. Introduction
1.1. Background and related work
Magnetic resonance imaging (MRI) is an outstanding medical
imaging modality but often suffers from noise pollution during
image acquisition and transmission. For single-coil MR acquisitions,
the magnitude MR images are usually modeled using a Rician
distribution. If multiple-coil acquisitions and no subsampling of the
k-space raw data are considered, the magnitude data will follow a
non-central Chi (nc-χ) distribution [1–3]. MRI denoising has
attracted tremendous attentions of researchers in the field of
biomedical imaging. Compared with non-medical imaging applications,
⁎ Correspondence to: Lin Shi, Department of Medicine and Therapeutics, The Chinese
University of Hong Kong, Shatin, New Territories, Hong Kong SAR, P.R. China. Tel.: +852
26322975; fax: +852 26487269.
⁎⁎ Correspondence to: Wenhua Huang, Institute of Clinical Anatomy, Southern Medical
University, Guangzhou, Guangdong, P.R. China. Tel.: +86 20 61648183.
E-mail addresses: [email protected] (L. Shi), [email protected]fimmu.com (W. Huang).
http://dx.doi.org/10.1016/j.mri.2014.03.004
0730-725X/© 2014 Elsevier Inc. All rights reserved.
MRI denoising methods should put more emphasis on preservation of
fine structures during noise reduction. In clinical settings, the fine
structures contain important medical information, which could assist
physicians with accurate diagnosis.
There are a variety of study methods that have been developed for
MRI denoising. Henkelman [4] presented the first attempt to estimate
the noiseless magnitude MR image from its noisy version. Recent
methods on MRI denoising have employed nonparametric statistical
techniques. For example, Awate and Whitaker [5] proposed a novel
MRI denoising method with nonparametric neighborhood statistics.
They also introduced a feature-preserving denoising approach by
inferring uncorrupted-signal Markov statistics as a prior under an
empirical Bayesian framework [6]. More recently, López-Rubio and
Florentín-Núñez [7] considered a nonparametric regression method
depending on a zeroth order three-dimensional (3D) kernel regression, which computed a weighted average of pixels over a regression
window. A novel 3D image denoising procedure was developed using
local smoothing and nonparametric regression [8]. This method could
preserve edges and major edge structures in the restored MR images.
R.W. Liu et al. / Magnetic Resonance Imaging 32 (2014) 702–720
Wavelet transform-based approaches have also been considered in
the literature [9,10]. The image in wavelet domain is decomposed into
a set of multiresolutional wavelet coefficients. The detail coefficients
are processed with hard or soft thresholding at some particular levels
to effectively remove noise from an image [11]. Wood and Johnson
[12] have used wavelet packets for Rician noise reduction in low
signal-to-noise ratio (SNR) MR images. Recently, a hybrid denoising
method which combined the 3D optimized blockwise nonlocal-means
(NLM) filter [13] with discrete wavelet transform was proposed based
on a multiresolution approach [14]. Anand and Sahambi [15] presented
a wavelet-based bilateral filtering scheme which has achieved
remarkable denoising performance. Besides, many other transforms
have also been successfully applied to reduce Rician noise. For instance,
Manjón et al. [16] introduced a new noise reduction method by using
an overcomplete local principal component analysis (PCA)-based
decomposition. Two-dimensional (2D) principal component decomposition and Haar transformation were incorporated into the
structure-adaptive sparse denoising scheme [17]. The discrete cosine
transform (DCT)-based MRI denoising methods have also attracted
much attention in recent years [18,19].
Other denoising methods have been presented based on partial
differential equations (PDE). The anisotropic diffusion (AD) filter was
first introduced by Gerig et al. [20] to MRI in 1992. By Basu et al. [21] the
Rician noise model was incorporated into the Perona-Malik filter [22]
for denoising of diffusion tensor MRI (DT-MRI) in the maximum a
posteriori (MAP) framework. However, the AD-based methods were
highly dependent on some parameters, such as edge enhancement
parameter, conductance parameter and iteration stopping time [23].
To overcome this disadvantage, Tong et al. [24] presented an automatic
parameter selection scheme for AD filter. In 2009, Krissian and AjaFernández [25] introduced an extension of the speckle reducing
anisotropic diffusion (SRAD) [26] to remove Rician noise. More
recently, Golshan et al. [27] proposed an effective denoising method
based on SNR-adapted nonlocal linear minimum mean square error
(SNLMMSE). A recursive version (RSNLMMSE) of SNLMMSE was also
addressed to further enhance denoising performance.
Following Buades et al. [28], the NLM-based denoising methods
have received considerable attention during recent years. In the
NLM-based methods, each noiseless pixel value is estimated by the
weighted average of pixels with related surrounding neighborhoods
[29]. But the methods may result in ineffective MRI denoising,
because they generate satisfactory results only under the assumption of Gaussian random noise. Based on the second-order moment
of a Rician distribution, the unbiased NLM (UNLM) was extended to
reduce Rician-distributed noise present in low-SNR MRI [30,31].
Furthermore, Coupé et al. [13] proposed an effective denoising
method based on a 3D optimized blockwise version of the NLM.
From a statistical point of view, nonlocal maximum likelihood
(NLML) estimation methods have been considered [32,33]. More
recently, Manjón et al. [18] presented two novel 3D MRI denoising
techniques based on sparseness and self-similarity, i.e., oracle-based
DCT filter (ODCT3D) and prefiltered rotationally invariant NLM filter
(PRI-NLM3D). MR images in practice tend to be influenced by
spatially varying Rician noise [34], the above-mentioned methods
easily lead to unsatisfactory denoising results. To overcome this
limitation, Manjón et al. [35] have developed an adaptive NLM-based
method to deal with the inhomogeneous distribution of noise.
Apart from the aforementioned methods, several denoising
methods have been considered based on total variation (TV)
regularizer. The TV-based method, first proposed by Rudin et al.
[36], is capable of preserving edges and removing image noise in
homogeneous regions. Drapaca [37] has investigated the effects of
regularization parameters on TV-based MRI denoising results. Wang
and Zhou [38] proposed a hybrid MRI denoising method by
combining the TV scheme with wavelet transform. However, these
703
TV-based models were presented based on the noise assumption of
Gaussian distribution, which was different from the Rician distribution present in degraded MR images. Following a Bayesian modeling
approach, the Rician-distributed-based TV models have been
considered for MRI denoising [39,40]. More recently, Riciandistributed-based second-order total generalized variation (TGV)
model and its modified versions have also been proposed to enhance
medical image quality [41,42]. We mainly focus on the TV-based MRI
denoising method in this paper.
1.2. Motivation and contributions
It is well known that the regularization parameter plays a critical
role in TV-based noise reduction. To achieve good denoising results,
the parameter should be calculated adaptively according to local image
features. In current literature, the optimal selection of regularization
parameter has attracted increasing amounts of attention for TV-based
additive/multiplicative noise removal. As discussed in Refs. [43,44], the
regularization parameter was inversely proportional to the scale of
image feature. Thus the parameters localized at image features of
different scales should be selected differently. In particular, the
parameters are large in detail regions with small-scale features and
small in homogeneous regions with large-scale features. Gilboa et al.
[45] developed an adaptive texture-preserving denoising method
which controlled the level of denoising by computing the spatially
varying constrains based on local variance measures. Dong et al. [44]
proposed a multi-scale TV-based image restoration model with an
automatic selection of spatially adaptive regularization parameter
scheme to enhance image details while removing noise in homogeneous regions. These TV-based methods have gained a lot of attention
due to their ability to effectively reduce additive Gaussian noise. To
remove multiplicative gamma and Poisson noise, Li et al. [46,47] and
Chen and Cheng [48] have also developed local feature-based methods
to estimate the adaptive regularization parameters to generate good
denoising results.
However, the magnitude MRI signal in practice is corrupted by a
Rician or nc-χ noise, not a Gaussian, gamma or Poisson one. Thus the
aforementioned TV-based methods are unable to effectively reduce the
noise present in magnitude MR images. There is a huge potential to
develop an automatic method to select the spatially adaptive
regularization parameter scheme for TV-based MRI denoising models
[39,40]. On the other hand, most existing denoising methods are based
on an assumption of uniform distribution of noise standard deviation
(NSD) across an image. In contrast, our proposed denoising method
will be implemented based on a more natural assumption of spatially
varying Rician noise map (i.e., spatially varying distributed NSD). In this
paper, a generalized total variation (GTV)-based MRI Rician denoising
method is proposed based on global hyper-Laplacian image prior and
Rician noise assumption. In what follows, a local variance estimatorbased method is developed to calculate the spatially adaptive
regularization parameters related to local image features and spatially
varying noise map. Our proposed MRI denoising framework significantly differs from previous works in the following aspects:
1. The denoising method proposed in this paper is implemented
based on an assumption of spatially varying Rician noise map.
Before doing noise reduction, a two-step wavelet-domain
estimation method is introduced to extract the spatially varying
noise map for each voxel individually.
2. Following a Bayesian modeling approach, a GTV-based MRI
denoising model is proposed based on global hyper-Laplacian
image prior and Rician noise assumption. The proposed model has
the properties of backward diffusion in local normal directions
and forward diffusion in local tangent directions to perform
feature-preserving denoising.
704
R.W. Liu et al. / Magnetic Resonance Imaging 32 (2014) 702–720
Fig. 1. Characteristics of Rician noise in MRI. From top-left to top-right: a 2D T1-weighted MR image corrupted by 15% of Rician noise and a mask image separating the
background, low-SNR and high-SNR regions. As shown, the Rician noise tends to be Gaussian distributed when SNR is high (bottom-left) and Rayleigh distributed when SNR
closes to zero (bottom-right). In low-SNR regions, the Rician noise tends to be neither Gaussian nor Rayleigh distributed (bottom-middle).
3. To further improve the denoising performance, a local variance
estimator-based method is introduced to calculate the spatially
adaptive regularization parameters according to local image
features and spatially varying noise map. In particular, the large
regularization parameters in texture regions result in fine detail
preservation; whereas small ones in homogeneous regions
perform well in noise reduction.
and estimation of spatially varying noise map in MRI. In Section 3, we
present the GTV-based MRI Rician denoising model with spatially
adaptive regularization parameters. Numerous experiments on both
synthetic and real MR data sets are performed in Section 4. Finally
we conclude this paper by summarizing our contributions and
discussing the future work in Section 5.
2. Noise estimation in magnitude MR images
The main benefit of our proposed method is that it takes full
advantage of the global MR image prior and local image features.
Thus it can effectively reduce noise levels while preserving edges and
fine scale details in practice.
1.3. Organization
The remainder of this paper is organized into several
sections. Section 2 briefly explains the Rician noise distribution
2.1. Noise characteristics in MRI
In MRI the measured signal magnitude consists of real and
imaginary parts. If both real and imaginary parts are corrupted
by the zero-mean uncorrelated Gaussian noise with equal
variance, noise in magnitude MRI data can no longer be modeled
by the Gaussian distribution via the complex Fourier transformation [34]. As mentioned in Introduction, the magnitude
Fig. 2. Left: graphical representation of the function F(θ) with different magnitudes of 〈M〉/σ. Right: the correction factor ξ(θ) as a function of SNR = θ.
R.W. Liu et al. / Magnetic Resonance Imaging 32 (2014) 702–720
signal from a multiple-coil MR acquisition system follows a nc-χ
distribution [49]
!
n−1
2
2
M
M
M þ A
MA
In−1
pðMjA; σ Þ ¼ 2 exp −
;
ð1Þ
2
2
A
σ
σ
2σ
where M is the value in the magnitude image, A is the amplitude of
the noiseless signal, σ 2 denotes the variance of Gaussian noise in
complex domain, n is the number of channels used for parallel data
acquisition, and In − 1(⋅) represents the (n − 1)th-order modified
Bessel function of the first kind. Without loss of generality, we only
consider the case of n = 1 in this paper, then (1) reduces to the
Rician distribution [50], i.e.,
! 2
2
M
M þ A
MA
pðMjA; σ Þ ¼ 2 exp −
:
ð2Þ
I0
2σ 2
σ
σ2
Rician distribution is essentially a special case of the nc-χ
distribution. In this paper, let us define SNR as A/σ. Once the
noiseless signal A tends to zero in image background (i.e., SNR → 0),
the Rician distribution (2) can simplify to a Rayleigh distribution. In
contrast, the noise distribution in high-SNR regions can be
approximated by a Gaussian distribution [51]. Fig. 1 shows the
characteristics of Rician noise in MRI. In low-SNR regions, the noise
distribution cannot be satisfactorily modeled by a Gaussian or
Rayleigh distribution. The signal-dependent Rician noise is neither
additive nor multiplicative in nature. Therefore, to achieve good
restoration results, the MRI denoising methods should take the
characteristics of Rician noise into account.
705
where θ is the estimated SNR of an image and the correction factor ξ(θ)
is defined by
!(
!
!)2
π
θ2
θ2
θ2
2
2
2
ξðθÞ ¼ 2 þ θ − exp −
: ð5Þ
2 þ θ I0
þ θ I1
8
2
4
4
Note that the θ itself can be estimated by satisfying the following
transcendental equation
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
θ ¼ F ðθÞ; F ðθÞ ¼ ζ ðθÞ 1 þ hMi2 =σ 2 −2;
ð6Þ
where σ is the first estimation result using the MAD-based method (3),
〈M〉 is the averaged signal for a given voxel and F(θ) is a function of the
SNR = θ and the given relationship 〈M〉/σ. Fig. 2 (left) shows such
a graphical representation of the function F(θ) with different values of
〈M〉/σ = 1.5 (red), 1.913 (green), 2.0 (blue) and 3.0 (black). The solution
of Eq. (6) is located at the crossing between F(θ) and linear function
θ. In particular, the solution is an exact zero if 〈M〉/σ = 1.913.
If 〈M〉/σ b 1.913, Eq. (6) leads to an imaginary solution. 〈M〉/σ = 1.5
visually illustrates the imaginary solution via the absolute values of F
(θ). Meanwhile, the correction factor ξ(θ) in Fig. 2 (right) is a
monotonically increasing function of SNR. If SNR approaches or
exceeds θ = 5, the correction factor ξ(θ) tends to unity and the
solution of Eq. (6) results in a Gaussian distribution approximation.
To achieve fast convergence at low SNR, we solve the transcendental
Eq. (6) using Newton's method [55], i.e., θt + 1 = θt − G(θt)/G'(θt) with
G(θ) = F(θ) − θ. Once the specified stopping criterion related to θ
update is achieved, we can obtain the correction factor ξ(θ) and spatially
varying noise map σSV accordingly. We refer the interested reader to
Refs. [34,55] for more details.
2.2. Estimation of spatially varying noise map
3. Proposed method
Before the noise reduction step, we first need to estimate the
noise parameters. Recently, many automatic methods have been
proposed for Rician noise estimation. The estimation process is
usually implemented using the properties of Rayleigh distribution in
image background or Gaussian distribution in high-SNR regions. We
refer the interested reader to see Refs. [2,3] and references therein
for more details. But most of the existing noise estimation methods
are based on an assumption of uniform distribution of noise standard
deviation (NSD) across an image. We consider in this paper a twostep wavelet-domain estimation method by assuming the spatially
varying noise map (i.e., spatially varying distributed NSD) [52,53]. In
the first step, the rough NSD is estimated using the Median Absolute
Deviation (MAD) estimator in wavelet domain [54]. To shorten the
implementation time, the rough estimate is a uniform result which is
different from the voxel-by-voxel estimation results in Ref. [34]. In
the second step, we obtain the corrected voxel-by-voxel results
using Koay and Basser's [55] correction scheme for different SNR
values. Under Gaussian assumption, the rough NSD in the first step is
estimated using the following MAD-based method
medianðjsjÞ
σ¼
;
0:6745
ð3Þ
with s being the wavelet coefficients of the highest sub-band HHH,
which is essentially composed of coefficients related to image noise
[54]. In high-SNR regions, the difference between Rician and Gaussian
NSD is negligible. However, the initial noise estimation (3) should be
corrected in low-SNR regions. Manjón et al. [35] and Coupé et al. [54]
have used Koay and Basser's correction scheme [55] to achieve a
uniformly unbiased estimation of σ for all the SNR values. Due to the
different SNR values, the corrected σ should be spatially varying across
the image. In order to obtain the spatially varying noise map σSV, we
correct the MAD estimate σ for each voxel individually as follows
1
σ SV ¼ pffiffiffiffiffiffiffiffiffi σ;
ξðθÞ
ð4Þ
3.1. TV-based MRI denoising
Let Ω be a bounded open subset of ℝd(d = 2, 3) defining the image
domain, f : Ω → ℝ be a degraded image and u : Ω → ℝ be the noiseless
image. Following a Bayesian modeling approach, the TV-based MRI
denoising is equivalent to the following minimization problem [39,40]
min fΦTV ðuÞ ¼ J ðuÞ þ λHðu; f Þg;
u∈BVðΩÞ
ð7Þ
where λ N 0 is a regularization parameter, BV(Ω) is the space of
functions with bounded variation on Ω equipped with the BV
seminorm |u|BV which is formally given by |u|BV = ∫ Ω|∇u| dx,
also referred to as the TV regularization term J(u) = ∫ Ω|∇u| dx
with x ∈ Ω [56]. According to the Rician distribution (2), the datafidelity term related to MRI denoising is given by
!
Z
f 2 þ u2
fu
H ðu; f Þ ¼ − logpð f ju; σ Þ∝
−
logI
dx:
ð8Þ
0
σ2
2σ 2
Ω
Therefore, the TV-based MRI denoising model [39,40] can be defined
as follows
(
)
!
Z
Z
2
2
f þu
fu
−
logI
dx
:
min ΦTV ðuÞ ¼
j∇uj dx þ λ
0
u∈BVðΩÞ
σ2
2σ 2
Ω
Ω
ð9Þ
However, the constant parameter λ over the entire image would limit
the denoising performance. To overcome this limitation, the regularization
parameter should be selected adaptively according to local image features.
In particular, the parameter should be large in texture regions and small in
homogeneous regions. From a statistical point of view, the TV-based
denoising model (9) is essentially based on an assumption of Laplacian
prior on MR image gradients, i.e., J(u) = − log p(u) ∝ ∫ Ω|∇u| dx.
Recently the hyper-Laplacian sparse prior on natural image gradients
has been successfully applied in many image processing applications
[57–59]. We believe there is a significant potential to use the hyperLaplacian prior to reduce Rician noise in MRI.
706
R.W. Liu et al. / Magnetic Resonance Imaging 32 (2014) 702–720
Fig. 3. Statistical analysis of MR brain images. From a statistical point of view, a hyper-Laplacian (blue) with exponent γ ∈ (0,1) is a better model of image gradients than a Gaussian (red) or
a Laplacian (green) for T1w, T2w and PDw images. In particular, the hyper-Laplacian can more closely fit the empirical distribution (black) and guide to feature-preserving MRI denoising.
3.2. Hyper-Laplacian MR image prior
Statistical image analysis has proven to be an enormously
powerful tool in image processing. Recent research shows that the
heavy-tailed marginal distribution of gradients in natural scenes
could be well modeled by the hyper-Laplacian prior [59].
As illustrated in Fig. 3, the heavy-tailed distribution is also found
in MR images. The corresponding hyper-Laplacian image prior can be
modeled as
γ
pðuÞ∝ exp −κ j∇uj ;
ð10Þ
where κ N 0 is a constant scale parameter, and γ ∈ (0, 1) denotes an
exponent parameter required to estimate. According to various MR
imaging features, we respectively collected 50 T1-weighted (T1w),
T2-weighted (T2w) and proton density-weighted (PDw) MR images
from the BrainWeb database 1 [60]. For the sake of simplicity, the
intensity values have been rescaled to [0, 255]. Table 1 summarizes
the estimated exponents γ related to MR image gradients in both x
and y directions for different types of MR images. The difference
1
http://brainweb.bic.mni.mcgill.ca/brainweb/.
R.W. Liu et al. / Magnetic Resonance Imaging 32 (2014) 702–720
Table 1
The estimated exponents γ related to hyper-Laplacian prior for different types of MR
images (T1w, T2w and PDw).
x Direction
Exponent γ
y Direction
T1w
T2w
PDw
T1w
T2w
PDw
0.81
0.61
0.61
0.75
0.63
0.60
between estimation results of T2w and PDw is basically negligible. In
contrast, the T1w leads to larger estimates in both x and y directions.
In practice, the small exponent yields visually sharp result.
However, the image noise may easily be magnified and seriously
degrades the image quality. In contrast, the large value can effectively
remove the noise, but may lead to excessive smoothing. Thus the
optimal exponent should be selected to keep a good balance between
detail preservation and noise reduction. This selection process will be
implemented in Section 4.3.2.
selected adaptively according to local image features. In particular,
the large parameters in texture regions (with small-scale features)
result in fine detail preservation; whereas small ones in homogenous
regions (with large-scale features) perform well in noise reduction.
In Refs. [44–48], the proposed models have taken full advantage of
the local image features and achieved remarkable denoising results.
To the best of our knowledge, no research has been conducted on TVbased MRI denoising models with spatially adaptive regularization
parameters thus far. In this section, we will investigate how to
calculate the adaptive parameters for our proposed GTV-based
denoising model (12).
For the sake of simplicity, the parameter selection method is
developed based on 2D MR image. In practice this method can be
ϖ
naturally generalized to 3D MRI volume. Let Ω(x,y)
⊂ Ω denote the
set of pixel-coordinates in a ϖ-by-ϖ region centered at (x, y) ∈ Ω
ϖ
Ωðx;yÞ ¼
3.3. GTV-based MRI denoising
We adopt the hyper-Laplacian model (10) for representation of
MR image sparse prior. In the MAP framework, the generalized total
variation (GTV) regularization term in our paper is defined by
Z
γ
J ðuÞ ¼ − logpðuÞ∝
j∇uj dx:
ð11Þ
ϖ−1
ϖ−1
≤e
x; e
y≤
;
ðx þ e
x; y þ e
yÞ : −
2
2
where ϖ is a positive odd integer. In our paper, the adaptive
parameter is calculated by considering a local data-fidelity term. To
achieve the term at (x, y) ∈ Ω, we introduce a ϖ × ϖ Gaussian
function K and generate the corresponding term as follows
(
)
2
2
f ðw; zÞuðw; zÞ
f ðw; zÞ þ u ðw; zÞ
K ðw−x; z−yÞ − logI0
dwdz
þ
σ2
2σ 2
Ωϖ
ðx;yÞ
(
)
Z
2
2
f ðw; zÞuðw; zÞ
f ðw; zÞ þ u ðw; zÞ
þ
K ðw−x; z−yÞ − logI 0
¼
dwdz;
σ2
2σ 2
Ω
Z
Ω
By taking the data-fidelity term H(u, f) (8), we can get the GTVbased MRI denoising model min fΦGTV ðuÞ ¼ J ðuÞ þ λHðu; f Þg, i.e.,
u
(
Z
min ΦGTV ðuÞ ¼
u
Ω
j∇uj
γ
Z
dx þ λ
Ω
where the regularization parameter λ N 0 achieves a balance
between the data-fidelity and regularization terms. We assume the
image as a function of space and time, the artificial time-marching
method is then used to seek the solution of (12)
∂u
¼
∂t
F local ðuÞðx; yÞ ¼
)
!
2
2
f þu
fu
dx
;
−
logI
0
2σ 2
σ2
ð12Þ
1
0
I1 fu=σ 2
γðγ−1Þ
γ
λ @
f A;
uNN þ
uTT − 2 u− σ
I0 fu=σ 2
j∇u þ εj2−γ
j∇u þ ε j2−γ
ð13Þ
where |∇u + ε| is a regularized version of |∇u| to reduce
degeneracies in homogeneous regions where |∇u| ≈ 0. We denote
by uNN and uTT the second derivatives of u in the normal and tangent
directions, respectively. The numerator γ(γ − 1) b 0 in (13) means
that the GTV-based denoising model has the properties of backward
diffusion in local normal directions and forward diffusion in local
tangent directions. Thus the model is capable of preserving
discontinuous image features such as edges, lines and other fine
details. In contrast, the TV-based model [39,40] with γ(γ − 1) = 0
can only restore images in the local tangent directions. It is clear that
the TV-based model easily results in an image with staircase effects
in homogeneous regions. As discussed beforehand, the adaptive
parameter λ plays an important role of improving denoising
performance. We will introduce a local variance estimator-based
method to adaptively calculate the parameter related to local image
features and spatially varying noise map.
707
ð14Þ
where the symmetric Gaussian kernel K meets the conditions of K
(w − x, z − y) = K(x − w, y − z) and ∫ K(x, y) dxdy = 1 .
Owing to the introduced local data-fidelity term, the new denoising
version proposed in this paper is defined as locally generalized total
variation (LGTV) model
min
u
Z
ΦLGTV ðu; λÞ ¼
Ω
j∇uj
γ
Z
dxdy þ
Ω
λ F local ðuÞ dxdy :
ð15Þ
If the block size ϖ → ∞, LGTV will simplify to the GTV model (12)
with constant regularization parameter. Inspired by the work in Refs.
[44], the parameter λ is defined as a spatially varying function of
pixel coordinate (x, y). By combining the local data-fidelity term
(14) with parameter λ (see Appendix A for more details), the LGTV
model (15) can be rewritten as follows
(
Z
min ΦLGTV ðu; λÞ ¼
u
γ
Ω
j∇uj
Z
dxdy þ
Ω
ðK⊗λÞ
)
!
2
2
f þu
fu
− logI 0
dxdy :
2
2
2σ
σ
ð16Þ
where ⊗ denotes a convolution kernel. Let ψ(∘) = I1(∘)/I0(∘), the
Euler–Lagrange equation associated with (16) is given by
−∇ γ∇u
K⊗λ
fu
þ
u−ψ
f ¼ 0;
j∇u þ εj2−γ
σ2
σ2
ð17Þ
3.4. Automated selection of spatially adaptive regularization parameters
with the Neumann boundary condition. In this paper, the spatially
adaptive regularization parameters are defined as follows
MR image contains multiple objects of different scales, thus the
constant parameter λ over the entire image is detrimental to
achieving satisfactory visual quality. The parameter should be
λðx; yÞ ¼
Q ðx; yÞ
;
Sðx; yÞ
ðx; yÞ∈Ω;
ð18Þ
708
R.W. Liu et al. / Magnetic Resonance Imaging 32 (2014) 702–720
where
4. Experiments and results
8
∇u
2
>
< Q ðx; yÞ ¼ γσ 2 ∇
f
u−ψ
fu=σ
j∇u þ εj2−γ
:
2
>
: Sðx; yÞ ¼ K⊗ u−ψ fu=σ 2 f
To evaluate and compare the competing denoising methods, a
set of experiments were performed on both synthetic and real MR
data sets.
ð19Þ
4.1. Comparison with other denoising techniques
Appendix B gives more details on automated selection of spatially
adaptive regularization parameters. For the sake of better reading,
we use ψ(u) instead of ψ( fu/σ 2) throughout the rest of this paper. If
ψ(u) → 1 in high-SNR regions, ψ(u)f simplifies back to the observed
degraded image, and the noise can be approximated by a Gaussian
distribution; whereas if ψ(u) → 0 in image background, the noise
can simplify to a Rayleigh distribution, and the magnitude of residual
image u − ψ(u)f tends to become quite negligible. Inspired by the
work in Refs. [44–48], the S(x, y) can be approximated by
Sðx; yÞ≈
σ 4SV ðx; yÞ
;
LV ðx; yÞ
ð20Þ
where σSV denotes the spatially varying noise standard deviation. If a
high-quality denoised version u is obtained, the corresponding local
variance of residual image ϕ ¼ u−ψðuÞf is given by
2
LV ðx; yÞ ¼ K⊗ ϕ−ϕ ;
where ϕ is the mean value of ϕ. For the sake of better understanding,
the noiseless image u is decomposed into two components: u =
uC + uT, where uC and uT denote the cartoon and texture regions,
respectively. If the image is cartoon-like (i.e., uT is close to zero),
many denoising methods could generate satisfactory results such
that u≈uC ¼ u and ϕ ≈ u − ψ(u)f. However, the magnitude MR
images commonly contain fine structural features, which include a
wealth of useful medical information. Until now, no technique can
completely extract the noiseless image from its noisy version. Some
fine texture details may be filtered out and included in the residual
image ϕ. Thus the local variance LV(x, y) should be approximated by
2
the sum of texture local variance LVtexture(x, y) and σSV
(x, y). In this
paper, Eq. (20) can be rewritten in the following form
4
Sðx; yÞ≈
4
σ SV ðx; yÞ
σ SV ðx; yÞ
2
≤σ SV ðx; yÞ:
2 ¼
LV texture ðx; yÞ þ σ 2SV ðx; yÞ
K⊗ ϕ−ϕ
ð21Þ
The local regularization parameters in texture regions are larger
than those in homogeneous regions. As discussed in Ref. [44], large
regularization parameters λ(x, y) in texture regions result in fine
detail preservation; whereas small ones in homogeneous regions
perform well in noise reduction. Thus the meaningful features could
be effectively preserved in MR images. In numerical implementation,
solution of the Euler–Lagrange Eq. (17) associated with LGTV is
achieved by iteratively evolving the following negative gradient flow
∂u
¼
∂t
γðγ−1Þ
γ
K⊗λ
− 2 ðu−ψðuÞf Þ in
u þ
u
2−γ NN
2−γ TT
j∇u þ εj
j∇u þ εj
σ
ð0; Τ Þ Ω;
and updating the spatially adaptive regularization parameters
defined by Eq. (18) until convergence. Theoretically the proposed
non-convex LGTV model is not well posed and has high computational complexity, but we can still achieve satisfactory denoising
performance in numerical experiments. In particular, this nonconvex model can enhance the structural differences between
various components in favor of sparse gradients.
Our proposed LGTV model was compared to six recently
developed methods as follows.
• RLMMSE: Recursive version of Linear Minimum Mean Square Error
Estimator [1]. This method presented by Aja-Fernández et al. [1]
has shown a good performance in both noise reduction and feature
preservation. In all cases, the search volume of 11 × 11 × 11 voxels, local neighborhood of 3 × 3 × 3 voxels and 5 iterations have
been used.
• RSNLMMSE: Recursive version of SNR-based Nonlocal LMMSE
[27]. RSNLMMSE was proposed based on image data redundancy and local SNR estimation. In the experiments, the same
parameters mentioned in RLMMSE were also adopted to
achieve a good balance between noise reduction and computational complexity.
• ODCT: Oracle-based Discrete Cosine Transform [18]. This method took
full advantage of the sparseness property of MR image, and was
developed based on a 3D moving-window DCT hard thresholding.
• UKR: Unbiased Kernel Regression filter [7]. The nonparametric
estimation method UKR was dependent on a zeroth order 3D
kernel regression, which computed a weighted average of pixels
over a regression window.
• WSM: Wavelet Subbands Mixing [14]. WSM was in essence a hybrid
denoising approach which combined the optimized blockwise
NLM filter [13] with 3D discrete wavelet transform (DWT).
• AONLM: Adaptive Optimized Nonlocal Means [35]. This method took
into consideration both the Rician nature of MR data and spatially
varying noise properties. It was the first approach developed for
denoising of MR images with spatially varying noise levels. The
search volume of 7 × 7 × 7 voxels and local neighborhood of
3 × 3 × 3 voxels were adopted in the experiments.
The noise reduction performance of these competing methods
was evaluated under both the objective criterion and subjective
criterion of visual quality.
4.2. Synthetic and in vivo real data sets
To compare the denoised results with a benchmark database, the
BrainWeb [60] was adopted to conduct synthetic experiments. We
considered the 3D PDw, T1w and T2w volumes of
181 × 217 × 181 voxels with zero noise and 1 × 1 × 1 mm 3 voxel
resolution. We employed the 12bit precision data where the original
values are in the range [0, 4095]. In the experiments, the entire range
[0, 4095] was scaled to [0, 255], and different levels of Rician noise
were added to the noiseless MRI data.
In vivo experiments were also performed to evaluate the
denoising performance of LGTV. Our DT-MRI data were acquired
using a clinical 3 T MRI scanner with an 8-channel SENSE head coil
(Achieva, Philips Medical Systems, Best, the Netherlands) at the
Prince of Wales Hospital in Hong Kong. DT-MRI was obtained using
single-shot echo planar imaging sequence with the following
parameters: 32 diffusion-weighted volumes (b = 1000 s/mm 2),
one diffusion-weighted volume (b = 0 s/mm 2), TR = 8667 ms,
TE = 60 ms, FOV = 224 × 224 mm 2 , NEX = 1, matrix =
112 × 109, slice = 70, slice thickness = 2 mm and gap = 2 mm.
After reconstruction, images were zero-padded and interpolated to
224 × 224 with spatial resolution at 1 × 1 × 1 mm 3.
R.W. Liu et al. / Magnetic Resonance Imaging 32 (2014) 702–720
709
Fig. 4. Effects of different exponents γ on denoising of different MR images corrupted with different levels of Rician noise (10% to 25%). The metrics PSNR and SSIM are used to
evaluate the denoising performance. In particular, a higher PSNR or SSIM value correlates to a higher-quality image.
4.3. Evaluation on synthetic data
Structural Similarity index (SSIM), which is more consistent with
human visual perception [61].
4.3.1. Image quality measures
In order to objectively evaluate the denoising performance, three
quality measures were used simultaneously. The first measure was
the Peak Signal-to-Noise Ratio (PSNR), i.e.,
PSNRðu; uÞ ¼ 10 log10 X
2552 M N
ðx;yÞ∈Ω
2
juðx; yÞ−uðx; yÞj
ðdBÞ;
ð22Þ
with the noiseless M × N monochrome image u, degraded image f
and its filtered version u. The second performance measure was the
ð2μ μ þ c1 Þð2σ uu þ c2 Þ
;
SSIMðu; uÞ ¼ 2 u 2 u
μ u þ μ u þ c1 σ 2u þ σ 2u þ c2
ð23Þ
where μu and μ u are the local mean values of images u and u, σu and
σ u represent the respective standard deviations, σ uu is the
covariance value, and c1, c2 are two constants to avoid instability.
Based on the assumption that a great amount of structural
information of an image was coded in its local variance distribution, Aja-Fernández et al. [62] originally developed a new image
710
R.W. Liu et al. / Magnetic Resonance Imaging 32 (2014) 702–720
Fig. 5. From top-left to bottom-right: the spatially varying Rician noise maps are respectively estimated from the PDw, T1w and T2w MR images corrupted with different levels of
Rician noise (10% to 25%).
quality assessment method, i.e., Quality Index based on Local
Variance (QILV).
2μ V ðuÞ μ V ðuÞ þ c1 2σ V ðuÞV ðuÞ þ c2
;
QILVðu; uÞ ¼ ð24Þ
μ 2V ðuÞ þ μ 2V ðuÞ þ c1 σ 2V ðuÞ þ σ 2V ðuÞ þ c2
where σ V ðuÞV ðuÞ is the covariance value, and μV(u) and σV(u) denote
the expected value and standard deviation of the local variances V
(u) and V ðuÞ, respectively. The QILV index takes values in [0, 1] and
increases as the image quality becomes better. In particular, the
QILV is more sensitive to the amount of blurring of image edges
caused by different denoising methods [62,63].
4.3.2. Optimal exponent selection
To guarantee high-quality denoising results, the exponent γ
associated with hyper-Laplacian prior should be selected properly. In
particular, the smaller exponent generates visually sharper image;
whereas, the larger one can remove noise and produce more smooth
version. If we directly use the exponent (γ b 1) estimated in
Section 3.2, the LGTV-based denoising result may suffer from
excessive spatial-sharpening and degrade image quality. In order
to select the optimal exponent, the effects of different exponents on
denoising results were investigated for different types of MR images.
In experiments, the PDw, T1w and T2w MR images were
corrupted with different levels of Rician noise (10% to 25%). Both
PSNR and SSIM were used to investigate the effects of different
exponents on denoising of different MR images. The experimental
results at different noise levels could be seen in Fig. 4. As can be
observed, the exponents should be different for different types of MR
images. Take PDw image as an example, the exponent γ = 0.9
outperformed other exponents under consideration in most of the
cases. To maintain a balance between PSNR and SSIM metrics, the
exponents γ = 0.8 and 0.7 were selected for the T1w and T2w MR
images, respectively.
Before doing noise reduction, the MAD-based method (3) was
first used to estimate the rough Rician noise. The final spatially
varying corrected Rician noise fields were achieved using the
correction scheme (4). As shown in Fig. 5, the Rician NSD tends to
Gaussian NSD in high-SNR regions, which is lower than Rician NSD in
low-SNR regions. In particular, the highest NSD occurs in homogeneous background and generates the smallest regularization parameters. This characteristic was confirmed by the estimated spatially
adaptive regularization parameters shown in Fig. 6. In particular,
small regularization parameters in homogeneous regions can
perform well in noise reduction; whereas large ones in texture
regions can result in fine detail preservation.
4.3.3. Quantitative comparison of denoising performance
This section is devoted to compare our proposed LGTV model
with some recently proposed related MRI denoising methods. The
experimental results were obtained using the simulated Brainweb
data sets. To evaluate the stability of our denoising method, the
simulated MRI volumes were corrupted with different levels of
Rician noise ranged from 5% to 25% of the maximum intensity with
5% in step. To evaluate the denoising performance of different
methods, three metrics (i.e., PSNR, SSIM and QILV) were used
simultaneously. Tables 2–4 depict the quantitative results with
R.W. Liu et al. / Magnetic Resonance Imaging 32 (2014) 702–720
711
Fig. 6. From top-left to bottom-right: the spatially adaptive regularization parameters are respectively estimated from the PDw, T1w and T2w MR images corrupted with different
levels of Rician noise (10% to 25%). In experiments, the optimal exponents are chosen as γ = 0.9, 0.8 and 0.7 for the PDw, T1w and T2w scenarios, respectively.
different image types and noise levels. LGTV outperformed other
methods under consideration in most of the cases. In contrast, UKR
slightly yielded better results when PDw volume was corrupted with
Rician noise levels of 5% and 10%. For noise levels of 15% and 20%, ODCT
generated the most satisfactory performance on the T2w volume in
terms of the PSNR metric. In terms of the QILV index, however, LGTV
showed improvements over all other methods in all cases.
4.3.4. Qualitative visual quality assessment
Visual quality comparison is an important criterion to judge
the performance of a denoising method. In practice, the denoised
version should contain noticeable geometrical structures, few or
(ideally) no visible artifacts. Some of the visual results are illustrated
in Figs. 7–11, which show the comparison between different denoising methods. The corresponding estimated noise maps and final
regularization parameters can be observed in Figs. 5 and 6,
respectively. As shown in Fig. 7, a simulated PDw MR volume was
corrupted with a Rician noise level of 10%. All methods mentioned in
this paper could effectively reduce noise in this case. From the
magnified 1D profiles shown in pink insets, it can be observed that
the intensity values of AONLM and LGTV are more structurally
similar to the original image.
We further compared the visual quality of the denoising results. A
simulated T1w MR volume was corrupted with a Rician noise level of
15%. The denoising results and their associated magnified views are
displayed in Figs. 8 and 9, respectively. The undesirable biases
Table 2
Quantitative evaluation of various denoising methods on a simulated PDw MR data set.
Noise level
Noisy data
RLMMSE
RSNLMMSE
ODCT
UKR
WSM
AONLM
LGTV
5%
10%
15%
PSNR
SSIM
QILV
PSNR
SSIM
QILV
PSNR
SSIM
QILV
25.822
32.229
33.749
36.066
36.472
35.287
32.194
36.047
0.6366
0.9062
0.9430
0.9597
0.9675
0.9143
0.8407
0.9655
0.9907
0.9981
0.9993
0.9996
0.9996
0.9990
0.9980
0.9997
18.276
22.613
27.613
28.602
29.013
28.197
25.529
28.974
0.4353
0.7835
0.7595
0.8153
0.8911
0.7701
0.7408
0.8983
0.9073
0.9832
0.9927
0.9957
0.9971
0.9942
0.9903
0.9975
15.058
25.165
26.901
28.027
27.532
27.032
22.782
28.242
0.3373
0.7326
0.7662
0.8008
0.8399
0.7282
0.6919
0.8736
0.7776
0.9813
0.9888
0.9956
0.9954
0.9943
0.9741
0.9965
The results are shown for different levels of Rician noise. The best value of each column is highlighted with bold.
712
R.W. Liu et al. / Magnetic Resonance Imaging 32 (2014) 702–720
Table 3
Quantitative evaluation of various denoising methods on a simulated T1w MR data set.
Noise level
Noisy data
RLMMSE
RSNLMMSE
ODCT
UKR
WSM
AONLM
LGTV
10%
15%
20%
PSNR
SSIM
QILV
PSNR
SSIM
QILV
PSNR
SSIM
QILV
17.013
22.248
25.378
26.246
26.414
25.505
22.934
26.925
0.4926
0.8091
0.7602
0.7894
0.8942
0.7456
0.7195
0.9084
0.8718
0.9816
0.9826
0.9885
0.9883
0.9865
0.9775
0.9915
14.171
21.742
24.527
25.329
25.019
24.245
21.229
26.311
0.3716
0.7131
0.7651
0.7669
0.8482
0.7048
0.6765
0.8883
0.6830
0.9710
0.9777
0.9830
0.9814
0.9842
0.9658
0.9911
12.316
17.242
22.506
24.290
22.422
21.886
19.903
25.789
0.2796
0.6667
0.6732
0.7203
0.6952
0.6628
0.6597
0.8677
0.4569
0.9147
0.9077
0.9683
0.9759
0.9648
0.9548
0.9902
The best value of each column is highlighted with bold.
estimated by RSNLMMSE on borders and homogeneous regions
resulted in visual quality degradation. As shown in Fig. 9, UKR, WSM
and AONLM were unable to completely reduce the noise in regions of
homogeneous. ODCT suffered from slightly oversmoothing of texture
regions. In contrast, LGTV reduced the image noise almost completely
in homogeneous regions and preserved the edges. It means that our
proposed model could keep a good balance between noise reduction
and detail preservation. The advantage of LGTV was further confirmed
by the surface plots of denoised images shown in Fig. 10. As shown by
the color arrows, LGTV performed very well on the T1w volume, and
produced the most similar surface plot to that of the original version.
Definition 1. Method error
Let u be the original image and f be the observed degraded
version. The method error Λ can be defined as the absolute difference
between the original and denoised images
Λ ¼ ju−ρð f Þj:
where ρ denotes a denoising operator. As discussed in Ref. [32], the
method error should contain as little structural information as possible.
If f is just the original image u, the method error Λ becomes equivalent
to method noise [29]. In this case, it also can make sense to evaluate any
denoising method without the traditional “add noise and then reduce
it” trick.
In Fig. 11, we compared our LGTV model with RSNLMMSE, ODCT,
UKR and WSM on a simulated T2w MR volume corrupted with a
Rician noise level of 20%. As can be observed, the methods
(RSNLMMSE, ODCT and LGTV) showed better performance in both
noise reduction and detail preservation. The corresponding method
errors appear more random and smaller compared with other
methods. Both UKR and WSM failed to preserve fine details, probably
because of the excessive smoothing over the texture regions. The
geometrical structures were noticeable in the method errors. It is
well known that both RSNLMMSE and ODCT take full advantage of
high degree of pattern redundancy in MR images, thus they can
generate high-quality results. In contrast, the superior performance
of LGTV benefits from the global hyper-Laplacian image prior and
spatially adaptive regularization parameters.
4.4. Evaluation on in vivo real data
4.4.1. Diffusion tensor metrics
In order to evaluate the consistency of the LGTV model on real
clinical data, an in vivo brain DT-MRI data set was used. Without loss
of generality, the exponent γ = 0.8 was selected in this experiment.
Several diffusion tensor metrics, such as Fractional Anisotropy (FA),
color-coded FA (CFA), the largest Eigenvalue (L1) and Apparent
Diffusion Coefficient (ADC), were simultaneously adopted to measure the DT-MRI denoising results. In particular, FA is a scalar value
between 0 and 1 that determines the fraction of the diffusion tensor.
CFA indicates the orientation of the principal eigenvector of the
diffusion tensor. L1 corresponds to the principal eigenvector. ADC
represents the trace of tensors [64]. The differences in FA, L1 or ADC
could reflect the denoising performance of all considered methods.
The rough Rician noise level was first estimated to be around 1.02%
of the maximum intensity using the MAD-based method (3). As shown
in Fig. 12, the correction scheme (4) was then adopted to achieve the
spatially varying noise maps for different diffusion weightings: b =
0 s/mm 2 (top-left) and b = 1000 s/mm2 (top-right). The spatially
adaptive regularization parameters (bottom row in Fig. 12) related to
LGTV-based denoised DT-MR images could be achieved accordingly. In
these two steps, the DT-MR images were scaled to a common range of
values [0, 255] and rescaled to original intensity range for diffusion
tensor metric measurements and fiber tracking.
Fig. 13 shows the scale maps of FA, CFA, L1 and ADC before and
after denoising by UKR, WSM, AONLM, RSNLMMSE and LGTV,
respectively. The FA and CFA maps computed from the original DTMRI had noticeable granular aspects. Both the original and WSM
suffered from the black spots in FA and CFA maps (shown by
the white circles), which correspond to tensors having negative
eigenvalues. As shown by the white arrows in local CFA, the
Table 4
Quantitative evaluation of various denoising methods on a simulated T2w MR data set.
Noise level
Noisy data
RLMMSE
RSNLMMSE
ODCT
UKR
WSM
AONLM
LGTV
15%
20%
25%
PSNR
SSIM
QILV
PSNR
SSIM
QILV
PSNR
SSIM
QILV
17.199
24.020
25.229
27.605
27.278
25.309
22.426
26.920
0.4906
0.8070
0.7567
0.8321
0.8139
0.7541
0.7538
0.8730
0.8665
0.9554
0.9892
0.9882
0.9877
0.9831
0.9763
0.9897
15.320
22.414
23.575
25.441
24.402
23.832
20.695
25.303
0.4177
0.7381
0.7242
0.7666
0.7509
0.7042
0.7095
0.8527
0.7613
0.9199
0.9647
0.9685
0.9773
0.9657
0.9614
0.9827
14.063
20.857
21.750
23.512
22.801
21.250
20.366
24.001
0.3544
0.6526
0.6591
0.7143
0.7115
0.6412
0.6643
0.8225
0.6354
0.8530
0.9146
0.9250
0.9709
0.9242
0.9419
0.9758
The best value of each column is highlighted with bold.
R.W. Liu et al. / Magnetic Resonance Imaging 32 (2014) 702–720
713
Fig. 7. 2D axial slices (top) and their associated 1D profiles (bottom). Top: Results of different denoising methods on an axial slice of the simulated PDw MR volume corrupted
with a Rician noise level of 10%. Bottom: The intensity values of 1D profiles horizontally through the original image, noisy image and denoised versions yielded by RSNLMMSE,
ODCT, UKR, WSM, AONLM and LGTV, respectively.
RSNLMMSE method showed some artifacts on borders and homogeneous areas due to inaccurate Rician noise bias correction. AONLM
could overcome this limitation but seemed to overcome some DTMRI details in L1 map. Compared to WSM, AONLM and RSNLMMSE,
both UKR and LGTV not only preserve the L1 of the diffusion tensor,
but also result in more regular FA and CFA maps without obvious
visual artifacts. It is also worth to note that the visual quality of ADC
map was significantly enhanced by all denoising methods.
Fig. 8. Example denoising results for an axial slice of the simulated T1w MR volume corrupted with a Rician noise level of 15%.
714
R.W. Liu et al. / Magnetic Resonance Imaging 32 (2014) 702–720
Fig. 9. Magnified views of the white square regions shown in Fig. 8.
Fig. 10. Surface plots of the magnified views shown in Fig. 9. The original image, noisy image and denoised versions generated by ODCT, WSM, AONLM and LGTV are displayed.
The qualitative analysis in Fig. 13 was confirmed by the quantitative
results shown in Table 5. It is proved by MonteCarlo simulation (MCS)
[65,66] that the eigenvalues are commonly biased by Rician noise.
Accordingly, the FA, L1 and ADC related to eigenvalues also suffer from
biases in practice. From the Table 5, we observe that denoising drastically
decreases the mean values of FA, L1 and ADC. LGTV outperforms all other
methods in terms of L1 and ADC. In contrast, UKR only performs slightly
better than LGTV in terms of FA. RLMMSE and ODCT lead to the least
improvement for all three diffusion tensor metrics.
4.4.2. Fiber tracking
It is well known that the fiber track propagates along the direction
of the principal eigenvector of the diffusion tensor from one position to
the next throughout the whole brain. Consequently, the tractography
Fig. 11. Qualitative comparison of different methods under consideration on a simulated T2w MR volume corrupted with a Rician noise level of 20%. The first row shows three
continuous slices of the T2w volume. For the sake of comparison, the other rows only illustrate the local magnification views shown in the white square regions. As can be
observed, the RSNLMMSE, ODCT and LGTV show better performance in both noise reduction and detail preservation.
R.W. Liu et al. / Magnetic Resonance Imaging 32 (2014) 702–720
715
716
R.W. Liu et al. / Magnetic Resonance Imaging 32 (2014) 702–720
Fig. 12. From top-left to bottom-right: spatially varying Rician noise maps and spatially adaptive regularization parameters estimated from diffusion images at b-values equal to 0 and
1000 s/mm2, respectively. In experiments, the diffusion images were scaled to a common range of values [0, 255] and rescaled to original intensity range after noise estimation and reduction.
experiments were also performed to evaluate the effects of different
denoising methods on the directional information of diffusion of water
molecules. Fig. 14 shows the fiber tracking results—using MedINRIA
software2—of the whole human brain before and after denoising with
the RLMMSE, RSNLMMSE, ODCT, UKR, WSM, AONLM and LGTV,
respectively. The fibers extracted from the original DT-MRI data set are
not well reconstructed and are somewhat irregular. In contrast, the
visual quality of the fibers reconstructed from the denoised tensors has
been significantly enhanced. However, the RLMMSE, RSNLMMSE and
AONLM still suffer from slightly irregular trajectories, as shown by the
white arrows in local magnification views. Fig. 15 shows a local fiber
bundle extracted from a manually selected region of interest (ROI)
area. The arrows point out erroneous and irregular trajectories
resulting from noise in the original DT-MRI. The tracking yielded
from the denoised tensors is smoother and more organized than that
obtained from the original tensors. In particular, the RLMMSE,
RSNLMMSE, UKR and WSM could considerably enhance the performance of fiber tracking, but somewhat fail in reconstruction of some
fine fiber bundles. Seemingly the ODCT, AONLM and LGTV produce the
most satisfactory results. As observed in Fig. 14, the AONLM yielded
irregular trajectories due to oversmoothing of the largest Eigenvalue
(shown in Fig. 13). The local magnification views in yellow boxes
illustrate that the tracking obtained from our proposed LGTV is much
more organized than that performed from the ODCT.
5. Discussion and conclusion
In this paper, we have proposed an LGTV-based MRI Rician
denoising model. The major contributions of this paper are as
2
http://www-sop.inria.fr/asclepios/software/MedINRIA/
follows: (1) The denoising method proposed in this paper was
implemented based on an assumption of spatially varying noise map.
A two-step wavelet-domain estimation method has been introduced
to extract the noise map, which played an important role in
automated selection of spatially adaptive regularization parameters.
(2) The hyper-Laplacian sparse prior on image gradients was
adopted to enhance the feature-preserving denoising property of
LGTV. In particular, LGTV has the compelling properties of backward
diffusion in local normal directions and forward diffusion in local
tangent directions. (3) To further improve the denoising performance, a local variance estimator-based method was developed to
calculate the spatially adaptive regularization parameters related to
local image features and spatially varying noise map. Thus the
proposed model can effectively reduce the noise level while
maintaining the image features.
It is well known that the nc-χ distribution is becoming more and
more common to model noisy magnitude MR signals in multiple-coil
acquisition systems. Our proposed Rician-distributed-based LGTV
denoising model could introduce unwanted bias in this case.
According to the noise characteristics in MRI, Rician distribution is
in essence a special case of the nc-χ distribution. Naturally the basic
idea behind our denoising approach can be generalized to the
reduction of nc-χ distributed noise. Recently the parallel MRI (pMRI)
techniques have been employed to accelerate the acquisition process
by subsampling k-space data. However these techniques can affect
the statistical distribution of magnitude signal and cause the noise to
be non-uniform across an image [2,3]. In future research, it would be
useful to estimate the non-uniform noise and correct it using Koay
and Basser's correction scheme [55]. From a practical point of view,
the proposed LGTV model can be further modified by taking into
account both nc-χ distribution and non-uniform nature of noise.
R.W. Liu et al. / Magnetic Resonance Imaging 32 (2014) 702–720
717
Fig. 13. Diffusion tensor metrics. From top to bottom: the original DT-MRI data yielded maps of fractional anisotropy (FA), color-coded FA (CFA), local magnification of CFA within a view
(Local CFA), the largest eigenvalue (L1) and apparent diffusion coefficient (ADC), which were improved by all the competing denoising methods (UKR, WSM, AONLM, RSNLMMSE and LGTV).
In addition, experimental results show that LGTV suffers from an
intrinsic limitation of high computational cost in practical cases.
Nevertheless, the proposed model is still worthy of consideration since
Table 5
Statistics of FA, L1 and ADC computed from the real brain DT-MRI and from the
restored results obtained by the competing denoising methods.
L1 (mm2)
FA
Original
RLMMSE
RSNLMMSE
ODCT
UKR
WSM
ANOLM
LGTV
ADC (mm2/s)
Mean
Std
Mean
Std
Mean
Std
0.3443
0.3246
0.3222
0.3352
0.3073
0.3174
0.3200
0.3094
0.1674
0.1611
0.1582
0.1660
0.1532
0.1587
0.1616
0.1545
1.4573
1.3932
1.2996
1.4132
1.2408
1.2860
1.4197
1.2352
0.6947
0.7287
0.5421
0.7458
0.4951
0.5423
0.8731
0.4841
3.2558
3.1141
2.9359
3.1660
2.8262
2.9127
3.1732
2.8081
1.7158
1.6246
1.3771
1.8340
1.2692
1.3916
1.9230
1.2404
The best value of each column is highlighted with bold.
it generates denoising results which are quantitatively and qualitatively comparable with some current state-of-the-art methods.
Graphics processing unit (GPU) has rapidly evolved into an acceleration technology for high-performance computing, especially in
medical image processing over the past few decades [67,68]. For
instance, Bui et al. [69] proposed an efficient GPU-based acceleration
technique and applied it on TV-based MRI Rician denoising [39].
Therefore, there is a strong incentive to accelerate LGTV for real-time
applications in the GPU-based parallel computation framework.
Acknowledgments
The authors would like to thank the editor and anonymous reviewers
for their valuable comments and suggestions that have led to a substantial
improvement in the paper. The work described in this paper was
supported by grants from the Research Grants Council of the Hong Kong
Special Administrative Region, China (Project No.: CUHK 475711, 416712
718
R.W. Liu et al. / Magnetic Resonance Imaging 32 (2014) 702–720
Fig. 14. From top-left to bottom-right: whole-brain fiber tracking obtained after estimating the diffusion tensor from original DT-MRI volume, filtered volumes generated by
RLMMSE, RSNLMMSE, ODCT, UKR, WSM, AONLM and LGTV, respectively. The fiber bundles have been colored according to the fractional anisotropy (normalized variance of the
eigenvalues of the diffusion tensor) at each location. The local magnification views are shown in the insets.
and 473012), grants from the National Natural Science Foundation of
China (Project No.: 81201157 and 81101111), a grant from BME-p2-13/
BME-CUHK of the Shun Hing Institute of Advanced Engineering, The
Chinese University of Hong Kong, a grant from the Science, Industry, Trade
and Information Commission of Shenzhen Municipality (Project No.:
JC201005250030A), the 863 Program of China (Project No.:
2012AA02A603), and the Research Fund for the Doctoral Program of Higher
Education of China (PhD supervisor grant) (Project No.: 20134433110012).
Fig. 15. Local-brain fiber tracking was achieved from a manually selected region of interest (ROI) area.
R.W. Liu et al. / Magnetic Resonance Imaging 32 (2014) 702–720
719
Appendix A. Locally generalized total variation model
By combining the local data-fidelity term Flocal(u)(x, y) (14) with spatially adaptive regularization parameters λ(x, y), we get
Z
λðx; yÞF local ðuÞðx; yÞ dxdy
(Z
"
)
#
Z
f 2 ðw; zÞ þ u2 ðw; zÞ
f ðw; zÞuðw; zÞ
¼
λðx; yÞ
K ðw−x; z−yÞ
−
logI
dwdz
dxdy
0
2σ 2
σ2
Ω
Ω
"
#
Z Z
f 2 ðx; yÞ þ u2 ðx; yÞ
f ðx; yÞuðx; yÞ
λðw; zÞK ðx−w; y−zÞ
− logI0
dxdydwdz
¼
2
2
2σ
σ
Ω
Ω
" 2
#
Z Z
2
f ðx; yÞ þ u ðx; yÞ
f ðx; yÞuðx; yÞ
K ðx−w; y−zÞλðw; zÞ dwdz
− logI0
dxdy
¼
2σ 2
σ2
Ω
Ω
"
#
Z
2
2
f ðx; yÞ þ u ðx; yÞ
f ðx; yÞuðx; yÞ
ðK⊗λÞðx; yÞ
− logI0
dxdy;
¼
2
2σ
σ2
Ω
Ω
where ⊗ denotes a convolution kernel. We can obtain a new version of the locally generalized total variation (LGTV) model (16) as follows
(
)
!
Z
Z
f 2 þ u2
fu
γ
j∇uj dxdy þ
ðK⊗λÞ
− logI0
dxdy :
min ΦLGTV ðu; λÞ ¼
u
2σ 2
σ2
Ω
Ω
Appendix B. Spatially adaptive regularization parameters
We herein analytically derive the mathematical equation for automated selection of spatially adaptive regularization parameters.
Let ψ(∘) = I1(∘)/I0(∘), the Gâteaux derivative of ΦLGTV(u, λ) (16) with respect to u can be calculated in the direction v using standard
variational method.
lim
α→0
1
ðΦ
ðu þ αv; λÞ−ΦLGTV ðu; λÞÞ
α LGTV
(
)
Z
Z
2
2
d
ðu þ αvÞ þ f
ðu þ αvÞf
γ
dxdy
j
¼
ðK⊗λÞ
−
logI
j∇u þ α∇vj dxdy þ
0
dα α¼0 Ω
σ2
2σ 2
Ω
Z
Z
γ∇u
1
fu
¼
∇v dxdy þ
ðK⊗λÞ u−ψ 2 f v þ uv dxdy
2−γ
2
j∇u
þ
εj
σ
σ
Ω
Ω
Z
Z γ∇u
K⊗λ
fu
γ∇u
!
−∇ þ
u−ψ
v
f
v
dxdy
þ
n dxdy;
¼
j∇u þ εj2−γ
σ2
σ2
j∇u þ εj2−γ
Ω
∂Ω
!
where n is the unit outward normal to ∂Ω. The Euler–Lagrange equation associated with LGTV (16) is given by
−∇ γ∇u
K⊗λ
fu
þ
u−ψ
f
¼ 0:
σ2
σ2
j∇u þ εj2−γ
ðB1Þ
To calculate the parameter λ, multiplying the Eq. (B1) by (u − ψ(fu/σ 2)f) and then integrating on image domain Ω, we can obtain
Z
γ∇uðx; yÞ
f ðx; yÞuðx; yÞ
∇
uðx; yÞ−ψ
f ðx; yÞ dxdy
2−γ
2
j∇uðx; yÞ þ εj
σ
Ω
2
Z Z
1
f ðx; yÞuðx; yÞ
¼
K
ð
x−w;
y−z
Þλ
ð
w;
z
Þ
dwdz
u
ð
x;
y
Þ−ψ
dxdy
f
ð
x;
y
Þ
σ2 Ω
σ2
Ω (
)
2
Z
Z
1
f ðw; zÞuðw; zÞ
¼
λðx; yÞ
K ðx−w; y−zÞ uðw; zÞ−ψ
f ðw; zÞ dwdz dxdy:
σ2 Ω
σ2
Ωðϖx;yÞ
Therefore,
γ∇uðx; yÞ
f ðx; yÞuðx; yÞ
uðx; yÞ−ψ
f ðx; yÞ
2−γ
2
j∇uðx; yÞ þ εj
(Z σ
)
2
1
f ðw; zÞuðw; zÞ
¼
λ
ð
x;
y
Þ
K
ð
x−w;
y−z
Þ
u
ð
w;
z
Þ−ψ
dwdz
:
f
ð
w;
z
Þ
σ2
σ2
Ωϖ
ðx;yÞ
The spatially adaptive regularization parameter λ(x, y) is achieved using λ(x, y) = Q(x, y)/S(x, y). The Q(x, y) and S(x, y) are respectively
defined as follows
∇
2
Q ðx; yÞ ¼ σ γ ∇
∇u
fu
u − ψ 2 f
;
2−γ
j∇u þ εj
σ
and
2
fu
Sðx; yÞ ¼ K⊗ u−ψ 2 f :
σ
720
R.W. Liu et al. / Magnetic Resonance Imaging 32 (2014) 702–720
References
[1] Aja-Fernández S, Alberola-Lopez C, Westin CF. Noise and signal estimation in
magnitude MRI and Rician distributed images: a LMMSE approach. IEEE Trans
Image Process 2008;17(8):1383–98.
[2] Aja-Fernández S, Tristán-Vega A, Alberola-Lopez C. Noise estimation in singleand multiple-coil magnetic resonance data based on statistical models. Magn
Reson Imaging 2009;27(10):1397–409.
[3] Aja-Fernández S, Brion V, Tristán-Vega A. Effective noise estimation and
filtering from correlated multiple-coil MR data. Magn Reson Imaging
2013;31(2):272–85.
[4] Henkelman RM. Measurement of signal intensities in the presence of noise in MR
images. Med Phys 1985;12(2):232–3.
[5] Awate SP, Whitaker RT. Nonparametric neighborhood statistics for MRI
denoising. Proc IPMI; 2005. p. 677–88.
[6] Awate SP, Whitaker RT. Feature-preserving MRI denoising: a nonparametric
empirical Beyes approach. IEEE Trans Med Imaging 2007;26(9):1242–55.
[7] López-Rubio E, Florentín-Núñez MN. Kernel regression based feature extraction
for 3D MR image denoising. Med Image Anal 2011;15(4):498–513.
[8] Mukherjee PS, Qiu P. 3-D image denoising by local smoothing and nonparametric regression. Technometrics 2011;53(2):196–208.
[9] Rabbani H, Nezafat R, Gazor S. Wavelet-domain medical image denoising using
bivariate Laplacian mixture model. IEEE Trans Biomed Eng 2009;56(12):2826–37.
[10] Nowak RD. Wavelet-based Rician noise removal for magnetic resonance
imaging. IEEE Trans Image Process 1999;8(10):1408–19.
[11] Castillo E, Morales DP, García A, Martínez-Martí F, Parrilla L, Palma AJ. Noise
suppression in ECG signals through efficient one-step wavelet processing
techniques. J Appl Math 2013;2013:763903.
[12] Wood JC, Johnson KM. Wavelet packet denoising of magnetic resonance images:
importance of Rician noise at low SNR. Magn Reson Med 1999;41(3):631–5.
[13] Coupé P, Yger P, Prima S, Hellier P, Kervrann C, Barillot C. An optimized blockwise
nonlocal means denoising filter for 3-D magnetic resonance images. IEEE Trans
Med Imaging 2008;27(4):425–41.
[14] Coupé P, Hellier P, Prima S, Kervrann C, Barillot C. 3D wavelet subbands mixing
for image denoising. Int J Biomed Imaging 2008;2008:590183.
[15] Anand CS, Sahambi JS. Wavelet domain non-linear filtering for MRI denoising.
Magn Reson Imaging 2010;28(6):842–61.
[16] Manjón JV, Coupé P, Concha L, Buades A, Collins DL, Robles M. Diffusion weighted
image denoising using overcomplete local PCA. PloS One 2013;8(9):e73021.
[17] Bao L, Robini M, Liu W, Zhu Y. Structure-adaptive sparse denoising for diffusiontensor MRI. Med Image Anal 2013;17(4):442–57.
[18] Manjón JV, Coupé P, Buades A, Louis Collins D, Robles M. New methods for
MRI denoising based on sparseness and self-similarity. Med Image Anal
2012;16(1):18–27.
[19] Hu J, Pu Y, Wu X, Zhang Y, Zhou J. Improved DCT-based nonlocal means filter for
MR images denoising. Comput Math Methods Med 2012;2012:232685.
[20] Gerig G, Kubler O, Kikinis R, Jolesz FA. Nonlinear anisotropic filtering of MRI data.
IEEE Trans Med Imaging 1992;11(2):221–32.
[21] Basu S, Fletcher T, Whitaker R. Rician noise removal in diffusion tensor MRI. Lect
Notes Comput Sci 2006;4190:117–25.
[22] Perona P, Malik J. Scale-space and edge detection using anisotropic diffusion.
IEEE Trans Pattern Anal Mach Intell 1990;12(7):629–39.
[23] Tsiotsios C, Petrou M. On the choice of the parameters for anisotropic diffusion in
image processing. Pattern Recogn 2013;46(5):1369–81.
[24] Tong CC, Sun Y, Payet N, Ong SH. A general strategy for anisotropic diffusion in
MR image denoising and enhancement. Magn Reson Imaging 2012;30
(10):1381–93.
[25] Krissian K, Aja-Fernández S. Noise-driven anisotropic diffusion filtering of MRI.
IEEE Trans Image Process 2009;18(10):2265–74.
[26] Yu YJ, Acton ST. Speckle reducing anisotropic diffusion. IEEE Trans Image Process
2002;11(11):1260–70.
[27] Golshan HM, Hasanzadeh RP, Yousefzadeh CS. An MRI denoising method using
image data redundancy and local SNR estimation. Magn Reson Imaging 2013;31
(7):1206–17.
[28] Buades A, Coll B, Morel JM. A review of image denoising algorithms, with a new
one. SIAM Multiscale Model Simul 2005;4(2):490–530.
[29] Buades A, Coll B, Morel JM. Image denoising methods. A new nonlocal principle.
SIAM Rev 2010;52(1):113–47.
[30] Manjón JV, Carbonell-Caballero J, Lull JJ, Garcia-Marti G, Robles M. MRI denoising
using non-local means. Med Image Anal 2008;12(4):514–23.
[31] Wiest-Daesslé N, Prima S, Coupe P, Morrissey SP, Barillot C. Rician noise removal
by non-local means filtering for low signal-to-noise ratio MRI: applications to
DT-MRI. Proc MICCAI; 2008. p. 171–9.
[32] He LL, Greenshields IR. A nonlocal maximum likelihood estimation method for
Rician noise reduction in MR images. IEEE Trans Med Imaging 2009;28
(2):165–72.
[33] Rajan J, Jeurissen B, Verhoye M, Van Audekerke J, Sijbers J. Maximum likelihood
estimation-based denoising of magnetic resonance images using restricted local
neighborhoods. Phys Med Biol 2011;56(16):5221–34.
[34] Maximov II, Farrher E, Grinberg F, Shah NJ. Spatially variable Rician noise in
magnetic resonance imaging. Med Image Anal 2012;16(2):536–48.
[35] Manjón JV, Coupé P, Marti-Bonmati L, Collins DL, Robles M. Adaptive non-local
means denoising of MR images with spatially varying noise levels. J Magn Reson
Imaging 2010;31(1):192–203.
[36] Rudin LI, Osher S, Fatemi E. Nonlinear total variation based noise removal
algorithms. Phys D 1992;60(1–4):259–68.
[37] Drapaca CS. A nonlinear total variation-based denoising method with two
regularization parameters. IEEE Trans Biomed Eng 2009;56(3):582–6.
[38] Wang Y, Zhou H. Total variation wavelet-based medical image denoising. Int
J Biomed Imaging 2006;2006:89095.
[39] Getreuer P, Tong M, Vese LA. A variational model for the restoration of MR images
corrupted by blur and Rician noise. Lect Notes Comput Sc 2011;6938:686–98.
[40] Martin A, Garamendi JF, Schiavi E. MRI TV-Rician denoising. Proc BIOSTEC; 2013.
p. 255–68.
[41] Valkonen T, Bredies K, Knoll F. Total generalized variation in diffusion tensor
imaging. SIAM J Imaging Sci 2013;6(1):487–525.
[42] Martin A, Schiavi E. Automatic total generalized variation-based DTI Rician
denoising. Lect Notes Comput Sc 2013;7950:581–8.
[43] Strong D, Chan T. Edge-preserving and scale-dependent properties of total
variation regularization. Inverse Probl 2003;19(6):165–87.
[44] Dong YQ, Hintermuller M, Rincon-Camacho MM. Automated regularization
parameter selection in multi-scale total variation models for image restoration.
J Math Imaging Vis 2011;40(1):82–104.
[45] Gilboa G, Sochen N, Zeevi YY. Variational denoising of partly textured images by
spatially varying constraints. IEEE Trans Image Process 2006;15(8):2281–9.
[46] Li F, Ng MK, Shen CM. Multiplicative noise removal with spatially varying
regularization parameters. SIAM J Imaging Sci 2010;3(1):1–20.
[47] Li F, Liu RH. Poisson noise removal with total variation regularization and local
fidelity. Chin J Electron 2012;21(2):304–8.
[48] Chen DQ, Cheng LZ. Spatially adapted total variation model to remove
multiplicative noise. IEEE Trans Image Process 2012;21(4):1650–62.
[49] Dietrich O, Raya JG, Reeder SB, Ingrisch M, Reiser MF, Schoenberg SO. Influence of
multichannel combination, parallel imaging and other reconstruction techniques on MRI noise characteristics. Magn Reson Imaging 2008;26(6):754–62.
[50] Koay CG, Ozarslan E, Pierpaoli C. Probabilistic identification and estimation of
noise (PIESNO): a self-consistent approach and its applications in MRI. J Magn
Reson 2009;199(1):94–103.
[51] Gudbjartsson H, Patz S. The Rician distribution of noisy MRI data. Magn Reson
Med 1995;34(6):910–4.
[52] Landman BA, Bazin PL, Smith SA, Prince JL. Robust estimation of spatially variable
noise fields. Magn Reson Med 2009;62(2):500–9.
[53] Landman BA, Bazin PL, Prince JL. Estimation and application of spatially variable
noise fields in diffusion tensor imaging. Magn Reson Imaging 2009;27(6):741–51.
[54] Coupé P, Manjon JV, Gedamu E, Arnold D, Robles M, Collins DL. Robust Rician
noise estimation for MR images. Med Image Anal 2010;14(4):483–93.
[55] Koay CG, Basser P. Analytically exact correction scheme for signal extraction
from noisy magnitude MR signals. J Magn Reson 2006;179(2):317–22.
[56] Shi J, Osher S. A nonlinear inverse scale space method for a convex multiplicative
noise model. SIAM J Imaging Sci 2008;1(3):294–321.
[57] Krishnan D, Fergus R. Fast image deconvolution using hyper-Laplacian priors.
Proc NIPS; 2009. p. 1033–41.
[58] Liu RW, Xu T. A robust alternating direction method for constrained hybrid
variational deblurring model. [arXiv, preprint] arXiv:1309.0123; 2013.
[59] Levin A, Fergus R, Durand F, Freeman WT. Image and depth from a conventional
camera with a coded aperture. ACM Trans Graphic 2007;26(3):1–9 [70].
[60] Collins DL, Zijdenbos AP, Kollokian V, Sled JG, Kabani NJ, Holmes CJ, et al. Design
and construction of a realistic digital brain phantom. IEEE Trans Med Imaging
1998;17(3):463–8.
[61] 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.
[62] Aja-Fernandez S, San Jose Estepar R, Alberola-Lopez C, Westin CF. Image quality
assessment based on local variance. Proc EMBS; 2006. p. 4053–6.
[63] Tristán-Vega A, García-Pérez V, Aja-Fernández S, Westin CF. Efficient and robust
nonlocal means denoising of MR data based on salient features matching.
Comput Methods Programs Biomed 2012;105(2):131–44.
[64] Fillard P, Souplet J, Toussaint N. Medical Image Navigation and Research Tool by
INRIA (MedINRIA), France; 2007.
[65] Pierpaoli C, Basser PJ. Toward a quantitative assessment of diffusion anisotropy.
Magn Reson Med 1996;36(6):893–906.
[66] Skare S, Li TQ, Nordell B, Ingvar M. Noise considerations in the determination of
diffusion tensor anisotropy. Magn Reson Imaging 2000;18(6):659–69.
[67] Shi L, Liu W, Zhang H, Xie Y, Wang D. A survey of GPU-based medical image
computing techniques. Quant Imaging Med Surg 2012;2(3):188–206.
[68] Pratx G, Xing L. GPU computing in medical physics: a review. Med Phys
2011;38(5):2685–97.
[69] Bui A, Cheng KT, Cong J, Vese L, Wang YC, Yuan B, et al. Platform characterization
for domain-specific computing. Proc ASP-DAC; 2012. p. 94–9.