UFDC Home |   Help  |   RSS

 Group Title: Department of Computer and Information Science and Engineering Technical Reports Title: A Novel tensor distribution model for the diffusion weighted MR signal
CITATION PDF VIEWER THUMBNAILS PAGE IMAGE ZOOMABLE
Full Citation
STANDARD VIEW MARC VIEW
 Material Information Title: A Novel tensor distribution model for the diffusion weighted MR signal Series Title: Department of Computer and Information Science and Engineering Technical Report ; 06-291 Physical Description: Book Language: English Creator: Jian, BingVemuri, Baba C.Ozarslan, EvrenCarney, PaulMareci, Thomas Publisher: Department of Computer and Information Science and Engineering, University of Florida Place of Publication: Gainesville, Fla. Copyright Date: 2006
 Record Information Bibliographic ID: UF00095643 Volume ID: VID00001 Source Institution: University of Florida Holding Location: University of Florida Rights Management: All rights reserved by the source institution and holding location.

2006291 ( PDF )

Full Text

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 06-291

Abstract

Diffusion MRI is a non-invasive 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 diffusion-weighted 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 Rigaut-type 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 voxel-by-voxel 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

Diffusion-weighted 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 relaxation-weighted 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 white-matter 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

A widely accepted method for diffusion-weighted imaging studies is to use appar-
ent difl,'itii coefficient (ADC) profiles which is governed by the Stejskal-Tanner
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 (DT-MRI 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 b-factor, 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 white-matter 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 DT-MRI 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 three-dimensional Cartesian lattice, the
q-space 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 time-intensive 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 Stejskal-Tanner 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 rank-2 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 q-ball imaging (QBI) method, in which the radial integral of
the displacement PDF is approximated by the spherical Funk-Radon transform. The
end result obtained by the radial projection of the three-dimensional 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 Funk-Radon 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 three-dimensional wave
vector by fitting higher-order 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 PAS-MRI was pro-
posed by Jansons and Alexander (2003). This method computed the so-called the
persistent angular structure (PAS), which is a function on a fixed spherical shell in
three-dimensions, by assuming its Fourier transform best fits the measurements and
incorporating a maximum-entropy 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 mono-exponential decay assumption of
the MR signal, and hence can be avoided by using the extension of the transform
to multi-exponential 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 multi-compartmental models or
multiple-fiber population models have been used to model the diffusion-attenuated
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 multi-compartment models. First, the
number of such compartments has to be pre-specified, presenting a model-selection
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 model-based approaches which extend to the multiple-fiber
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
model-selection and fitting problems as the mixture models discussed above.

To address the model-selection 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 real-valued 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 multiple-compartment mod-
els (Tuch et al., 2002; Assaf et al., 2004) since both previous single-tensor and
multiple-compartment 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 Rigaut-type 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 positive-definite 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
ill-posed 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 q-space. 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 positive-definite
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 Rigaut-type 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 matrix-valued 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 = 2-p/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(E-ZY)) 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, + OE-P' 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 Rigaut-type 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 power-law 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 diffusion-weighted MR data before (K6pf et al., 1998), until now, there was
no rigorous derivation of the Rigaut-type expression used to explain the MR sig-
nal behavior as a function of b-value. Therefore, this derivation may be useful in
understanding the apparent fractal-like 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 non-integer, 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 log-log plots-hence the phrase "Rigaut-type".

analog of Debye-Porod law of diffraction (Sen et al., 1995) ensures that in three-
dimensional space the signal should have the asymptotic behavior, S(q) ~ q-4
Since b oc q2 a reasonable choice for p is 2.
1 .0000 .

0 000i ',. .. .. "' ""*.' .-.-.: rigaul-.ype
0,1000 i '. t a urnp-tli 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 Rigaut-type 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 rank-2 tensors,
the corresponding multi-exponential 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 mono-exponential model where the diffu-
sion tensor is taken to be the expected value of the Wishart distribution. The right
extreme solid curve is the Rigaut-type decay derived from (11). Note that the tail of
the solid curve is linear indicating the power-law 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 Rigaut-type 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 Gaussian-like, but as the number of tensors increases, the at-
tenuation curve converges to a Rigaut-type 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 mono-exponential 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 bi-exponential and multi-exponential
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 diffusion-weighted 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 Multi-fiber reconstruction using deconvolution

However, the single Wishart model can not resolve the intra-voxel 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 Debye-Porod law of diffraction (Sen et al., 1995) in three-dimensional 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 white-matter 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 under-determined 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 one-time 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 ill-conditioned 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 real-valued 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 F-test 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 single-fiber 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 Stejskal-Tanner 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 Rician-distributed 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 Stejskal-Tanner
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
3-fiber 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 Q-ball 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 gradient-based 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 3-fibers (b = 1500s/mm2). Orientations: azimuthal angles
01 = 30, 2 = {20,100}, 03 = {20, 75,135}; polar angles were all 90. Top: Q-ball 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 Rician-distributed noise into the
magnitude of the MR signal.
* For all 1-, 2- and 3-fiber 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 Q-ball 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 Quasi-Newton 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, perfusion-fixed rat ner-

J1~-~.~ l,,s 4k4'~i k4
4^ <'M^4r
S-i rA ^

(a) ODF from QBI

