A novel tensor distribution model for the diffusion
weighted MR signal1
Bing Jian a, Baba C. Vemuri a, Evren Ozarslan b, Paul Carney c and
Thomas Mareci d
aDepartment of Computer and Information Science and F;l i,i;i i, ii. University of
Florida, PO.Box 116120, Gainesville, FL 32611, USA
bSection on Tissue Biophysics and Biomimetics, LIMB, NICHD, National Institutes of
Health, Bethesda, MD 20892, USA
CDepartment of Pediatrics, University of Florida, Gainesville, FL 32610, USA
dDepartment of Biochemistry and Molecular Biology, University of Florida, P.O.Box
100245, Gainesville, FL 32610, USA
CISE Technical Report 06291
Abstract
Diffusion MRI is a noninvasive imaging technique that allows the measurement of water
molecule diffusion through tissue in vivo. The directional features of water diffusion allow
one to infer the connectivity patterns prevalent in tissue and possibly track changes in this
connectivity over time for various clinical applications. In this paper, we present a novel sta
tistical model for diffusionweighted MR signal attenuation which postulates that the water
molecule diffusion can be characterized by a continuous mixture of diffusion tensors. An
interesting observation is that this continuous mixture and the MR signal attenuation are
related through the Laplace transform of a probability distribution over symmetric positive
definite matrices. We then show that when the mixing distribution is a Wishart distribution,
the resulting closed form of the Laplace transform leads to a Rigauttype a uipl ni fractal
expression, which has been phenomenologically used in the past to explain the MR signal
decay but never with a rigorous mathematical justification until now. Our model not only
includes the traditional diffusion tensor model as a special instance in the limiting case,
but also can be adjusted to describe complex tissue structure involving multiple fiber pop
ulations. Using this new model in conjunction with a spherical deconvolution approach,
we present an efficient scheme for estimating the water molecule displacement probability
functions on a voxelbyvoxel basis. Experimental results on both simulations and real data
are presented to demonstrate the robustness and accuracy the proposed algorithms.
1 This research was in part supported by NIH EB007082, NIH NS42075 to BCV and NIH
EB004752 to PC & TM.
Preprint submitted to Elsevier Science
17 December 2006
1 Introduction
Diffusionweighted imaging (DWI) is a magnetic resonance (MR) technique ex
ploiting the sensitivity of the MR signal to the Brownian motion of water molecules.
It adds to the conventional relaxationweighted magnetic resonance imaging (MRI)
the capability of measuring the water diffusion characteristics in local tissue, which
may be substantially altered by diseases, neurologic disorders, and during neurode
velopment and aging. DWI has steadily evolved into an important clinical tool since
its sensitivity to the evaluation of early ischemic stages of the brain was shown
(Moseley et al., 1990b). The directional dependence of water diffusion in fibrous
tissues, like muscle and whitematter in the brain, provides an indirect but power
ful means to probe the anisotropic microstructure of these tissues (Cleveland et al.,
1976; LeBihan et al., 1986; Moseley et al., 1990a). As of today, DWI is the unique
noninvasive technique capable of quantifying the anisotropic diffusion of water
molecules in tissues allowing one to draw inference about neuronal connections
between different regions of the central nervous system.
In diffusion MRI, most applications rely on the fundamental relationship between
the MR signal measurement S(q) and the average particle displacement probability
density function (PDF) P(r) which is given by the following Fourier transform
(Callaghan, 1991):
S(q) =So J P(r) eiqrdr (1)
where So is the signal in the absence of any diffusion weighting gradient, r is the
displacement vector and q = 7'Gg, where 7 is the gyromagnetic ratio, 6 is the
diffusion gradient duration, and G and g are the magnitude and direction of the
diffusion sensitizing gradients respectively.
A widely accepted method for diffusionweighted imaging studies is to use appar
ent difl,'itii coefficient (ADC) profiles which is governed by the StejskalTanner
equation (Stejskal and Tanner, 1965):
S(q) = So exp(bDapp) (2)
where b is the diffusion weighting factor depending on the strength as well as the
effective time of diffusion and Dapp is the so called apparent diffhit i coefficient.
Diffusion tensor MRI (DTMRI or DTI), introduced by Basser et al. (1994), pro
vides a relatively simple way of quantifying diffusional anisotropy as well as pre
dicting the local fiber direction within the tissue from multidirectional diffusion
weighted MRI data. DTI assumes a displacement probability characterized by an
oriented Gaussian probability distribution function. Substituting the oriented Gaus
sian model for P(r) in Eq.(1) leads to a signal decay given by
S(q) So exp (tqDq) So exp (bgDg) (3)
where b = llql12t is the bfactor, t is the effective diffusion time and D is the
apparent dif/,'it ini tensor.
Despite its modest requirements, the DTI model has been shown to be quite suc
cessful in regions of the brain and spinal cord with significant whitematter coher
ence and has enabled the mapping of anatomical connections in the central nervous
system (Conturo et al., 1999; Mori et al., 1999; Basser et al., 2000). However, the
major drawback of diffusion tensor MRI is that it can only reveal a single fiber ori
entation in each voxel and fails in voxels with orientational heterogeneity (IVOH)
(von dem Hagen and Henkelman, 2002; Tuch et al., 2002), which makes DTI an
inappropriate model for use in the presence of multiple fibers within a voxel.
This limitation of DTMRI has prompted interest in the development of both im
proved image acquisition strategies and more sophisticated reconstruction meth
ods. By sampling the diffusion signal on a threedimensional Cartesian lattice, the
qspace imaging (QSI) technique, also referred to as diffusion spectrum imaging
(DSI) (Wedeen et al., 2000), measures the microscopic diffusion function directly
based on the Fourier relation between the diffusion signal and the diffusion func
tion, without recourse to a model of the diffusion process. However, the sampling
burden of QSI makes the acquisition timeintensive and limits the wide spread ap
plication of QSI. Tuch et al. (1999) developed a clinically feasible approach called
high angular resolution diffusion imaging (HARDI), in which apparent diffusion
coefficients are measured along many directions distributed almost isotropically on
a spherical shell in the diffusion wave vector space.
Various models and numerical fitting procedures have been proposed to relate the
diffusion weighted MR signal to the underlying diffusivity function. The first class
of approaches studies the diffusivity profile assuming an exponential attenuation
based on the StejskalTanner expression (Stejskal and Tanner, 1965). It has been
shown that the diffusivity function has a complicated structure in voxels with ori
entational heterogeneity (von dem Hagen and Henkelman, 2002; Tuch et al., 2002).
Several studies proposed to represent the diffusivity function using the spherical
harmonic expansion (Frank, 2002; Alexander et al., 2002). Ozarslan and Mareci
(2003) showed that the ADC profiles can be equivalently modeled using Cartesian
tensors of rank higher than 2, giving rise to a generalization of DTI's rank2 tensor
model.
Since the peaks of the diffusivity profile do not necessarily yield the orientation
of the distinct fiber populations, a second class of approaches attempt to transform
the multidirectional signals into a probability function describing the probability of
water molecular diffusion over a sphere of directions. Recently, Tuch et al. (2003)
proposed the so called qball imaging (QBI) method, in which the radial integral of
the displacement PDF is approximated by the spherical FunkRadon transform. The
end result obtained by the radial projection of the threedimensional displacement
PDF via a line integral is a convolution of the real probability values with a 0
order Bessel function (Tuch, 2004) and not the probability values themselves. This
convolution gives rise to an undesirable "contamination" of the probability along
one direction with probabilities from other directions. More recent studies have
expressed QBI's FunkRadon transform in a spherical harmonic basis (Anderson,
2005; Hess et al., 2006; Descoteaux et al., 2006).
Ozarslan et al. (2004) generated signal values from the threedimensional wave
vector by fitting higherorder tensor models to the measurements from a spherical
acquisition scheme, and then evaluated the displacement probabilities by employ
ing an FFT. Another reconstruction algorithm referred to as PASMRI was pro
posed by Jansons and Alexander (2003). This method computed the socalled the
persistent angular structure (PAS), which is a function on a fixed spherical shell in
threedimensions, by assuming its Fourier transform best fits the measurements and
incorporating a maximumentropy condition. Although the results are promising,
both of these schemes are computationally expensive. More recently, a robust and
fast method, called the diffusion orientation transform (DOT), was introduced in
(Ozarslan et al., 2006b). By expressing the Fourier relation in spherical coordinates
and evaluating the radial part of the integral analytically, DOT is able to transform
the diffusivity profiles into probability profiles either directly or parametrically in
terms of a spherical harmonic series. In this case, the estimated probability func
tion is also "impure" in the sense that the end result is the true probability values
convolved with a function which cannot be specified analytically. However, much
of the blurring in the DOT is due to the monoexponential decay assumption of
the MR signal, and hence can be avoided by using the extension of the transform
to multiexponential attenuation as described in (Ozarslan et al., 2006b). However,
this extension would necessitate collecting data on several spherical shells in the
wave vector space.
There exists a third class of methods in which some multicompartmental models or
multiplefiber population models have been used to model the diffusionattenuated
MR signal. Tuch et al. (2002) proposed to model the diffusion signal using a mix
ture of Gaussian densities:
S(q) = So w exp (b gTDg) (4)
J
where wj is the apparent volume fraction of the compartment with diffusion ten
sor Dj. Two problems accompany the use of multicompartment models. First, the
number of such compartments has to be prespecified, presenting a modelselection
problem. Second, the nonlinear fitting procedure is unstable and heavily depends
on the choice of the starting point because of the local minima in the objective
function. Other similar modelbased approaches which extend to the multiplefiber
case by including multiple restricted populations have been proposed in literature
(Behrens et al., 2003; Assaf et al., 2004). However, both methods have the same
modelselection and fitting problems as the mixture models discussed above.
To address the modelselection problem, Tournier et al. (2004) employed the spher
ical deconvolution method, assuming a distribution, rather than a discrete number,
of fiber orientations. Under this assumption, the diffusion MR signal is the convo
lution of the fiber orientation distribution (FOD), which is a realvalued function
on the unit sphere, with some kernel function representing the response derived
from a single fiber. A number of spherical deconvolution based approaches have
been proposed (Anderson, 2005; Alexander, 2005; Tournier et al., 2006) with dif
ferent choices of FOD parameterizations, deconvolution kernels and regularization
schemes.
In jhi'% paper, we present a novel probabilistic model based on the traditional dif
fusion tensor model (Basser et al., 1994). First, we assume that each voxel is asso
ciated with an underlying probability distribution defined on the space of diffi ,'itl
tensors (the manifold of 3 x 3 symmetric pt",iiii~ d filniiil matrices). Conceptually,
our model can be viewed as a natural extension of the multiplecompartment mod
els (Tuch et al., 2002; Assaf et al., 2004) since both previous singletensor and
multiplecompartment models can be treated as special cases of our generalized
model when the underlying distribution is represented by a Dirac delta function
or a linear combination of Dirac delta functions. From hii' point of view, the con
ventional tensor fitting can be interpreted as the location parameter 2 estimation
problem solved by maximizing the likelihood function. Our continuous probabilis
tic structure naturally leads to the Bayesian formulation of the location parameter
estimation. Moreover, by extending the discrete mixture to the continuous mixture,
we make an interesting observation that the resulting continuous mixture model and
MR signal attenuation are related ;r,ln',,og a Laplace tinlrm''iii defined for matrix
variate functions. It is worth noting that the Laplace tir tril'i'i can be evaluated
in closed form for the case when the mixing distribution is a Wishart distribution.
We show that the resulting closed form leads to a Rigauttype function which has
been used in the past to explain the MR signal decay (Kipf et al., 1998) but never
with a rigorous mathematical juiifi, Olitil until now. This theoretical result depicts
surprising consistency with the experimental observations.
The rest of the paper is organized as follows: Section 2 presents the technical de
tails of our new model. In Section 3.1, we introduce a new diffusion tensor imaging
framework as a direct application of the proposed model and compare it with the
traditional diffusion tensor imaging method. A method for resolving multiple fiber
orientations based on the proposed model is also developed in Section 3.2. Section
4 contains results of application of our new methods to synthetic and real diffu
sion weighted MRI data and comparisons with results from representative existing
methods. We draw conclusions in Section 5.
2 Here the domain of the location parameter is the manifold of 3 x 3 symmetric positive
matrices.
2 Theory
As mentioned in the Section 1, we assume that each voxel is associated with an
underlying probability distribution defined on the space of diffusion tensors. For
mally speaking, by assumption, at each voxel there is an underlying probability
measure induced on the manifold of n x n symmetric positivedefinite matrices,
denoted by P,. 3 Let F be the underlying probability measure, then we can model
the diffusion signal by:
S(q) So exp[bgTDg] dF So f (D) exp[bgTDg] dD (5)
where f(D) is the density function of F with respect to some carrier measure dD
on P,. Note that Eq.(5) implies a more general form of mixture model with f(D)
being a mixing density over the covariance matrices of Gaussian distributions.
Clearly, our model simplifies to the diffusion tensor model when the underlying
probability measure is the Dirac measure.
Since bgTDg in Eq.(5) can be replaced by trace(BD) where B bggT, the diffu
sion signal model presented in the form of (5) can be expressed as follows
S(q) = So exp (bgTDg) dF = So exp (trace(BD)) dF, (6)
which is exactly the Laplace transform of the probability measure F on P, (For
definition of Laplace transform on P, see Appendix A.1).
This expression naturally leads to an inverse problem: recovering of a distribution
F defined on P, that best explains the observed diffusion signal S(q). This is an
illposed problem and in general is intractable without prior knowledge on the prob
abilistic structure. In conventional DTI, the diffusion tensor is usually estimated by
solving a linear or nonlinear least squares problem, which amounts to applying
the maximum likelihood estimator. Since our approach allows us to model the un
certainty of the diffusion tensor estimation, one can relate the proposed model to
a Bayesian inference problem which views the diffusion tensor as random vari
able(matrix) having some known prior distribution. Note that in DTI the diffusion
tensor can be interpreted as the concentration matrix (inverse of the covariance ma
trix) of the Gaussian distribution in the qspace. It is a common practice to put a
Wishart distribution, on the concentration matrix in multivariate analysis. In fact,
the Wishart distribution is the standard conjugate prior for the concentration matrix
estimation problem. Based on this motivation, in this paper, we propose to model
the underlying distribution through some parametric probability family on P,, in
particular, the Wishart distribution or the mixture of Wishart distributions.
3 Throughout this paper, P, is by default the manifold of 3 x 3 symmetric positivedefinite
matrices.
In the following, we first briefly introduce the definition of Wishart distribution as
well as its relevant properties, then we analytically derive that when the mixing
distribution in the proposed continuous mixture tensor model is a Wishart distri
bution, the Laplace transform leads to a Rigauttype asymptotic fractal law for the
MR signal decay which has been phenomenologically used previously to explain
the MR signal decay but never with a rigorous mathematical justification until now.
As one of the most important probability distribution families for nonnegative
definite matrixvalued random variables ("random matrices"), the Wishart distribu
tion (Wishart, 1928) is most typically used when describing the covariance matrix
of multinormal samples in multivariate statistics (Murihead, 1982).
Usually, the probability density function of Wishart distribution with respect to the
Lebesgue measure dY is defined as follows (Murihead, 1982):
Definition 1 A random matrix YE cP, is said to have the (central) Wishart distri
bution W,(p, E) with scale matrix E and p degrees of freedom, n < p, if the joint
distribution of the entries of Y has the density function:
f(Y) c Y 1)/2 / exp( trace(E Y) (7)
with E cE P, and c = 2p/2F,(p/2)1 where F, is the multivariate gamma func
tion:
F,(p) = exp (trace(Y)) ylp(+)/2 dY
and I denotes the determinant of a matrix.
Recently Letac and Massam (1998) showed Wishart distribution can be viewed as
a natural generalization of the gamma distribution by introducing the definition in
Eq.(8)
Definition 2 (Letac and Massam, 1998)
For andfor p in A 1,..., U ) the Wishart distribu
ForY c PT,, andforp in A 21 .. U 1 1
tion ,7,y with scale parameter E and shape parameter p is defined as
d7p,s(Y) = F,(p)1 lYl(n+1)/2 p exp(trace(EZY)) dY, (8)
where T, is the multivariate gamma function:
r,(p) J exp (trace(Y)) yip(n+1)/2 dY
and  denotes the determinant of a matrix.
Note that the above definition differs slightly from the traditional notation W,(p, E)
for Wishart distribution (e.g., in Anderson (1958); Murihead (1982)) and the corre
spondence between the two notations is simply given by 7p/2,2 = W, (p, E). In the
rest of this paper, we will use the notation q7,y as provided in Letac and Massam
(1998).
Clearly, the definition in Eq.(8) leads to a natural generalization of the gamma
distribution. It can be further shown that the Wishart distribution preserves the fol
lowing two important properties of the gamma distribution:
Theorem 1 The expected value of a random variable(matrix) following distribu
tion 7p,s is pE.
Theorem 2 The Laplace tinrcu'i',i, of the (generalized) gamma distribution 7,Y
is
/exp(trace(eu)) 7p,s(du) 1, + OEP' where ( + 1) Pz. (9)
Substituting the general probability measure F in (6) by the Wishart measure 7p,e
and noting that B bg gT, we have
S(q)/So IF, + BE' = (1 + (bgTEg)) (10)
Consider the family of Wishart distributions 7,y, with fixed expected value D pE.
In this case, the above expression takes the form:
S(q) So (1 + (bgTDg)/p) (11)
This is a familiar Rigauttype asymptotic fractal expression (Rigaut, 1984) when
the argument is taken to be the ADC associated with the expected tensor of the
Wishart distributed diffusion tensors4 The important point is that this expression
implies a signal decay characterized by a powerlaw in the large q hence large
b region exhibiting asymptotic behavior. This is the expected asymptotic behavior
for the MR signal attenuation in porous media (Sen et al., 1995). Note that al
though this form of a signal attenuation curve had been phenomenologically fitted
to the diffusionweighted MR data before (K6pf et al., 1998), until now, there was
no rigorous derivation of the Rigauttype expression used to explain the MR sig
nal behavior as a function of bvalue. Therefore, this derivation may be useful in
understanding the apparent fractallike behavior of the neural tissue in diffusion
weighted MR experiments (Kipf et al., 1998; Ozarslan et al., 2006a). Note that in
Eq.(11) the value of p depends on the dimension of the space in which diffusion
is taking place. Although for fractal spaces this exponent can be a noninteger, the
4 Note that the form of 11 is slightly different from Rigaut (1984)'s own formula; how
ever, it possesses the desirable properties of the original formula such as concavity and the
a,\ mpl nlic linearity in the loglog plotshence the phrase "Rigauttype".
analog of DebyePorod law of diffraction (Sen et al., 1995) ensures that in three
dimensional space the signal should have the asymptotic behavior, S(q) ~ q4
Since b oc q2 a reasonable choice for p is 2.
1 .0000 .
0 000i ',. .. .. "' ""*.' ...: rigaul.ype
0,1000 i '. t a urnptli i fracal
Gaussian '. . ,
all 10000 tenors ', : '
0.0001 L.j
0.01. 0.10 1.00
(7 (,1" 1)
Fig. 1. The Wishart distributed tensors lead to a Rigauttype signal decay
To empirically validate (11), the following experiment is designed. We first draw a
sequence of random samples of increasing sample size from a Wishart distribution
with p = 2, and then for each random sample {D1,..., D,} of rank2 tensors,
the corresponding multiexponential signal decay can be simulated using a discrete
mixture of tensors as follows:
E(q) (q)/So exp ( bgTDg). (12)
To illustrate the relation between signal decay behavior and the sample size, we plot
the signal decay curves for different sample sizes in Fig. 1, by fixing the direction
of diffusion gradient q and increasing the strength q = q The left extreme dotted
curve depicts the signal decay from a monoexponential model where the diffu
sion tensor is taken to be the expected value of the Wishart distribution. The right
extreme solid curve is the Rigauttype decay derived from (11). Note that the tail of
the solid curve is linear indicating the powerlaw behavior. The dotted curves be
tween these two extremes exhibit the decay for random samples of increasing size
but smaller than 10,000. The dashed curve uses a random sample of size 10,000 and
is almost identical to the expected Rigauttype function. As shown in the figure, a
single tensor gives a Gaussian decay, and the sum of a few Gaussians also produces
a curve whose tail is Gaussianlike, but as the number of tensors increases, the at
tenuation curve converges to a Rigauttype asymptotic fractal curve with desired
linear tail and the expected slope in the double logarithmic plot.
It is well known that the gamma distribution 7p, with integer p is also the distri
bution of the sum of p independent random variables following exponential dis
tribution with parameter a. It follows from the central limit theorem that if p (not
necessary integer) is large, the gamma distribution 7,, can be approximated by the
normal distribution with mean pa and variance pa2. More precisely, the gamma
distribution converges to a normal distribution when p goes to infinity. The similar
behavior applies to the Wishart distribution. Note when p tends to infinity, we have
S(q) So exp(bgTDg) (13)
which implies that the monoexponential model can be viewed as a limiting case
(p + o) of our model. Hence Eq.(11) can be seen as a generalization of Eq.(3).
By the linearity of the Laplace transform, the biexponential and multiexponential
models can be derived from the Laplace transform of the discrete mixture of Wishart
distributions as well.
3 Applications
3.1 A new framework for difli ii ,n tensor estimation
The model in Eq.(6) also suggests a new method for the estimation of diffusion
tensors from diffusionweighted images. We first consider a set of diffusion mea
surements performed at a voxel containing a single fiber bundle. In this case, it is
natural to use the Wishart distribution 7y, as the mixing distribution in Eq.(6) and
thus the following equation is obtained:
( ) trace(BZ) 1
So "
or in the matrix form:
(S1) Bxzz 2B2 (So)P 1
(S2) B ** *.. 2B,z (XX 1
(14)
(SK) p Bxx, 2B, \ \ 1
where K is the number of measurements at each voxel and Bij and Zj are the six
components of the matrices B and Y, respectively. Note that in the above expression
the components of the matrices B and a should be ordered consistently.
Hence the fitting problem can be reformulated as a linear system problem. As a re
sult, the So and the six components in a can be estimated efficiently by using linear
regression as in the traditional DTI method. Since this new framework of diffu
sion tensor estimation is not the focus of this paper, we simply formulate a linear
system which shows the analogy to tensor estimation using the Wishart model and
the traditional DTI model respectively. Many complicated methods which involve
nonlinear optimization and enforce the positivity constraint on the diffusion tensor,
for example, (Chefd'Hotel et al., 2002; Wang et al., 2004), can be applied to the
Wishart model proposed here. Similarly, the resulting diffusion tensor field repre
sented by the estimated p a at each voxel and can be then analyzed by numerous
existing diffusion tensor image analysis methods (Weickert and H.Hagen, 2005).
3.2 Multifiber reconstruction using deconvolution
However, the single Wishart model can not resolve the intravoxel orientational het
erogeneity due to the single diffusion maximum in the model. Actually, the Laplace
transform relation between the MR signal and the probability distributions on P,
naturally leads to an inverse problem: to recover a distribution on P, that best ex
plains the observed diffusion signal. In order to make the problem tractable, several
simplifying assumptions are made as follows.
We first propose a discrete mixture of Wishart distribution model where the mix
ing distribution in Eq.(6) is expressed as dF = Zi dE, In this model
(pi, Ei) are treated as basis and will be fixed as described below. This leaves us
with the weights w, as the unknowns to be estimated. Note the number of compo
nents in mixture, N, only depends on the resolution of the manifold discretization
and should not be interpreted as the expected number of fiber bundles. We assume
that all the pi take the same value p = 2 based on the analogy between the Eq.(11)
and DebyePorod law of diffraction (Sen et al., 1995) in threedimensional space
as discussed in Section 2. Since the fibers have an approximately cylindrical geom
etry, it is reasonable to assume that the two smaller eigenvalues of diffusion ten
sors are equal. In practice, we fix the eigenvalues of D = pE, to specified values
(Al, A2, A3) (1.5, 0.4, 0.4)p2/ms consistent with the values commonly observed
in the whitematter tracts (Tuch et al., 2002). This rotational symmetry leads to a
tessellation where N unit vectors evenly distributed on the unit sphere are chosen
as the principal directions of Es. In this way, the distribution can be estimated us
ing a spherical deconvolution scheme (Toumier et al., 2004). For K measurements
with qj, the signal model equation:
N
S(q) So (1 + trace(Bij))p (15)
i= 1
leads to a linear system Aw = s, where s = (S(q)/So) is the vector of normalized
measurements, w = (,,), is the vector of weights to be estimated and A is the
matrix with Aj = (1 + trace(BZji))p. It is worth noting that if we take p = oo,
the deconvolution kernel becomes the Gaussian function and the resulting w resem
bles very closely the fiber orientation estimated using continuous axially symmetric
tensors (FORECAST) method proposed in Anderson (2005).
To avoid an underdetermined system, K > N is required without an interpolation
on the measurements or exploring the sparsity constraints on w. Since the matrix A
only depends on the sampling scheme and therefore needs only onetime compu
tation, the computational burden of this method is light and comparable to that of
the traditional linear least squares estimation of diffusion tensors from the Stejskal
Tanner equation. However, the induced inverse problem can be illconditioned due
to the possible singular configurations of the linear system. In practice, the damped
least squares (DLS) (Wampler, 1986) inverse is employed to overcome the insta
bility problem. Instead of inverting small singular values, the damped least squares
technique builds a smooth function converging to zero when the singular value
tends towards zero. Note that the DLS is the closed form solution to one special
case of the Tikhonov regularization used in Toumier et al. (2006).
Ozarslan et al. (2006b) have shown that the distinct fiber orientations can be esti
mated by computing the peaks of the probability profiles, i.e. the probabilities for
water molecules to move a fixed distance along different directions. Using a simi
lar idea, we now present a way to compute the displacement probabilities using the
proposed continuous distribution of tensors model.
First, the displacement probabilities can be approximated by the Fourier transform
P(r) = E(q) exp(iq r) dq where E(q) = S(q)/So is the MR signal atten
uation. Assuming a continuous diffusion tensor model as in Eq.(5) with mixing
distribution F(D) = .I" .: ,we have
P(r)= J exp(qTDqt) dF(D) exp(iq r) dq
N (16)
D exp(rTD, r/4t)
ii (47wt)3 D
where D, pE, are the expected values of :
After the displacement probability profile is obtained as as a realvalued function
on the sphere can be computed, one could simply follow the common strategy
used in (Tournier et al., 2004; Ozarslan et al., 2006b) to estimate the distinct fiber
orientations by first fitting the probability profile using spherical harmonics and
then finding the peaks of the resulting spherical function. Note here the number of
fiber populations can be obtained as the byproduct of this procedure. An alternative
approach is applying statistical hypothesis testing like the Ftest used in (Alexander
et al., 2002) to determine the appropriate number of fiber populations.
4 Experiments
4.1 Difl,\ii1 tensor field estimation using Wishart model
We first test the method proposed in section 3.1 on a synthetic data set represent
ing singlefiber diffusion with sinusoidally varying orientations as shown in Figure
2. The simulations employed the exact form of the MR signal attenuation from
particles diffusing inside cylindrical boundaries (Siderman and J6nsson, 1995).
Fig. 2. Diffusion tensor fitting of a simulated data set. Left: the tensor field obtained from
fitting the linearized StejskalTanner equation. Right: the tensor field using the equation in
section 2.2 with p = 2.
To quantitatively assess the performance of the proposed model, we compare the
accuracy of the tensor fits with those computed from the monoexponential Stejskal
Tanner equation by quantifying the angular deviation of the dominant eigenvectors.
Let 0 be the angle (in degrees) between the dominant eigenvector of the estimated
diffusion tensor field and the dominant eigenvector of the original tensor field. Ta
ble 1 shows the mean (p ) and standard deviation(o0) of 0 using different meth
ods for the synthetic data with different levels of Riciandistributed noise into the
magnitude of the MR signal approximated by additive Gaussian noise on real and
imaginary part. From this table, we can see that the deviation angles obtained from
our model are lower than those obtained using monoexponential StejskalTanner
equation under all noise levels.
DTI model Our model
SNR mean std. dev. mean std. dev.
No noise 11.25 7.29 11.25 7.08
25dB 11.70 7.63 11.60 7.52
20dB 14.44 8.27 14.00 7.85
15dB 15.00 8.92 14.62 8.42
Table 1
Comparison of the accuracy of the estimated dominant eigen vectors using different meth
ods under different noise levels.
4.2 Multi fiber reconstruction on simulated data
We have applied the scheme described in section 3.2 to the simulations of single
fiber and crossing fiber systems. As in previous experiments, the simulations were
performed following the exact form of the MR signal attenuation from particles
diffusing inside cylindrical boundaries (Siderman and Jinsson, 1995).
Fig. 3 shows the probability surfaces for a simulated image of two crossing fiber
bundles computed using the proposed method. The surfaces are consistent with
the underlying known fibrous structure. The curved and linear fiber bundles were
chosen so that a distribution of crossing angles is achieved across the region with
orientational heterogeneity. We notice that distinct fiber orientations are better re
solved when the different fiber bundles make larger angles with each other.
\\\\\\ / /
\\\\ \\\y / /// / / ///\ \\\ ///////////!////
\\\\\\\. /t ///////:/' ///////// \ \\\vyvs / /////////
\\\ \\v\ .\ '///'///////// \\\xxvi'. ///////// ///////
SXXX~ xxxyt~ X ^///////
(a) (b)
Fig. 3. Probability maps from a simulated image of two crossing fiber bundles computed
using (a) DOT (b) proposed method.
We also empirically investigated the performance of our reconstruction method.
Of special interest is its accuracy in fiber orientation detection in the presence of
noise. To study this problem, we first took the HARDI simulations of 1,2 and
3fiber profiles with known fiber orientations and computed the probability profiles
which are depicted in Fig.4.
In the case of the noiseless signal, the proposed method as well as QBI and DOT
are all able to recover the fiber orientation quite accurately. The Qball orientation
distribution functions (ODF) is computed by using the spherical harmonic expan
sion formula given in Anderson (2005, Eq.(21)). Since our method and the non
parametric DOT method compute the probability values directly, we fit the result
ing probability profiles from proposed method and DOT using spherical harmonics
basis for the purpose of better surface display via rendering. The existence of ana
lytical angular derivatives of spherical harmonic functions also enables us to apply
fast gradientbased numerical optimization routines to find the peaks of the proba
bility surfaces.
To provide a more quantitative assessment of the proposed method and its sensitiv
ity to noise, we performed a series of experiments as follows:
Fig. 4. Simulations of 1, 2 and 3fibers (b = 1500s/mm2). Orientations: azimuthal angles
01 = 30, 2 = {20,100}, 03 = {20, 75,135}; polar angles were all 90. Top: Qball ODF
surfaces; Bottom: Probability surfaces computed using proposed method.
* Let Sj,cean be the simulated signal magnitude along diffusion gradient gj. The
noise enters the system via Sj,noisy = Sj,cean + eRe + iIm where both cRe
and Cem are independently drawn from a normal distribution N(0, o). Note that
this scheme is commonly used to introduce a Riciandistributed noise into the
magnitude of the MR signal.
* For all 1, 2 and 3fiber systems as shown in Fig.4, the noise profile was created
as described above with increasing noise levels (a = .02, .04, .06, .08).
* For each corrupted fiber system, we recomputed the probability profiles by using
the proposed method and DOT. Similarly, the Qball ODF were computed using
the formula given by Anderson (2005, Eq.(21)). All the resulting surfaces were
represented by spherical harmonics coefficients.
* Then we estimated the fiber orientations of each system by numerically finding
the maxima of the probability surfaces with a QuasiNewton numerical optimiza
tion algorithm and computed the deviation angles b between the estimated and
the true fiber orientations.
* For each system and each noise level, the above steps were repeated 100 times to
provide a distribution of deviation angles. Table 2 reports the mean and standard
deviation of these distributions in degrees.
As expected, the deviation angles between the recovered and the true fiber orienta
tions increase with increasing noise levels and it is more challenging to accurately
resolve the distinct orientations when there are more fiber orientations. The statis
tics reported in Table 2 also indicate that the proposed method has stronger resis
tance to the noise than the DOT and the QBI methods respectively. Fig.5 and Fig.6
also illustrate this finding.
From proposed method
((a = 0) 0(a = .02) 0(a = .04) 0(a = .06) 0(a = .08)
1 fiber { 0.243} 0.65 0.39 1.19 0.65 1.66 0.87 2.19 1.27
2 fibers .74} 1.18 0.66 2.55 1.29 3.85 2.12 4.91 3.26
2 fibers
{0.69} 1.30 0.66 2.76 1.34 3.63 1.91 5.11 2.65
{1.02} 4.87 3.23 8.59 5.82 11.79 6.86 13.84 8.73
3 fibers {0.97} 5.81 3.61 7.70 5.02 11.27 6.36 12.54 7.48
{1.72} 4.92 3.32 7.94 4.59 12.57 7.09 14.27 7.66
From DOT
(o = 0) (a = .02) (a = .04) (a = .06) (a = .08)
1 fiber {0.414} 0.71 0.35 1.08 0.58 1.84 0.88 2.20 1.28
2 fibers {1.55} 1.97 0.96 3.37 1.90 5.39 2.99 7.00 4.25
2 fibers
{1.10} 1.73 1.00 3.28 1.87 4.78 2.37 6.29 3.19
{4.11} 7.89 5.71 10.82 6.66 14.56 8.74 16.68 10.21
3 fibers {3.46} 6.94 3.70 11.28 5.98 16.92 10.36 17.02 10.95
{1.68} 6.76 5.21 10.90 5.63 14.08 9.05 13.99 9.74
From QBI
(o = 0) (a = .02) (a = .04) (a = .06) (a = .08)
1 fiber {0.089} 1.28 0.75 3.34 1.97 5.94 3.19 7.67 4.16
2 fibers .45} 2.39 1.26 4.82 2.44 7.95 4.45 8.91 4.64
2 fibers
{0.42} 2.30 1.10 4.94 2.15 7.49 3.88 9.34 4.45
{0.90} 10.80 5.59 12.15 4.42 20.21 11.10 18.78 11.39
3 fibers {0.90} 11.59 5.44 13.07 4.74 19.54 11.80 20.79 10.81
{0.19} 11.66 5.18 12.25 4.93 20.36 11.50 19.10 10.18
Table 2
Mean and standard deviation values for the deviation angles < between the computed and
true fiber orientations after adding Rician noise of increasing noise levels a. Note, in all
cases, we discarded very large deviation angles that are greater than 300 when a = .02, 400
when a = .04, 500 when a = .06 or a = .08 since these large errors are mostly due to the
failure of the numerical optimization routine.
4.3 Multifiber reconstruction on real data
The rat optic chiasm is an excellent experimental validation of our approach due
to its distinct myelinated structure with both parallel and descussating optic nerve
fibers. HARDI data from optic chiasm region of excised, perfusionfixed rat ner
J1~~.~ l,,s 4k4'~i k4
4^ <'M^4r
Si rA ^
(a) ODF from QBI
kk)V _1'kA.t
4( 41 A., A
4 i k41 +k8~
(b) Proposed method
Fig. 5. Resistance to noise (2fibers, a = 0.08)
^**^'K **4*J*
^ (a) DF ^from QIt
^^ %^^^~~b ~ ~
3 k 8sr^ ,p3i^ ^
^ ^ l^^i^9t ~ ~~
(a) ODF from QBI
's dW~ P' ~ J 3 t4
(b) Proposed method
(b)Prpoedmeho
Fig. 6. Resistance to noise (3fibers, o = 0.04)
vous tissue was acquired at 14. IT using Bruker Avance imaging systems. A diffusion
weighted spin echo pulse sequence was used. Diffusionweighted images were
acquired along 46 directions with a bvalue of 1250s/mm2 along with a single
image acquired at b Os/mm2. Echo time and repetition time were 23ms and
0.5s respectively; A value and 6 value were 12.4ms and 1.2ms; bandwidth was
set to 35kHz; signal average was 10; matrix size was 128 x 128 x 5 and resolu
tion was 33.6 x 33.6 x 200/m3. The optic chiasm images were signal averaged to
67.2 x 67.2 x 200pm3 resolution prior to probability calculations.
Figure 7 shows the displacement probabilities computed from the optic chiasm
image. For the sake of clarity, we excluded every other pixel and overlaid the prob
ability surfaces on generalized anisotropy (GA) maps (Ozarslan et al., 2005). As
evident from this figure, our method is able to demonstrate the distinct fiber orien
tations in the central region of the optic chiasm where ipilateral myelinated axons
from the two optic nerves cross and form the contralateral optic tracts.
"r'
Y.I/
/t/l/ //////
/j / / / optic nerve
^J //.i /1' f4J I
., f~ ~ ~I~ Y If Yf j / optic chiasm
optic tract
Fig. 7. Probability maps computed from a rat optic chiasm data set overlaid on axially
oriented GA maps. The decussations of myelinated axons from the two optic nerves at the
center of the optic chiasm are readily apparent. Decussating fibers carry information from
the temporal visual fields to the geniculate body. Upper left corner shows the corresponding
reference (SO) image.
To investigate the capability of diffusion weighted imaging in revealing the effects
in local tissue caused by diseases or neurologic disorders, further experiments were
carried out on two data sets collected from a pair of epileptic/normal rat brains.
Under deep anesthesia, a Sprague Dawley rat was transcardial exsanguinated and
then perfused with a fixative solution of 4% paraformaldehyde in phosphate buffered
saline (PBS). The corpse was stored in a refrigerator overnight then the brain was
extracted and stored in the fixative solution. For MR measurements, the brain was
removed from the fixative solution then soaked in PBS, without fixative, for ap
proximately 12 hours (overnight). Prior to MR imaging, the brain was removed
from the saline solution and placed in a 20 mm tube with fluorinated oil (Fluorinert
FC43, 3M Corp., St. Paul, MN) and held in place with plugs. Extra care was taken
to remove any air bubbles in the sample preparation.
The multipleslice diffusion weighted image data were measured at 750 MHz us
ing a 17.6 Tesla, 89 mm bore magnet with Bruker Avance console (Bruker NMR
Instruments, Billerica, MA). A spinecho, pulsedfieldgradient sequence was used
for data acquisition with a repetition time of 1400 ms and an echo time of 28 ms.
The diffusion weighted gradient pulses were 1.5 ms long and separated by 17.5 ms.
A total of 32 slices, with a thickness of 0.3 mm, were measured with an orientation
parallel to the longaxis of the brain (slices progressed in the dorsalventral direc
tion). These slices have a fieldofview 30 mm x 15 mm in a matrix of 200 x 100.
The diffusion weighted images were interpolated to a matrix of 400 x 200 for each
slice. Each image was measured with 2 diffusion weightings: 100 and 1250s//mm2.
Diffusionweighted images with 100s/mm2 were measured in 6 gradient direc
tions determined by the vertices of an icosahedron in one of the hemispheres. The
images with a diffusionweighting of 1250s/mm2 were measured in 46 gradient
directions, which are determined by the tessellations of the icosahedron on the same
hemisphere. The 100s/mm2 images were acquired with 20 signal averages and the
1250s/mm2 images were acquired with 5 signal averages in a total measurement
time of approximately 14 hours.
Figure 9 shows the displacement probabilities calculated from excised coronal rat
brain MRI data in a, (a) control and (b) an epileptic rat. (See Fig.8 for the cor
responding SO images.) The hippocampus and entorhinal cortex is expanded and
depicts the orientations of the highly anisotropic and coherent fibers. Note voxels
with crossing orientations located in the dentate gyms (dg) and entorhinal cortex
(ec). The region superior to CA1 represent the stratum lacunosummoleculare and
statum radiatum. Note that in the control hippocampus, the molecular layer and
stratum radiatum fiber orientations paralleled the apical dendrites of granule cells
and pyramidal neurons respectively. In the epileptic hippocampus, the CA1 subfield
pyramidal cell layer is notably lost relative to the control. The architecture of the
dentate gyrus is also notably altered with more evidence of crossing fibers. Future
investigations employing this method should improve our understanding of normal
and pathologically altered neuroanatomy in regions of complex fiber architecture
such as the hippocampus and entorhincal cortex.
5 Conclusion
In this paper, we presented a new mathematical model for the diffusionweighted
MR signals obtained from a single voxel. According to our model, the signal is gen
erated by a continuous distribution of diffusion tensors, where the relevant distribu
tion is a Wishart distribution. In this case, the MR signal was shown to be a Laplace
transform of this distribution defined on the manifold of symmetric positivedefinite
tensors. We presented an explicit form of the expected MR signal attenuation given
by a Rigauttype asymptotic fractal formula. This form of the signal attenuation
has the correct asymptotic dependence of the signal values on the diffusion gra
dient strength. Moreover, the angular dependence of the expected MR signal is
different from the angular dependence implied by traditional DTI. The simulations
of diffusion inside cylindrical boundaries suggested that the principal eigenvec
tors of the diffusion tensors obtained from the proposed model are more accurate
than those implied by traditional DTI. Using this new model in conjunction with a
deconvolution approach, we presented an efficient estimation scheme for the dis
tinct fiber orientations and the water molecule displacement probability functions
(a)
(b)
Fig. 8. Reference (SO) images of: (a) control rat brain (b) epileptic rat brain. The rectangle
regions are the hippocampi.
at each voxel in a HARDI data set. Both synthetic and real data sets were used to
depict the performance of the proposed algorithms. Comparisons with competing
methods from literature depicted our model in a favorable light.
Appendix
A.1 Laplace trnur'i'ilm on P,
For definition of Laplace transforms on P,, we follow the notations in (Terras,
1985).
Hilus
CA1
wif I ** $4146i ft/B a#a
I I II AM Li L
Fig. 9. Probability maps of coronally oriented GA images of a control (a) and epileptic
hippocampus (b). In the control hippocampus, the molecular layer and stratum radiatum
fiber orientations paralleled the apical dendrites of granule cells and pyramidal neurons
respectively, whereas in the stratum lacunosum, molecular orientations paralleled Schaffer
collaterals from CA1 neurons. In the epileptic hippocampus, the overall architecture is
notably altered; the CA1 subfield is lost, while an increase in crossing fibers can be seen in
the hilus and dentate gyrus (dg). Increased crossing fibers can also be seen in the entorhinal
cortex (ec). Fiber density within the statum lacunosum molecular and statum radiale is
also notably reduced, although fiber orientation remains unaltered.
Definition 3 The Laplace ticn'i'liim off : P, C, denoted by J2f, at the sym
metric matrix Z C is defined by Herz (1955):
tff(Z) J f (Y) exp [trace(YZ)] dY, where dY = 7 dy~, 1 < i < j < n.
1..
(17)
Foir oa mrientatioi nice function f, the integral above converges in the right half
plane, Re(Z) > Xo (Re() denotes the real part of Z), meaning that Re(Z) cXo he
collateral from CAI neurons. In the epileptic hippocampus, the overall architecture is
cortex (ec). Fiber density within the status lacunosum molecular and status radiale is
Definition 3 The Laplace tj%,i m of f : ',,, C, denoted by f, at the sym
plane,/Re (Z) > Xo (Re (Z) denotes the real Part of Z), meaning that Rle(Z) Xo G
PP, and the inversion formula for ;hii Laplace tiro\'i'i,n is:
(2ri)'n(n+1)/2 [ Yf(Z) exp [trace(YZ)] dZ f (Y), for Ye c2,
JReZ=Xo 0, otherwise.
(18)
Here dZ = n dzi and the integral is over symmetric matrices Z with fixed real
part.
If f is the density function of some probability measure F on 'P, with respect to
the dominating measure dY, i.e. dF(Y) = f(Y)dY, then Eq.(17) also defines the
Laplace transform of the probability measure F on P,< which is denoted by YF.
References
Alexander, D. C., 2005. Maximum entropy spherical deconvolution for diffusion MRI. In:
Christensen, G. E., Sonka, M. (Eds.), IPMI. Vol. 3565 of Lecture Notes in Computer
Science. Springer, pp. 7687.
Alexander, D. C., Barker, G. J., Arridge, S. R., August 2002. Detection and modeling of
nonGaussian apparent diffusion coefficient profiles in human brain data. Magn. Reson.
Med. 48 (2), 331340.
Anderson, A. W., 2005. Measurement of fiber orientation distributions using high angular
resolution diffusion imaging. Magn. Reson. Med. 54 (5), 11941206.
Anderson, T. W., 1958. An Introduction to Multivariate Statistical Analysis. John Wiley
and Sons.
Assaf, Y., Freidlin, R. Z., Rohde, G. K., Basser, P. J., 2004. New modeling and experimental
framework to characterize hindered and restricted water diffusion in brain white matter.
Magn. Reson. Med. 52 (5), 965978.
Basser, P. J., Mattiello, J., LeBihan, D., 1994. Estimation of the effective selfdiffusion
tensor from the NMR spin echo. J. Magn. Reson. B 103, 247254.
Basser, P. J., Pajevic, S., Peierpaoli, C., J, D., Aldroubi, A., October 2000. In vivo fiber
tractography using DTMRI data. Magn. Reson. Med. 44 (4), 625632.
Behrens, T., Woolrich, M., Jenkinson, M., JohansenBerg, H., Nunes, R., Clare, S.,
Matthews, P., Brady, J., Smith, S., 2003. Characterization and propagation of uncertainty
in diffusionweighted mr imaging. Magn. Reson. Med. 50 (2), 10771088.
Callaghan, P., 1991. Principles of Nuclear Magnetic Resonance Microscopy. Clarendon
Press, Oxford.
Chefd'Hotel, C., Tschumperl6, D., Deriche, R., Faugeras, O. D., 2002. Constrained flows
of matrixvalued functions: Application to diffusion tensor regularization. In: ECCV (1).
pp.251265.
Cleveland, G. G., Chang, D. C., Hazlewood, C. E, Rorschach, H. E., 1976. Nuclear mag
netic resonance measurement of skeletal muscle: anisotropy of the diffusion coefficient
of the intracellular water. Biophys J 16 (9), 10431053.
Conturo, T. E., Lori, N. E, Cull, T. S., Akbudak, E., Snyder, A. Z., Shimony, J. S., McK
instry, R. C., Burton, H., Raichle, M. E., 1999. Tracking neuronal fiber pathways in the
living human brain. Proc Natl Acad Sci 96, 1042210427.
Descoteaux, M., Angelino, E., Fitzgibbons, S., Deriche, R., 2006. A fast and robust odf es
timation algorithm in qball imaging. In: International Symposium on Biomedical Imag
ing: From Nano to Macro. pp. 8184.
Frank, L., 2002. Characterization of anisotropy in high angular resolution diffusion
weighted MRI. Magn. Reson. Med. 47 (6), 10831099.
Herz, C. S., 1955. Bessel functions of matrix argument. The Annals of Mathematics 61 (3),
474523.
Hess, C. P., Mukherjee, P., Han, E. T., Xu, D., Vigneron, D. B., 2006. Qball reconstruction
of multimodal fiber orientations using the spherical harmonic basis. Magn. Reson. Med.
56 (1), 104117.
Jansons, K. M., Alexander, D. C., 2003. Persistent angular structure: new insights from
diffusion MRI data. Inverse Problems 19, 10311046.
K6pf, M., Metzler, R., Haferkamp, O., Nonnenmacher, T. F., 1998. NMR studies of anoma
lous diffusion in biological tissues: Experimental observation of L6vy stable processes.
In: Losa, G. A., Merlini, D., Nonnenmacher, T. E, Weibel, E. R. (Eds.), Fractals in Biol
ogy and Medicine. Vol. 2. Birkhiiuser, Basel, pp. 354364.
LeBihan, D., Breton, E., Lallemand, D., Grenier, P., Cabanis, E., LavalJeantet, M., 1986.
MR imaging of intravoxel incoherent motions: Application to diffusion and perfusion in
neurologic disorders. Radiology 161, 401407.
Letac, G., Massam, H., 1998. Quadratic and inverse regressions for Wishart distributions.
The Annals of Statistics 26 (2), 573595.
Mori, S., Crain, B. J., Chacko, V. P., van Zijl, P. C. M., 1999. Threedimensional tracking
of axonal projections in the brain by magnetic resonance imaging. Ann Neurol 45, 265
269.
Moseley, M. E., Cohen, Y., Kucharczyk, J., Mintorovitch, J., Asgari, H. S., Wendland,
M. F, Tsuruda, J., Norman, D., 1990a. Diffusionweighted MR imaging of anisotropic
water diffusion in cat central nervous system. Radiology 176 (2), 439445.
Moseley, M. E., Cohen, Y., Mintorovitch, J., Chileuitt, L., Shimizu, H., Kucharczyk, J.,
Wendland, M. F, Weinstein, P. R., 1990b. Early detection of regional cerebral ischemia
in cats: comparison of diffusion and T2weighted MRI and spectroscopy. Magn Reson
Med 14, 330346.
Murihead, R. J., 1982. Aspects of multivariate statistical theory. John Wiley & Sons.
Ozarslan, E., Basser, P. J., Shepherd, T. M., Thelwall, P. E., Vemuri, B. C., Blackband,
S. J., 2006a. Observation of anomalous diffusion in excised tissue by characterizing the
diffusiontime dependence of the MR signal. J Magn Reson 183 (2), 315323.
Ozarslan, E., Mareci, T. H., 2003. Generalized diffusion tensor imaging and analytical rela
tionships between diffusion tensor imaging and high angular resolution diffusion imag
ing. Magn. Reson. Med. 50 (5), 955965.
Ozarslan, E., Shepherd, T. M., Vemuri, B. C., Blackband, S. J., Mareci, T. H., 2006b. Res
olution of complex tissue microarchitecture using the diffusion orientation transform
(DOT). Neurolmage 31, 10861103.
Ozarslan, E., Vemuri, B. C., Mareci, T., 2004. Fiber orientation mapping using generalized
diffusion tensor imaging. In: International Symposium on Biomedical Imaging: From
Nano to Macro. pp. 10361038.
Ozarslan, E., Vemuri, B. C., Mareci, T. H., 2005. Generalized scalar measures for diffusion
MRI using trace, variance, and entropy. Magn. Reson. Med. 53 (4), 866876.
Rigaut, J. P., 1984. An empirical formulation relating boundary lengths to resolution in
specimens showing 'nonideally fractal' dimensions. J Microsc 133, 4154.
Sen, P. N., Hiirlimann, M. D., de Swiet, T. M., 1995. DebyePorod law of diffraction for
diffusion in porous media. Phys Rev B 51 (1), 601604.
Siderman, O., J6nsson, B., 1995. Restricted diffusion in cylindirical geometry. J. Magn.
Reson. B A (117), 9497.
Stejskal, E. O., Tanner, J. E., 1965. Spin diffusion measurements: Spin echoes in the pres
ence of a timedependentfield gradient. J. Chem. Phys. 42 (2), 288292.
Terras, A., 1985. Harmonic analysis on symmetric spaces and applications. Springer.
Tournier, J.D., Calamante, F, Gadian, D. G., Connelly, A., November 2004. Direct esti
mation of the fiber orientation density function from diffusionweighted MRI data using
spherical deconvolution. Neurolmage 23 (3), 11761185.
Tournier, J.D., Calamente, F, Connelly, A., 2006. Improved characterisation of crossing
fibres: spherical deconvolution combined with Tikhonov regularization. In: Proceedings
of the ISMRM 14th Scientific Meeting and Exhibition. Seattle, Washington.
Tuch, D. S., 2004. Qball imaging. Magn. Reson. Med. 52 (6), 13581372.
Tuch, D. S., Reese, T. G., Wiegell, M. R., Makris, N., Belliveau, J. W., Wedeen, V. J.,
2002. High angular resolution diffusion imaging reveals intravoxel white matter fiber
heterogeneity. Magn. Reson. Med. 48 (4), 577582.
Tuch, D. S., Reese, T. G., Wiegell, M. R., Wedeen, V. J., December 2003. Diffusion MRI
of complex neural architecture. Neuron 40, 885895.
Tuch, D. S., Weisskoff, R. M., Belliveau, J. W., Wedeen, V. J., 1999. High angular resolution
diffusion imaging of the human brain. In: Proc. of the 7th ISMRM. Philadelphia, p. 321.
von dem Hagen, E., Henkelman, R., 2002. Orientational diffusion reflects fiber structure
within a voxel. Magn. Reson. Med. 48 (3), 454459.
Wampler, C. W., 1986. Manipulator inverse kinematic solution based on damped least
squares solutions. IEEE Trans. Systems, Man and Cybernetics 16 (1), 93101.
Wang, Z., Vemuri, B. C., Chen, Y., Mareci, T. H., 2004. A constrained variational principle
for direct estimation and smoothing of the diffusion tensor field from complex DWI.
IEEE Trans. Med. Imaging 23 (8), 930939.
Wedeen, V. J., Reese, T., Tuch, D. S., Weigl, M. R., Dou, J.G., Weiskoff, R., Chesler, D.,
2000. Mapping fiber orientation spectra in cerebral white matter with fourier transform
diffusion mri. In: Proc. of the 8th ISMRM. Denver, p. 82.
Weickert, J., H.Hagen (Eds.), 2005. Visualization and Processing of Tensor Fields.
Springer.
Wishart, J., 1928. The generalized product moment distribution in samples from a normal
multivariate population. Biometrika 20, 3252.
