Denoising and Visualization of HARDI Data
Tim McGraw, Evren Ozarslan, Baba Vemuri, Fellow, IEEE, Yunmei Chen, and Thomas Mareci,
AbstractDespite its apparent success, Diffusion Tensor MRI
(DTMRI) has significant shortcomings when the tissue of interest
has a complicated structure. This is due to the relatively simple
tensor model that assumes a unidirectional if not isotropic local
structure. As a more viable alternative Tuch et. al. have proposed
to do the data acquisition such that the diffusion sensitizing
gradient directions sample the surface of a sphere. In this high
angular resolution diffusion imaging (HARDI) method, one does
not have to be restricted to the tensor model and instead, it
is possible to calculate diffusion coefficients independently along
many directions. This imaging technique can reveal white matter
fiber crossings which would not be apparent in DT images.
In this paper, we present a novel variational formulation for
restoring the HARDI data and visualizing the fibers from
this restored data. This formulation involves smoothing signal
measurements over the spherical domain and across the 3D image
lattice. The smoothing on the spheres at each lattice point is
achieved using first and second order smoothness constraints, and
across the lattice via a total variation norm based scheme. For
the smoothing problem on the sphere, we use the finite element
method (FEM). Unlike the reported work on spherical harmonic
basis expansion of the diffusivity function on the sphere, the
FEM basis functions have local support and are therefore more
stable with respect to local perturbations caused by noise in
the data. To visualize the fiber paths, the probability values for
water molecules to move a particular distance along different
orientations were calculated using a Laplace series expansion of
these probabilities. We compute the Shannon as well as the Renyi
entropies of this distribution to characterize the anisotropy of
diffusion, with higher entropy corresponding to lower anisotropy.
These local distributions are also used to compute a vector
quantity called expected direction. Surfaces rendered using colors
corresponding to expected direction reveal anisotropy and fiber
direction in the imaged tissue. Further, examples are presented
to depict the performance of the HARDI data restoration and
visualization schemes on synthetic data, rat brain and spinal cord
data sets respectively.
Index Terms high angular resolution diffusion image, diffusion
tensor magnetic resonance imaging
I. INTRODUCTION
Fundamental advances in understanding living biological sys
tems require detailed knowledge of structural and functional
organization. This is particularly important in the nervous
system where anatomical connections determine the infor
mation pathways and how this information is processed. In
addition, understanding fundamental structural relationships is
essential to the development and application of therapies to
treat pathological conditions (e.g. disease or injury).
Observing the directional dependence of water diffusion can
allow us to infer structural information about the surrounding
tissue. White matter fiber bundles present a barrier to diffusion,
causing relatively high diffusivity along the fiber direction, and
lower diffusivity across the fiber.
Most imaging methods give only an anatomically isolated rep
resentation of living tissue because the images do not contain
connectivity information. Such information would allow the
identification and correlation of system elements responding
during function. For example in brain trauma, the relationship
between anatomy and behavior will only become apparent
when we are able to resolve the afferent nerve fiber pathways
transmitting the sensation from a stimulus to the brain or the
efferent pathways transmitting impulses from the brain area
controlling behavior.
MR measurements can be made sensitive to the translational
diffusion of water molecules by the utilization of magnetic
field gradients [1]. In general, the signal acquired depends on
the strength and the direction of these diffusion sensitizing
gradients. Repeated measurements of water diffusion in tissue
with varying gradient directions provide a means to quantify
the level of anisotropy as well as to determine the local
fiber orientation within the tissue. In a series of publications,
Basser and colleagues [2], [3], [4] have formulated a new
imaging modality called "diffusion tensor MRI (DTMRI)"
that employs a second order, positive definite, symmetric
diffusion tensor to represent the local tissue structure. They
have proposed several rotationally invariant scalar indices
that quantify different aspects of water diffusion observed in
tissue, similar to different "stains" used in histological studies
[5]. Under the hypothesis that the direction along which
the diffusion coefficient is largest will yield the local fiber
orientation, one can determine the directionality of neuronal
fiber bundles. This fact has been exploited to generate fiber
tract maps that yield information on structural connections in
human [4], [6], [7], [8] as well as rat brains [9], [10], [11],
[12], [13] and spinal cords [14].
Fig. 1. Orientational heterogeneity in DTI (left), and HARDI (right).
Despite its apparent success, DTMRI has significant short
comings when the tissue of interest has a complicated geom
etry. This is due to the relatively simple tensor model that
assumes a unidirectional if not isotropic local structure. In
the case of orientational heterogeneity, DTMRI technique is
likely to yield incorrect fiber directions, and artificially low
anisotropy values. This is due to violation of the assumption
of Gaussian probability model characterizing the diffusion
implicit in DTI. In order to overcome these difficulties several
approaches have been taken. Qspace imaging, a technique
commonly used to examine porous structures [15], has been
suggested as a possible solution [16]. However this scheme
requires strong gradient strengths and long acquisition times
[17], or significant reduction in the resolution of the images.
As a more viable alternative Tuch et. al. have proposed to
do the acquisition such that the diffusion sensitizing gradients
sample the surface of a sphere [18], [19]. In this high angular
resolution diffusion imaging (HARDI) method, one does not
have to be restricted to the tensor model and instead, it
is possible to calculate diffusion coefficients independently
along many directions. This method does not require more
powerful hardware systems than that required by DTMRI.
Several groups have already performed HARDI acquisitions
in clinical settings and have reported 43 to 126 different
diffusion weighted images acquired in 20 to 40 minutes of
total scanning time [20], [19], [21] indicating the feasibility
of the high angular resolution scheme as a clinical diagnostic
tool. In Figure 1, we present a matrix of voxels from a rat
brain data set showing, renderings of DTIbased estimates
of orientation and HARDIbased orientation estimates. The
orientation heterogeneity is evident from the HARDIbased
renderings at each voxel since HARDI measurements can
resolve multiple dominant directions of anisotropy at a voxel,
a feature which is lacking in the DTI. Since the HARDI data
acquisition is very nascent, not many techniques of processing
the HARDI data have been reported in literature. In the
following, we will briefly review the few very recently reported
techniques of HARDI data denoising, which must be done
prior to further analysis or visualization.
A. Restoration
Processing of HARDI data sets has received increased atten
tion lately and a few researchers have reported their results in
literature. The use of spherical harmonic expansions have been
quite popular in this context since the HARDI data primarily
consists of scalar signal measurements on a sphere located
at each lattice point on a 3D image grid. Tuch [18], [19]
developed the HARDI acquisition and processing and later
Frank [20] showed that it is possible to use the spherical
harmonics expansion of the HARDI data to characterize the
local geometry of the diffusivity profiles. In his work however,
there is no discussion of denoising/restoring the HARDI signal
measurements, which is essential for subsequent processing
and interpretation. Chen et al. [22] find a regularized spherical
harmonic expansion by solving a constrained minimization
problem. However the expansion is a truncated spherical
harmonic expansion of order four, and hence the solution can
at best only represent two fiber directions within a voxel. In
[21], Jansons and Alexander described a new statistic which
was called persistent angular structure that was obtained from
the samples of a 3D function, in this case the displacement
of water molecules in each direction. The goal in their work
was to resolve voxels containing one or more fibers. However,
there was no discussion on how to restore the noisy HARDI
data prior to resolution of the fiber paths.
In contrast to HARDI denoising, DTMRI denoising has been
more popular and numerous techniques exist in literature. Most
of the techniques in literature have employed vectorvalued
(See [23]) smoothing to tensor field data by using spectral
decomposition of the tensors and smoothing the eigen vectors
and eigen values using techniques introduced in [24]. Recent
work by Tschumperle et al., [25] deals with smoothing the
eigenvector field of the diffusion tensors computed from the
raw echo intensity image data. Several other methods have
been developed for restoring the DTMRI data sets but most
of them use existing restoration schemes for scalar or vector
valued functions from image processing literature. More re
cently, matrix valued function restoration was introduced in
[26] and applied to the restoration of noisy diffusion tensor
fields. An interesting alternative to the variational principle
approach was taken by Weickert and Brox [27] wherein,
they developed an anisotropic diffusion filter that achieves
an "edge" preserving smoothing of the positive semidefinite
tensorvalued image. All of the matrixvalued smoothing meth
ods share one commonality i.e., they use a linearized Stejskal
Tanner equation as the data acquisition model which does
not accurately reflect the physics of the data acquisition. In
recently reported work, by Wang et al., [28], [29], an alter
native approach which overcomes this weakness by directly
estimating a smooth positive semidefinite tensor field from the
raw data using the actual nonlinear monoexponential model
characterized by the StejskalTanner equation [1] was reported.
B. Tractography
Water in the brain preferentially diffuses along white mat
ter fibers. By tracking the direction of fastest diffusion, as
measured by MRI, noninvasive fiber tracking of the brain
can be accomplished. In the context of DTI data, fiber tracks
estimated in reported literature were obtained by repeatedly
stepping in the direction of fastest diffusion. The direction
along which the diffusion is dominant corresponds to the
direction of eigenvector corresponding to the largest eigen
value of the tensor D. This approach was taken by many
researchers [30], [9], [31], [32], [23], [33], [14], [13]. Most of
these techniques do incorporate some regularization in their
stream line estimation schemes in order to generate the fiber
pathways. Techniques that are quite distinct from the idea of
stream line generation have also been reported in literature.
Batchelor et. al., [34] reported a fiber tract mapping scheme
wherein, they produce a map indicating the probability of a
fiber passing through each location in the field. However, no
discussion on how to estimate the actual fibers was described.
An alternative approach based on sequential importance sam
pling and regularization techniques was proposed in [35],
which allowed paths to originate from a single location and
branch out and produced a probability distribution of the
paths. O'Donnell et. al., [36] describe ways to estimate the
connectivity from the given tensor field. One approach they
suggested was to estimate the geodesics in the locally warped
space where the warping is derived from the local tensor.
Behrens et al. [37] estimate the local probability distribution
of fiber directions at each voxel and proceed to compute global
connectivity from these local PDFs. Unlike previous work, the
model parameters are described using distributions, and the
PDF at each voxel are estimated using Markov Chain Monte
Carlo techniques.
In the context of HARDI processing, ,..... ,i,.,,0 for fiber
path computation have been reported in [20], [21], [38]. Frank
[20] described a spherical harmonic transform representation
of the HARDI data and pointed out that, due to the antipodal
symmetry on the sphere, only the even terms in the spherical
harmonic transform contributed toward the representation. Any
nonzero odd terms would be due to artifacts in imaging
such as noise. The fiber pathways were not computed in
his reported work, however it was suggested that one could
estimate the direction of the crossing fibers by using a multi
tensor model and finding the principal diffusion directions of
this multitensor expansion. The approach of using HARDI data
and then reverting to a multitensor approach seems somewhat
awkward. The approach described by Ozarslan et al. [38]
expresses the diffusivity function, a function defined on the
sphere, as a generalized (higher rank) Cartesian tensor and
then estimates the probability distribution of water molecule
displacement over all directions using the FFT (fast Fourier
transform) of the signal measurements either on the sphere or
an interpolated Cartesian grid. These distribution profiles are
then sharpened to display the possible orientation of the fibers
with complex local geometries. In this paper we will not be
(., iijt ' i the fiber tracts but will present ways of visualizing
the probability maps of water di/ttr.\in via color displays of
measures of anisotropy.
C. Anisotropy
In DTI there are many scalar measures which characterize
the anisotropy of the diffusion phenomenon. Most useful are
the rotationally invariant measures, those independent of the
laboratory reference frame, such as volume ratio, rational
anisotropy and fractional anisotropy [5]. These can be com
puted in terms of the eigenvalues of the diffusion tensor. These
scalar indices alone have been useful in clinical studies [39],
so it is important to find analogous measures of anisotropy
for HARDI. The measures for HARDI should overcome the
weaknesses of the tensor model in regards to regions of
crossing fibers. In the case of orientational heterogeneity, DT
MRI technique yields incorrect fiber directions, and artificially
low anisotropy values [40].
Frank [20] used the coefficients of the spherical harmonic
expansion of the apparent diffusion coefficient to quantify
anisotropy. The Oth, 2nd and 4th order coefficients describe
isotropic diffusion, singlefiber diffusion and twofiber diffu
sion respectively. Voxels can be classified using the relative
magnitudes of these coefficients. In [38] Ozarslan et al. use
a high rank tensor to describe the diffusion process and
generalize the rank2 tensor trace to higher rank tensors. The
variance of the diffusivity is then expressed in terms of the
generalized trace, and a transfer function maps this value to the
[0, 1] range. This quantity is called "generalized anisotropy"
(GA). It was shown to be useful in detecting anisotropy at
fiber crossings which are not detectable using FA.
D. Overview of Our Modeling Scheme
Our goal here is to denoise the HARDI data, develop a fast
computation scheme for estimating the molecular displace
ment distribution, and compute a novel index of anisotropy
from these distributions along with a visualization scheme for
the same. We present a novel and effective variational formu
lation that will directly estimate a smooth signal S(O, 0) and
the probability distribution of the water molecule displacement
over all directions p(O, 0), given the noisy measurement
S(0, ) = Soexp(bd(O, 4))+ (0, 0) (1)
where S is the signal measurement taken on a sphere of
constant gradient magnitude over all (0, 0), b is the diffusion
weighting factor, d(O, 0) is the diffusivity as a function of the
direction expressed by the elevation and azimuth angles on
the sphere and r1(0, 0) is the Gaussian noise. The variational
principle involves smoothing S values over the sphere and
across the 3D image lattice. The key factor that complicates
this problem is that the domain of the data at each voxel in the
voxel lattice is a sphere. One may use the levelset techniques
developed by Tang et al., [41] to achieve this smoothing
however, when data sets are large, it becomes computationally
impractical to apply the levelset technique at each voxel
independently to restore these scalarvalued measurements on
the sphere. We arrive at a computationally efficient solution
to this problem by using the finite element method (FEM)
on the sphere and choosing local basis functions for the data
restoration. Unlike the reported work on spherical harmonic
basis expansion of the diffusivity function on the sphere [42],
[38], [22], the FEM basis functions have local support and are
more stable to perturbations due to noise in the data.
From the denoised data we will compute a probability,
pt(0, 0), of molecular diffusion over a sphere of directions.
The Shannon entropy as well as Renyi entropy of this distri
bution will then be used to quantify anisotropy.
II. METHODS
The HARDI processing proceeds by acquiring diffusion
weighted images with many diffusion encoding gradient di
rections, effectively sampling a spherical shell of the qspace.
It is desired that this sampling be uniform, or nearly so.
The gradient direction for each image is usually chosen to
correspond to the vertices of an icosahedron which has been
repeatedly subdivided. Since the process of diffusion is known
to be symmetric, we need only sample one hemisphere of
the qspace. In the case of our data, we consider 81 or 46
gradient directions. The gradient directions correspond to the
vertices of the surface shown in Figure (2). In addition a low
bvalue (small diffusion encoding gradient magnitude) image
is acquired in order to estimate So.
The random process of diffusion of water molecules is de
scribed by the diffusion displacement PDF pt(r). This is the
probability that a given molecule has a diffusion displacement
of r after time t. The relation between the measured image,
and the diffusion displacement PDF is given by Callaghan in
[15] as
t(r)j S( exp(2riq. r)dq (2)
so
where So is the image acquired with no diffusion encoding
gradient applied. This is simply the Fourier transform of s(q)
It is the modes of pt(r) that are taken to be the underlying
fiber directions.
Fig. 2. HARD gradient directions (0,4) correspond to the vertices of a
subdivided icosahedron.
III. VARIATIONAL PRINCIPLE FORMULATION
We propose a membranespline deformation energy mini
mization for smoothing the measured image S(x, 0, 0). The
variational principle for estimating a smooth S(x, 0, 0) is given
by
minE(S)= II S S(x, 0,) S(x, 0,)12dSdx
SjV(e,0)Sl2dS+ g(x)V Sl. (3)
where Q is the domain of the image lattice and S2 is the
sphere on which the signal measurements are specified at
each voxel. The first term of Equation (3) is a data fidelity
term which makes the solution to be close to the given data.
The degree of data fidelity can be controlled by the input
parameter p. The second term is a regularization constraint
enforcing smoothness of the data over the spherical domain at
each voxel. The third term is another regularization term which
causes the solution to be smooth over the spatial domain (the
3D voxel lattice).
A. Finite Element ,u.'. .i rii of S(O, )
We will consider a deformation energy functional which
is a weighted combination of thinplate spline energy and
membrane spline energy, that is commonly used in computer
vision literature for smoothing scalarvalued data in 913 (see
McInerney and Terzopoulos [43]).
The diffusionencoding gradient directions are taken as the
vertices of a subdivided icosahedron, to acheive a nearly uni
form sampling of spherical direction. We map this piecewise
linear approximation of a sphere to the plane by taking the
spherical coordinates (0, 0) of the imaging gradient direction
as the planar global coordinates (u, v). The gradient directions,
and their embedding in the plane are shown in Figures (2) and
(3).
Fig. 3. 2 dimensional (u,v) domain for FEMbased smoothing of HARDI
data.
Note that, the data can be seen as height values defined on the
sphere at the coordinates specified parametrically by (u,v).
The smoothing will be applied to these height measurements
z(u,v) using the following smoothing functional.
/' /i(a(zu2+ zv2) +/3(zu12 +2zv12 + 2))dudv
(4)
The weight on the membrane term is a and the weight on
the thinplate term is /3. Once we have computed a smooth
z(u,v), the result will then be mapped back to the image on
the sphere, S(O, 4).
The data energy due to virtual work of the data forces is
,= z(u,v)f(u,v)dudv (5)
The restoration at each voxel is formulated at the energy
minimization
min (S) = p(S)+ Sd(S)
S
We use polynomial shape functions, N, as a basis for the data
over the domain of spherical directions.
z(u, v) qN, (u, v) = Nq (7)
Where the N is a (1 x n) row vector, and q is a column vector
of nodal variables.
1) Element Matrices: The domain is subdivided into ele
ments, each with their own local shape functions. For each
element j, we have,
z(u,v) NJ (u,v)q (8)
........, ,,
for (u,v) E Q,. The local energy function for each element is
given by
:= (azi12+a l2+ Ul12+2Plz,2+ lz12)dudv
(9)
Then, the global energy is the sum of the energies of each
element,
(10)
The element strain vector (given by Dhatt and Touzot [44]) is,
(N1),
(N1)v
(N1),
(Na>
qJ Bq
(11)
(12)
TABLE I
GAUSSRADAUWEIGHTS
1 r, wi, a,
1 0.0469100770 0.1184634425 0.0398098571 0.1007941926
2 0.2307653449 0.2393143353 0.1980134179 0.2084506672
3 0.5 0.2844444444 0.4379748102 0.2604633916
4 0.7692346551 0.2393143353 0.6954642734 0.2426935942
5 0.9530899230 0.1184634425 0.9014649142 0.1598203766
2) Local Element Coordinates: We now present the coor
dinate system for the local elements. For local elements,
triangular elements are used with a barycentric coordinate
system (y, r) so that each coordinate is in [0,1] and y
15r.
S0=0
Fig. 4. Mapping to barycentric coordinates
(13)
(14)
where we have defined B as the (n x 5) matrix of derivatives
of the nodal basis functions. We can then write the element
strain energy as:
JJE"JDEJdudv
where we define
a 0 0 0 0
0 a 0 0 0
D = 0 0o 0 0
0 0 0 2 0
0 0 0 0 3
Since qJ is constant over each element we can derive the
element stiffness matrix in terms of D and B as follows:
= J qJTBJTDBJqJdudv= qJTK'q'. (15)
We will model the data constraint as springs pulling z(u,v)
toward some given value zo(u, v). The force at each point will
obey f k(z zo), where k is the spring constant.
NJ'qk(N'q zo)dudv (16)
kqJ' j N N'Jqdudv+kqJi NJTzodudv
(17)
FJ and fJ are defined such that the first term of Equation (17)
is qIFJFqJ and the second term is qTfj. We can then balance
deformation energy and data energy by solving the following
linear system:
(KJ +FJ)qJ f. (18)
u(n,r7) =(1 l)uo +U l+u2
v(, 7r) (1  )vo V1+ V2
u 1 l
vj [ v1
(19)
uo u2 o 0 ] [
VO V2 V ]o
(20)
(21)
Su1 au u
du du d
dvJ v v d J
Integrals over the (u, v) domain can be converted to integrals
over the local (, n) domain in the following way:
f f(u,v)dudv f f(u(dv),v( ,))det(J)d dn
(22)
Using GaussRadau quadrature rules, we can approximate the
integral in Equation (22) by
Swi,wjjf(u(, ,,j), v(g,, n,,,)) det(J) (23)
1j 1
where ,,, = r,(1 sj), wjj, a,(1 ), j, and wi, are given
in Table I. Derivatives over (u,v) become
dN dN d N dcr
dcu 9 cu c rd cu
dN dN d dN d2
dv d 9v d 9 (24)
C
I
The partial derivatives of and rn with respect to u and v can
be computed by inverting the Jacobian
dn
det(J) [
Sdu i du
dv I dv
V2 V (u22 o)
(vi o) u uo
which we can relate to each other by
r
(25) qu,v
0
0
0
(furlv + rlv)
2ov
0
0
0
2
lu
rlurlv
2
nV
(29)
(26) 3) Global Matrices: We now wish to construct global matri
ces so that the energy balance over the entire FEM mesh is
given by the linear system
We use the fifth order element shape functions given by Dhatt
and Touzot [44]. This element guarantees C1 (surface normal)
continuity across triangles. The basis functions are given by:
N= 2(10A 152 + 613 +301q( + 1))
N2 2 (3 2 3 +641)
N3 = qA2(32 3132 +6)
1
N4 242(1 +271)
AN5 (L2
1
NA6 12x2(1 +24 )
2
N7 = 2(104 152 +6 3+1512)
N8 2(8 +142 63 15 12)
1
N9 = 42 (6 44 31 312 +341)
2
NO 1 2(2 4(1 _)2 +512)
4
=4 1 ^2(634 41 3 +3{1)
1
N15 2( + 142 + + 12 541)
2
12 2+ 13 2
N12  271 4 q*
4 2
N13 12q2(101 152 +63+15 12L)
I
N114 1 _12 (6 34 41 3 2 +34q)
2
N16 127 12x q 2 q2
N17 1412(2 + (1 ) +2 25_2 )
2
4112(1
Kq f
(30)
where K is a (6n x 6n) matrix since we have 6 variables per
node.
We will consider the simple case of 2 elements. Expanding the
element Equation(18) in terms of nodal variables for element
0, we get
K o[
Ko2, 0
K,2
K01,2
KO,2
KO2,2
(31)
and likewise for element 1
K3,3
K1,3
K,2
K1
1,2
(32)
K1
1,1
where each qJ is a (6 x 1) column vector of nodal variables.
We expand each KJ to be (6n x 6n) by inserting rows and
columns of zeros corresponding to each node of the mesh.
Also expand f' to (6n x 1). The global K and q are obtained
by summing the expanded matrices from each element in the
mesh. For our 2 element example we have
0,0 Kj,1 K ,2 0 q0 1
'1 0
K0 KI0 +1 l K,2 + K1,2 K11,3 q f +fl
K02' K2' + K2,l KO, + K 2 K2 3 q1 1+
0 K, 1 1K q i
) K3, K3 K3 3 q
(33)
The global linear system Equation (30) is symmetric, and has
a sparse banded structure with 18 nonzero diagonal bands. An
efficient solution to q is obtained via Cholesky factorization.
(27) B. Spatial ,,,. *. ri,,, of S(x)
The quintic shape functions have nodal variables which can
be written in terms of local or global coordinates as,
Z
z
Zr q
Z41
Zrjr
Z111
z
Zu
Zv
Zuu
Zuv
Zvv
Smoothing the raw vectorvalued data, S(x), is posed as a
variational principle involving a first order smoothness con
straint on the solution to the smoothing problem. Note that
the data at each voxel is a large set of S measurements over a
sphere of directions and can be assembled into a vector after
the smoothing on the spherical coordinate domain has been
accomplished. We propose a weighted TVnorm minimization
for smoothing this vectorvalued image S. The variational
principle for estimating a smooth S(x) is given by
im 11 m
m ,,S) /(g(x) VSl+ P I(S 2)dx
s Q 1=1 2 ,=1
(34)
where, Q is the image domain and p is a regularization factor.
The first term here is the regularization constraint on the solu
tion to have a certain degree of smoothness. The second term
in the variational principle makes the solution faithful to the
data to a certain degree. We have used the coupling function
g(x) 1/ (1+ GA (x)) for smoothing HARDI, where GA is the
generalized anisotropy defined in Ozarslan et al., [40] and is
computed from the variance of normalized diffusivity. For a
more detailed discussion on GA, we refer the reader to [40].
This selection criterion preserves locations of dominant GA
while smoothing the rest of the data. Note that since we are
interested in regions of dominant anisotropy, it is apt to choose
such a selective term.
Here we have used a different TVnorm than the one used by
Blomgren and Chan [45]. The TVn,m norm is an L2 norm of
the vector of TVn,1 norms (fQ VS,(,> 1. / I for each channel.
We use the L1 norm instead, which is known to have better
discontinuity preservation properties.
The gradient descent form of the above minimization is given
by
as'
dt
ds,
dn
div VS, L) (us I S) i=1,... ,m
^ II }V 'I
2,
0 and S(x,t 0) S(x)
The use of a modified TVnorm in equation (34) results in a
looser coupling between channels than when using the TVn,m
norm. This reduces the numerical complexity of Equation (35)
and makes solution for large data sets feasible.
The gradient descent of the vectorvalued image smoothing
using the TVn,mnorm TVn,m(S(x)) XI[TV, 1(S,)]2 pre
sented in [45] is given by,
dS,(x,t) TV, 1(S,) VS,
dt TV,,.(S) IIVS,I
S(x,0) So(x).
(36)
Note that the TVn,m norm appears in the gradient descent so
lution of the vectorvalued minimization problem. Considering
that our data sets consist of up to 82 images, corresponding
to (magnetic field) gradient directions, calculating the TVn,m
norm by numerically integrating over the 3dimensional data
set at each step of an iterative process would be prohibitively
expensive. In contrast using our modified TVnorm suggested
earlier leads to a more efficient solution.
1) FixedPoint LaggedD,im,,, ,i Since the m Equations(35)
are coupled only through the function g, we can drop the
subscript on S with no ambiguity (later the subscript will
refer to spatial coordinates.) In this section we will discuss
the numerical solution for each channel, S, of the vector
valued image S. Equation (35) is nonlinear due to the presence
of IIVS, in the denominator of the first term. We linearize
Equation (35) by using the method of "laggeddiffusivity"
presented by Chan and Mulet [46]. By considering IIVSI to
be a constant for each iteration, and using the value from the
previous iteration we can instead solve
1
VS (Vg. VSt +gV2St+l)+ (St+l S) 0 (37)
Here the superscript denotes iteration number. First, rewrite
Equation (37) with all of the S4+1 terms on the lefthand side
_g2s + P IIVlII S 1 VS + Vg. V
g g
(38)
2) Discretized Equations: To write Equation (38) as a linear
system (ASt+l f), discretize the Laplacian and gradient
terms. Using central differences for the Laplacian we have
V2St1+
1t+l it c+1 i 1^+l
x 1,y,z x ,y 1,z x,y,z 1
69 + t+1 + 1 S 1 (39)
,z x+l,y,z x,y+ l,z x,y,z+
Define the standard central differences to be
1
1
AyS 2 (Sxy+1z Sx,y1,z)
AzS = (Sxy+l Sxy,z 1)
(40)
We can rewrite Equation (38) in discrete form using the
definitions in Equation (40)
(35)
Sx 1,y,z Sx,y l,z Sx,y,z 1
Vg/( A S ((S '
+(6+ / )2+(ASg)2(Az)2 )Sx,y,z
x+l,y,z x,y+ 1,z x,y,+ 1
1 (S v(AxSt)2 + (ASt)2 + (Az)2
+AgAxS + AygAyS + AzgAzS)
This results in a sparse banded linear system with 7 nonzero
coefficients per row.
6, 1ii1 0
0 1 1 1 0 g2
0+6+ 1 1
n~3
(42)
where the righthand side of Equation (41) has been replaced
with A. The matrix in Equation (42) is symmetric and diago
nally dominant. We have successfully used conjugate gradient
descent to solve this system.
The solution of Equation (42) represents one fixedpoint
iteration. This iteration is continued until IS' S+11 < c, where
c is a small constant.
C. C ,,j,ia, Probabilities
The probability in Equation 2 can now be evaluated by
computing the quantity S(q)/So and performing the FFT Since
we know that the signal, S, decays exponentially from the
origin of qspace (where S(0) = So) we can interpolate signal
values for arbitrary q. We resample from spherical coordinates
to cartesian and perform the FFT on the resampled data. The
result is a probability of water molecule displacement over a
small time constant. We are interested in only the direction of
(41)
water displacement, so we integrate out the radial component
of pt(r) to get pt(O, 0). This is commonly referred to as the
diffusion orientation distribution function (ODF). Computing
the ODF with this method is computationally expensive since
it requires a 3D FFT at each voxel, and then a numerical
integration for each direction. In this section we present an
alternative method which makes computing the ODF for large
datasets extremely efficient.
The Fourier transform that relates the signal attenuation to
the water displacement probability (Eq. 2) can be written in
spherical coordinates. This is a consequence of the pointwise
convergent expansion of the plane wave in spherical coordi
nates [47] given by
e 2'iq.r = 41 r (i)'lj(2rqr)Y1,(O, J)*Ylm,(O, 0,)
I=0m=l
(43)
where ji is the lth order spherical Bessel function, Yim is the
spherical harmonics and Or and 0r are the spherical coordinates
depicting the direction of the vector r. Inserting this expression
into Eq. 2 and setting r = ro, we get
S=0m=
pt(ro, Ori r) S I (i)I1lYi(Or,4Jr) jsYi(0, 4)*Il(0, 4) dS
(44)
where 1/(0, 0) are evaluated from the radial part of the integral
in qspace. If we assume that along each direction in qspace
the signal values attenuate according to Equation (1), then the
radial integrals yield
Sr F(+3)
l 2+3 7 3/2 (D(0, O)t)(l+3)/2 (l+ 3/2)
F 1 1+3 ; 3 r ) (45)
F 2 2' 4D(O, )t
where 1F1 is the confluent hypergeometric function [48].
Note that the function pt(ro, Or, Or) is the probability of the
water molecule to move from the origin to a point ro away
from the origin along the direction specified by the spherical
coordinates Or and 0r, i.e., we will be interested in the
probability values on a sphere of radius ro. In order to estimate
the probability on the surface of a sphere of radius ro, we go
back to Eqs. 44 and 45, and expand 1l(0, 0) in Laplace series,
i.e.,
where
I'
T(o,0), E aIM e (
'=0m'=I'
al', Y"m(O(0,0)*ii(0,0)dS
which is just a Laplace series expansion of pt(ro, 0, 0). Note
that coefficients of this Laplace series for some 1 value come
from the lth order spherical harmonic transform of I (0, 0).
In summary, the estimation of the angular probability dis
tribution at a distance ro away from the origin involves the
following steps:
Given HARDI data including NG images acquired with
gradient directions along different orientations, calculate
the diffusivity D(O, 0) along each direction.
Calculate 11(0, 0) using Eq. 45.
Calculate allm, the lth order spherical harmonic trans
form of 1/(0, 0) for each 1.
Evaluate Eq. 48, which is just
Pt(ro, 0r, Or)
tkooo Cx22m42m(Orr)
=000 2
= a 22m 2m(nO,,r)
4
+ I a44mY4nm(Or, r)
m= 4
6
Y 0666mnY6nm(0r, r) +
m= 6
(49)
Note that only even order 1 terms are retained which is
a consequence of the antipodal symmetry of D(O, 0). This
scheme provides a fast way to calculate the orientation profiles.
In our implementation we have evaluated the series given in
Eq. 49 up to 1 6 terms, and ro was set to 17.5pm.
D. Visualizing Probabilities
In diffusion tensor image processing, the scalar quantity
known as fractional anisotropy (FA) is often considered a
useful quantity to visualize. An example of an FA image
is shown in Figure (5). The FA value at each voxel can be
computed from the eigenvalues of the tensor at that voxel. The
values of FA range from 0, for completely isotropic diffusion,
to 1, for completely anisotropic diffusion. The bright regions
of Figure (5) correspond to high FA values. In rank2 tensor
images FA can indicate the presence of a single fiber direction
within a voxel.
(46)
Comparing the integration over the sphere in Eq. 44 with the
expression in Eq. 47, it is easy to see that
pt(ro, Or, Or)= (i)'alilIml(Or, Or) (48)
l=0m=
Fig. 5. FA image (left), generalized anisotropy (center), Shannon anisotropy
(HA) (right), from coronal slice of rat brain.
Since we have a distribution at each voxel, we can compute
the Shannon entropy value at each voxel as given by,
H(p) pi(O, ,)logp(O, ). (50)
I=1
Considering the entropy of several example distributions,
we can get a feel for the interpretation of entropy in the
context of HARDI. Entropy attains its maximum value for a
uniform distribution. In our case, this corresponds to isotropic
diffusion. The entropy of a Gaussian distribution decreases
as the variance decreases. A voxel with this distribution has
oriented diffusion greatest in the direction of the mode of the
distribution, implying the presence of fibrous structure. Tuch
[49] defined a scalar measure of anisotropy called "normalized
entropy", which was also based on Shannon entropy.
The anisotropy images we present are mapped using Equation
(51) such that the black color is high entropy isotropicc dif
fusion) and higher intensity grey colors represent low entropy
(anisotropic diffusion). We denote the anisotropy measure
computed from Shannon entropy as HA.
HA(p) = (1.0+ H ) (51)
log n
This allows the image to be interpreted in the same way as a
fractional anisotropy image where the white color corresponds
to white matter, grey corresponds to grey matter and black
corresponds to cerebrospinal fluid.
The HA index is not a generalization of FA, and their values
cannot be compared in a meaningful way. In order to highlight
the difference in anisotropy measures for HARDI and DTI,
a Gaussian ODF was computed from the tensor data, and
the Shannon entropy of these distributions was computed.
The anisotropy values for the Gaussian ODFs can then be
compared to the anisotropy values for the true ODFs. The
difference image is shown in Figure (6). It can be seen that
there is a structure to the difference image. This ditlierenl e is
more pronounced in regions where the tensor model predicts
diltti.\ion, is more isotropic than it actually is.
Fig. 6. Shannon anisotropy difference between Gaussian ODF and general
ODF.
Shannon entropy is not the only measure of uncertainty in
a random variable. Another definition of entropy, called the
R6nyi entropy can be seen as a generalization of Shannon
entropy, [50]. R6nyi entropy of order a is given by
Ha.(= log(p(O,,,)a) (52)
1 a =1
The parameter, a, of this entropy formulation has several
interesting properties:
lima, Ha(p) =H(p) (Shannon entropy)
Ho(p) = number of nonempty bins in the histogram of
P.
H_(p)= log(max,(p,)) (depends only on the mode of
P)
One may interpret the order, a, as a parameter which changes
the shape of the distribution, p. For a > 1, small values of
p(x) will shrink closer to zero. As a increases, these small
probabilities approach zero more quickly. For high a, events
with high probability influence the entropy more. This can
be seen as controlling the contrast between white and gray
matter. We can formulate an anisotropy index based on R6nyi
entropies just as we did for Shannon entropy. Anisotropy
images computed from R6nyi entropies of different orders are
presented in Figure (7).
Fig. 7. Anisotropy computed from H2 (topleft), H5 (topright), Hlo (bottom
left), H2o (bottomright).
In statistical physics, a quantity called structural entropy has
useful physical interpretation [51]. This entropy is computed as
the difference H\ H2. It is unclear whether there is any such
physical meaning for this quantity in the context of diffusion,
but it does provide a means of applying a transfer function to
entropy images which has a firm information theoretic footing.
We present entropy difference images for several values of a
in Figure (8).
The scalar entropy value has no directional information how
ever. In DTI, there is precedence for using color values to
represent direction. Color FA images are a mapping of the
principal diffusion direction (the dominant eigenvector of the
tensor) to a hue, and the fractional anisotropy value to an
intensity. These images are useful for distinguishing between
adjacent anisotropic regions which differ in direction. Since we
may have high diffusion coefficients in several directions, we
integrate over the sphere to determine a representative color
for each voxel. We denote this value as ED, for "expected
direction".
Fig. 10. ED colors for synthetic distribution field.
Fig. 8. R6nyi entropy differences H1 H2 (topleft), H1 H5 (topright),
H1 Hlo(bottom left), H1 H20 (bottomright).
ED
1=1
cos 0, sin 01
sin0,sin (i (p(O,,1)) Pm1n)
cos 0I
The resulting vector is interpreted as red, green, and blue color
components. The directions with the highest probability of
diffusion will have their corresponding color contribute most
to the resulting color. Figure (10) depicts the colors computed
for a synthetic ODF field directed in a curve. Figure (11) shows
the colors computed for ODFs oriented in various directions,
both with and without crossings.

Fig. 9. Original ODF
sharpened ODF (right).
(left), Minimum probability sphere (center), and
To enhance the visual impact of the ODF we apply a sharpen
ing transform to the ODF by subtracting a uniform distribution
(sphere) from each distribution, as shown in Figure (9). The
radius of the sphere is the minimum of the probability over
all directions. By performing this operation the direction of
maximum probability becomes more apparent. This causes
regions of isotropic diffusion to have lower intensity in the
color expected direction images.
IV. EXPERIMENTAL RESULTS
The denoising and rendering techniques described in the pre
vious section were first applied to a synthetic HARDI dataset.
This dataset was generated using the technique described by
Ozarslan et al. in [40]. The dataset was designed to depict
a region of curving fibers, a region of straight fibers, and
Fig. 11. ED colors for synthetic distribution field.
a crossing between the two. Representative sharpened ODF
profile are shown overlaid on the HA image on the left side
of Figure (12). The expected direction visualization shown on
the right side of Figure (12).
Y AI
Fig. 12. ODF profiles overlaid on HA image (left), and ED (right) of the
synthetic dataset.
A small sample of the synthetic data, taken from near the
crossing region, is shown in Figure (13). The synthetic data
was corrupted with Gaussian noise of mean zero, and variance
a2 0.005. The noisy data is shown in Figure (14). The same
voxels are shown after smoothing over the spherical manifold
in Figure (15), after smoothing over the cartesian image
domain in Figure (16) and after both techniques have been
used in Figure (17). The righthand side of each figure shows
the sharpened ODF computed from the S values depicted on
the lefthand side.
' t. L
i C
1: ~
^IQ & W)'
L.
K\\N fcC
9.: JC
Fig. 13. Synthetic S (left), resulting ODF (right).
lb I
tv. ?4 k
tcasi '7
I 'i )I
.^,^ ^
< '
i
'.E
, ,
Fig. 14. Synthetic S with noise added (left), resulting ODF (right).
Fig. 17. Manifold and lattice smoothed S (left), resulting ODF (right).
TABLE II
ERROR BETWEEN GROUNDTRUTH AND RESTORED SYNTHETIC DATA
Method
p = No Restoration
p = FEM Restoration
p = TV Restoration
p = FEM + TV Restoration
I. '))
1.0409
0.9088
0.7420
0.6576
0.0173
0.0125
0.0119
0.0139
variational formulation, the spikes of noise present in the raw
data have been smoothed while preserving the overall shape
S the S profile. This smoothness is propagated to the computed
S ODF profiles.
We can compare the resulting distributions with the ground
truth by using the square root of Jdivergence (symmetric KL
divergence) as a measure. This divergence is defined as
d(p, q)= J(pq)
t. I"
(54)
where
J(p, q) = p(O,, ,) log + qP(,, 0,) log
2=1 (g ,p ( 55)
(55)
Sr In Table (II) we compare the distances, d(P,p), between
the original synthetic data, (f), and the unrestored data, the
S C ,4i data restored only using the FEM method, the data restored
 using only the TVnorm minimization, and the data restored
using both techniques. For each technique, the mean distance
Fig. 15. Manifold (FEM) smoothing results for S (left), resulting ODF (right). p(d(, p)), as well as the variance o2(d(P, p)), between the
distributions in corresponding voxels is presented.
~ eg
~ \i
S As shown in Table (II), the TV restoration had superior
performance to the FEM technique, both in terms of the mean
S. error and variance of the error. The combination of techniques
Shad a lower mean error than either the FEM or TV restoration,
However the variance of the error was higher than that with
Either technique alone.
LJ .4
k
Fig. 16. Lattice (TV) smoothing results for S (left), resulting ODF (right).
From Figure (14), it can be seen that the noise has a large
influence on the smoothness of the ODF. As expected from the
The denoising algorithm was applied to a dataset consisting of
47 diffusion weighted images of a rat spinal cord. Axial slices
of one such image, before and after denoising are shown in
Figure (18).
The ringing artifacts visible near the sample boundary in the
raw DWI in Figure (18) have been noticeably decreased. Note
that the edges in the image have been well preserved. From
the expected direction images in Figure (19), it is clear that
. //


; ~s
i
. a L.. _:., L tL
L,
Fig. 18. Original diffusionweighted image (left), and denoised (right) from
spinal cord data.
the white matter fiber tracts are predominantly in the axial
direction, which is represented by the blue color.
Fig. 19. Original expected direction image (left), and denoised (right) from
spinal cord data.
Figures (20) and (21) show the HA and expected direction
visualization results respectively for various coronal slices of a
rat brain dataset. The chosen slices show the corpus callosum,
an anatomical structure known to have longrange white matter
tracts.
The HA images show a visible distinction between grey
and white matter. In addition, the distinction between fiber
directions is evident in the expected direction images. Red
corresponds to the left right direction and green corresponds
to updown and blue is inout of the page.
V. CONCLUSION
In this paper, we presented a new variational formulation for
restoring HARDI data, an FEM technique for implementing
the restoration, a novel and efficient technique for computing
the water molecule diffusion probability over a sphere of
directions, and two novel techniques for visualizing the water
molecule diffusion probability fields (computed from this
restored data). To the best of our knowledge, there are no
reported results on efficient computation and visualization
of these probability fields. Our formulation of the HARDI
restoration involves two types of smoothness constraints. The
first is smoothness over the spherical domain of acquisition
directions, and the second is smoothness between neighboring
voxels in the cartesian domain. The smoothing technique is
capable of preserving discontinuities in the data. This was
demonstrated on synthetic and real anatomical data. By using
Jdivergence as a measure of distance between distribution,
we were able to show quantitatively that the combination of
restoration techniques performs better than either technique
alone.
A statistical property of the ODF, namely entropy, was shown
to be a useful indicator of anisotropy computed from HARDI
data. Novel results showing how expected direction images
computed from the ODF could be used to visualize diffusion
direction and anisotropy in HARDI data were presented. The
effect of the restoration on these measures was shown to
improve the clarity of the images.
ACKNOWLEDGEMENT
This research was supported in part by the grant NIHNS42075
and by Siemens Corporate Research (Princeton, NJ). We wish
to thank Sara Berens and Robert Yezierski for providing the
spinal cord sample and Ron Hayes for the brain.
REFERENCES
[1] E. O. Stejskal and J. E. Tanner, "Spin diffusion measurements: Spin
echoes in the presence of a timedependent field gradient," J Chem.
Phys., vol. 42, no. 1, pp. 288292, 1965.
[2] P. J. Basser, J. Mattiello, and D. LeBihan, "MR diffusion tensor
spectroscopy and imaging," Biophys. J., vol. 66, pp. 259267, 1994.
[3] P. J. Basser, J. Mattiello, and D.LeBihan, "Estimate of the effective self
diffusion tensor from the NMR spin echo," J. Magn. Reson., vol. B, no.
103, pp. 247254, 1994.
[4] P. J. Basser, S. Pajevic, C. Pierpaoli, J. Duda, and A. Aldroubi, "In vivo
fiber tractography using DTMRI data," Magn. Reson. Med., vol. 44, pp.
625 632, 2000.
[5] P. J. Basser, "New histological and physiological stains derived from
diffusiontensor MR images," Ann. N YAcad. Sci., vol. 820, pp. 123
138, 1997.
[6] D. K. Jones, A. Simmons, S. C. R. Williams, and M. A. Horsfield, "Non
invasive assessment of axonal fiber connectivity in the human brain via
diffusion tensor MRI," Magn. Reson. Med., vol. 42, pp. 3741, 1999.
[7] N. Makris, A. J. Worth, A. G. Sorensen, G. M. Papadimitriou, O. Wu,
T. G. Reese, V. J. Wedeen, and et al, "Morphometry of in vivo human
white matter association pathways with diffusionweighted magnetic
resonance imaging," Ann. Neurol., vol. 42, pp. 951962, 1999.
[8] T. E. Conturo, N. F. Lori, T. S. Cull, E. Akbudak, A. Z. Snyder, J. S.
Shimony, R. C. McKinstry, H. Burton, and M. E. Raichle, "Tracking
neuronal fiber pathways in the living human brain," Proc. Natl. Acad.
Sci., vol. 96, pp. 1042210427, 1999.
[9] S. Mori, B. J. Crain, V. P. Chacko, and P. C. M. van Zijl, "Three
dimensional tracking of axonal projections in the brain by magnetic
resonance imaging," Ann. Neurol., vol. 45, pp. 265269, 1999.
[10] R. Xue, P. C. M. van Zijl, B. J. Crain, H. Solaiyappan, and S. Mori, "In
vivo threedimensional reconstruction of rat brain axonal projections by
diffusion tensor imaging," Magn. Reson. Med., vol. 42, pp. 11231127,
1999.
[11] B. C. Vemuri, Y. Chen, M. Rao, Z. Wang, T. McGraw, T. Mareci, S. J.
Blackband, and P. Reier, "Automatic fiber tractography from DTI and
its validation," in IEEE Intl. on Biomedical Imaging, 2002, pp.
501504.
[12] T. McGraw, B. Vemuri, Z. Wang, Y. Chen, M. Rao, and T. Mareci, "Line
integral convolution for visualization of fiber tract maps from DTI," in
Fifth Intl. Conf on MICCAI, Tokyo, Japan, 2002, pp. 61522.
[13] T. McGraw, B. C. Vemuri, Y. Chen, M. Rao, and T. Mareci, "DTMRI
denoising and neuronal fiber tracking," Medical Image Analysis, vol. 8,
no. 2, pp. 95111, June 2004.
[14] B. C. Vemuri, Y. Chen, M. Rao, T. McGraw, Z. Wang, and T. Mareci,
"Fiber tract mapping from diffusion tensor MRI," in VLSM '01: Pro
ceedings of the IEEE Workshop on Variational and Level Set Methods
(VLSM'O1). Vancouver, Canada: IEEE Computer Society, 2001, p. 81.
[15] P Callaghan, Principles of Nuclear Magnetic Resonance Microscopy.
Oxford: Clarendon Press, 1991.
[16] V J. Wedeen, T. G. Reese, D. S. Tuch, M. R. Weigel, J. G. Dou,
R. Weiskoff, and D. Chessler, "Mapping fiber orientation spectra in
cerebral white matter with fourier transform diffusion MRI," in Proc. of
the 8th Annual Meeting ofISMRM, Denver, CO, 2000, p. 82.
[17] P J. Basser, "Relationships between diffusion tensor and qspace MRI,"
Magn. Reson. Med, vol. 47, pp. 392397, 2002.
[18] D. S. Tuch, R. M. Weisskoff, J. W. Belliveau, and V. J. Wedeen, "High
angular resolution diffusion imaging of the human brain," in Proc. of
the 7th Annual Meeting ofISMRM, Philadelphia, PA, 1999, p. 321.
[19] D. S. Tuch, T. G. Reese, M. R. Wiegell, and V J. Wedeen, "Diffusion
MRI of complex neural architecture," Neuron, vol. 40, pp. 885895,
2003.
[20] L. R. Frank, "Characterization of anisotropy in high angular resolution
diffusionweighted MRI," Magn Reson Med, vol. 47, no. 6, pp. 1083
1099, 2002.
[21] K. M. Jansons and D. C. Alexander, "Persistent angular structure:
new insights from diffusion magnetic resonance imaging data," Inverse
problems, vol. 19, pp. 1031 1046, 2003.
[22] Y. Chen, W. Guo, Q. Zeng, X. Yan, F. Huang, H. Zhang, G. He, B. C.
Vemuri, and Y. Liu, "Estimation, smoothing, and characterization of
apparent diffusion coefficient profiles from high angular resolution dwi."
in CVPR (1), 2004, pp. 588593.
[23] O. Coulon, D. C. Alexander, and S. R. Arridge, "A regularization scheme
for diffusion tensor magnetic resonance images," in Proc. of the 17th
Intl. Conf on Information Processing in Medical Imaging, Davis, CA,
2001, pp. 92105.
[24] T. F. Chan and J. Shen, "Restoration of nonflat image features: models
and algorithms," SIAM J. Appl. Math., vol. 61, no. 4, pp. 13381361,
2000.
[25] D. Tschumperle and R. Deriche, "Variational frameworks for DTMRI
estimation, regularization and visualization," ICCV, pp. 116121, 2003.
[26] C. Chefd'hotel, D. Tschumperl, R. Deriche, and O. Faugeras, "Regular
izing flows for constrained matrixvalued images," Journal ofMathemat
ical Imaging and Vision, vol. 20, no. 12, pp. 147162, 2004. [Online].
Available: ftp://ftpsop.inria.fr/odyssee/Publications/2004/chefdhotel
tschumperleetal:04.pdf
[27] J. Weickert and T. Brox, "Diffusion and regularization of vector and
matrixvalued images," Contemporary Mathematics, vol. 313, pp. 251
268, 2002.
[28] F. Wang, B. C. Vemuri, M. Rao, and Y. Chen, "A new and robust
information theoretic measure and its application to image alignment,"
in Proc. of the 18th Intl. Conf on Information Processing in Medical
Imaging, 2003, pp. 388400.
[29] "Cumulative residual entropy, a new measure of information and
its application to image alignment," ICCV, pp. 548553, 2003.
[30] T. E. Conturo, R. C. McKinstry, E. Akbudak, and B. H. Robinson,
"Encoding of anisotropic diffusion with tetrahedral gradients : A general
mathematical diffusion formalism and experimental results," Magn.
Reson. Med, vol. 35, pp. 399412, 1996.
[31] C. F. Westin, S. E. Maier, B. Khidir, P Everett, F. A. Jolesz, and
R. Kikinis, "Image processing for diffusion tensor magnetic resonance
imaging," in Second Intl. Conf on MICCAI, 1999, pp. 441452.
[32] C. Poupon, C. A. Clark, V Frouin, J. Regis, I. Bloch, D. L. Bihan, and
J.F. Mangin, "Regularization of diffusionbased direction maps for the
tracking of brain white matter fascicles," Neurolmage, vol. 12, no. 2,
pp. 184195, 2000.
[33] G. J. M. Parker, C. A. M. WheelerKinghtshott, and G. J. Baker,
"Distributed anatomical brain connectivity derived from diffusion tensor
imaging," in Proc. of the 17th Intl. Conf on Information Processing in
Medical Imaging, Davis, CA, 2001, pp. 106120.
[34] P G. Batchelor, D. L. G. Hill, F. Calamante, and D. Atkinson, "Study
of connectivity in the brain using the full diffusion tensor from MRI," in
Proceedings of the 17th Intl. Conf on Information Processing in Medical
Imaging, Davis, CA, 2001, pp. 121133.
[35] M. Bj6memo, A. Brun, R. Kikinis, and C. F. Westin, "Regularized
stochastic white matter tractography using diffusion tensor MRI," in
Fifth Intl. Conf on MICCAI, vol. 1, 2002, pp. 435442.
[36] L. O'Donnell, S. Haker, and C.F. Westin, "New approaches to estima
tion of white matter connectivity in diffusion tensor MRI: Elliptic PDEs
and geodesics in a tensorwarped space," in .. Intl. Conf on MICCAI,
Tokyo, Japan, 2002, pp. 459466.
[37] T. Behrens, H. JohansenBerg, M. Woolrich, S. Smith, C. Wheeler
Kingshott, P Boulby, G. Barker, E. Sillery, K. Sheehan, O. Ciccarelli,
A. Thompson, J. Brady, and P Matthews, "Noninvasive mapping
of connections between human thalamus and cortex using diffusion
imaging," Nature Neuroscience, vol. 6, no. 7, pp. 750757, July 2003.
[38] E. Ozarslan and T. H. Mareci, "Generalized diffusion tensor imaging
and analytical relationships between diffusion tensor imaging and high
angular resolution diffusion imaging," Magn Reson Med, vol. 50, pp.
955965, 2003.
[39] P Sundgren, Q. Dong, D. G6mezHassan, S. Mukherji, P Maly, and
R. Welsh, "Diffusion tensor imaging of the brain: review of clinical
applications," Neuroradiology, vol. 46, pp. 339 350, May 2004.
[40] E. Ozarslan, B. C. Vemuri, and T. H. Mareci, "Generalized scalar
measures for diffusion MRI using trace, variance and entropy," Magn
Reson Med, vol. 53, no. 4, pp. 866 876, 2005.
[41] B. Tang, G. Sapiro, and V. Caselles, "Diffusion of general data on non
flat manifolds via harmonic maps theory: The direction diffusion case,"
IJCV, vol. 36, no. 2, pp. 149 161, 2000.
[42] L. R. Frank, "Anisotropy in high angular resolution diffusionweighted
MRI," Magn Reson Med, vol. 45, no. 6, pp. 935939, 2001.
[43] T. Mclnemey and D. Terzopoulos, "A dynamic finite element surface
model for segmentation and tracking in multidimensional medical
images with application to cardiac 4D image analysis," Comp. Med.
Imag. Graph., vol. 19, no. 1, pp. 6983, 1995. [Online]. Available:
citeseer.ist.psu.edu/mcinemey95dynamic.html
[44] G. Dhatt and G. Touzot, The finite element method displayed. New
York: Wiley, 1984.
[45] P Blomgren and T. F. Chan, "Color TV: total variation methods
for restoration of vectorvalued images," IEEE Transaction on Image
Processing, vol. 7, no. 3, pp. 304309, 1998.
[46] T. Chan and P Mulet, "On the convergence of the lagged diffusivity
fixed point method in total variation image restoration," SIAM Journal
on Numerical Analysis, vol. 36, no. 2, pp. 354367, 1999.
[47] F. Schwabl, Quantum Mechanics. Berlin: SpringerVerlag, 1989.
[48] G. Arfken and H. Weber, Mathematical Methods for Physicists. Lon
don: Academic Press, 1995.
[49] D. S. Tuch, "Qball imaging," Magn Reson Med., vol. 52, no. 6, pp.
13581372, 2004.
[50] T. M. Cover and J. A. Thomas, Elements of Information Theory. John
Wiley & sons, 1991.
[51] K. Zyczkowski, "Renyi extrapolation of shannon entropy," Open
Inf Dyn., vol. 10, pp. 297310, 2003.
Fig. 20. Original (left), and denoised (right) HA images for coronal slices Fig. 21. Original (left), and denoised (right) ED images for coronal slices
of rat brain. of rat brain.