kk)V _1'kA.t
-4( 4-1 A., A-

4 i -k41 +k8~

(b) Proposed method

Fig. 5. Resistance to noise (2-fibers, 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 (3-fibers, o = 0.04)
vous tissue was acquired at 14. IT using Bruker Avance imaging systems. A diffusion-
weighted spin echo pulse sequence was used. Diffusion-weighted images were
acquired along 46 directions with a b-value 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

FC-43, 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 multiple-slice 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 spin-echo, pulsed-field-gradient 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 long-axis of the brain (slices progressed in the dorsal-ventral direc-

tion). These slices have a field-of-view 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.
Diffusion-weighted 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 diffusion-weighting 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 lacunosum-moleculare 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 diffusion-weighted
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 positive-definite
tensors. We presented an explicit form of the expected MR signal attenuation given
by a Rigaut-type 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. 76-87.
Alexander, D. C., Barker, G. J., Arridge, S. R., August 2002. Detection and modeling of
non-Gaussian apparent diffusion coefficient profiles in human brain data. Magn. Reson.
Med. 48 (2), 331-340.
Anderson, A. W., 2005. Measurement of fiber orientation distributions using high angular
resolution diffusion imaging. Magn. Reson. Med. 54 (5), 1194-1206.
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), 965-978.
Basser, P. J., Mattiello, J., LeBihan, D., 1994. Estimation of the effective self-diffusion
tensor from the NMR spin echo. J. Magn. Reson. B 103, 247-254.
Basser, P. J., Pajevic, S., Peierpaoli, C., J, D., Aldroubi, A., October 2000. In vivo fiber
tractography using DT-MRI data. Magn. Reson. Med. 44 (4), 625-632.
Behrens, T., Woolrich, M., Jenkinson, M., Johansen-Berg, H., Nunes, R., Clare, S.,
Matthews, P., Brady, J., Smith, S., 2003. Characterization and propagation of uncertainty
in diffusion-weighted mr imaging. Magn. Reson. Med. 50 (2), 1077-1088.
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 matrix-valued functions: Application to diffusion tensor regularization. In: ECCV (1).
pp.251-265.
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), 1043-1053.
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, 10422-10427.
Descoteaux, M., Angelino, E., Fitzgibbons, S., Deriche, R., 2006. A fast and robust odf es-
timation algorithm in q-ball imaging. In: International Symposium on Biomedical Imag-
ing: From Nano to Macro. pp. 81-84.
Frank, L., 2002. Characterization of anisotropy in high angular resolution diffusion
weighted MRI. Magn. Reson. Med. 47 (6), 1083-1099.
Herz, C. S., 1955. Bessel functions of matrix argument. The Annals of Mathematics 61 (3),
474-523.
Hess, C. P., Mukherjee, P., Han, E. T., Xu, D., Vigneron, D. B., 2006. Q-ball reconstruction
of multimodal fiber orientations using the spherical harmonic basis. Magn. Reson. Med.
56 (1), 104-117.
Jansons, K. M., Alexander, D. C., 2003. Persistent angular structure: new insights from
diffusion MRI data. Inverse Problems 19, 1031-1046.
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. 354-364.
LeBihan, D., Breton, E., Lallemand, D., Grenier, P., Cabanis, E., Laval-Jeantet, M., 1986.
MR imaging of intravoxel incoherent motions: Application to diffusion and perfusion in
Letac, G., Massam, H., 1998. Quadratic and inverse regressions for Wishart distributions.
The Annals of Statistics 26 (2), 573-595.
Mori, S., Crain, B. J., Chacko, V. P., van Zijl, P. C. M., 1999. Three-dimensional 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. Diffusion-weighted MR imaging of anisotropic
water diffusion in cat central nervous system. Radiology 176 (2), 439-445.
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 T2-weighted MRI and spectroscopy. Magn Reson
Med 14, 330-346.
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
diffusion-time dependence of the MR signal. J Magn Reson 183 (2), 315-323.
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), 955-965.
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, 1086-1103.
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. 1036-1038.
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), 866-876.
Rigaut, J. P., 1984. An empirical formulation relating boundary lengths to resolution in
specimens showing 'non-ideally fractal' dimensions. J Microsc 133, 41-54.
Sen, P. N., Hiirlimann, M. D., de Swiet, T. M., 1995. Debye-Porod law of diffraction for
diffusion in porous media. Phys Rev B 51 (1), 601-604.
Siderman, O., J6nsson, B., 1995. Restricted diffusion in cylindirical geometry. J. Magn.
Reson. B A (117), 94-97.
Stejskal, E. O., Tanner, J. E., 1965. Spin diffusion measurements: Spin echoes in the pres-
ence of a time-dependentfield gradient. J. Chem. Phys. 42 (2), 288-292.
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 diffusion-weighted MRI data using
spherical deconvolution. Neurolmage 23 (3), 1176-1185.
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. Q-ball imaging. Magn. Reson. Med. 52 (6), 1358-1372.
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), 577-582.
Tuch, D. S., Reese, T. G., Wiegell, M. R., Wedeen, V. J., December 2003. Diffusion MRI
of complex neural architecture. Neuron 40, 885-895.
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), 454-459.
Wampler, C. W., 1986. Manipulator inverse kinematic solution based on damped least-
squares solutions. IEEE Trans. Systems, Man and Cybernetics 16 (1), 93-101.
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), 930-939.
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, 32-52.