Resolution of Complex Tissue Microarchitecture using the
Diffusion Orientation Transform (DOT)*
Evren O. ,ii ,,' Timothy M. Shepherd, Baba C. V\ iii',
Stephen J. Blackbandb,c, Thomas H. M .i i.
a Department of Computer and Inf.i ,iin:ii.. Science and Egiyii iiiy.
Uo,;.',. .:7,/ of Florida, Gainesville, FL 32611, USA
b Department of Neuroscience, U',.::' i .'/; of Florida, Gainesville, FL 32610, USA.
c National High Magnetic Field Lil,..,'l..I '. Tallahassee, FL 32310, USA.
d Department of Biochemistry and Molecular B.: /...;~ .
U',.:: .:// of Florida, Gainesville, FL 32610, USA.
*E. Ozarslan, T. M. Shepherd, B. C. Vemuri, S. J. Blackband, T. H. Mareci, Resolution of Complex Tis
sue Microarchitecture using the L'D i.. . Orientation Transform (DOT), University of Florida, Gainesville, FL,
Technical Report TR05004, July 2005.
tCorresponding author. Email: evren@cise.ufl.edu
Abstract
This article describes an accurate and fast method for fiber orientation mapping using mul
tidirectional diffusionweighted magnetic resonance (MR) data. This novel approach utilizes
the Fourier transform relationship between the water displacement probabilities and diffusion
attenuated MR signal expressed in spherical coordinates. The radial part of the Fourier integral
is evaluated analytically under the assumption that MR signal attenuates exponentially. The
values of the resulting functions are evaluated at a fixed distance away from the origin. The
spherical harmonic transform of these functions yields the Laplace series coefficients of the
probabilities on a sphere of fixed radius. Alternatively, probability values can be computed
nonparametrically using Legendre polynomials. Orientation maps calculated from excised rat
nervous tissue data demonstrate this technique's ability to accurately resolve crossing fibers in
anatomical regions such as the optic chiasm. This proposed methodology has a trivial extension
to multiexponential diffusionweighted signal decay. The developed methods will improve the
reliability of tractography schemes and may make it possible to correctly identify the neural
connections between functionally connected regions of the nervous system.
Keywords: MRI, tensor, anisotropy, HARDI, Fourier, spherical harmonics
1 Introduction
The diffusional attenuation of MR signal in pulsed field gradient experiments (Sj' i1:. and
Tanner, 1'.".) has been exploited to characterize diffusional anisotropy in fibrous tissues like
muscle (e.g. C'I 1. I !1.I et al., 1976) and whitematter in animal (e.g. Moseley et al., 1990) and
human (e.g C'!,. x. i. it et al., 1990) nervous tissue. When the narrow pulse condition is met,
i.e. the duration of the applied diffusion sensitizing gradients (6) is much smaller than the time
between the two pulses (A), the fundamental relationship between the MR signal attenuation
E(q) and average displacement probabilities P(R) is given by a Fourier integral (Callaghan,
1991):
P(R) E(q) exp(27iq R) dq (1)
where R is the displacement vector and q is the reciprocal space vector defined by q = 76G/2r,
where 7 is the gyromagnetic ratio and G is the gradient vector.
Diffusional anisotropy is wellreflected in the water displacement probabilities and it is ex
pected that in fibrous tissues, the orientations specified by large displacement probabilities will
coincide with the fiber orientations. One could in principle estimate these displacement proba
bility functions by using Eq. 1 and the fast Fourier transform (FFT), however this would require
data points all across the space spanned by the diffusion gradients (or q vectors). This qspace
approach would require very high gradient strengths and long acquisition times that are difficult
to achieve in clinical settings (Basser, 2002). Although attempts have been made to acquire such
datasets in vivo (Wedeen et al., 2000), the results typically suffer from undersampled qspace
and sacrificed spatial resolution.
More than a decade ago, Basser et al. (1994a,b) introduced an imaging method called diffu
sion tensor imaging (DTI) that replaced the apparent diffusion coefficients that had been cal
culated in diffusionweighted imaging studies with a symmetric, positivedefinite, secondorder
tensor. This model required only 7 diffusionweighted images with clinically feasible diffusion
gradient strengths. This approach enabled simple estimation of diffusional anisotropy, and pre
dicted a fiber orientation specified by the principal eigenvector of the diffusion tensor. Despite
its modest requirements, the results achieved using DTI have been very successful in regions of
the brain and spinal cord with substantial whitematter coherence and have enabled the map
ping of some anatomical connections in the central nervous system (e.g. Conturo et al., 1999;
Mori et al., 1999; Basser et al., 2000).
DTI assumes a displacement probability characterized by an oriented Gaussian probability
distribution function (PDF) whose covariance matrix is proportonal to the diffusion tensor.
Such a PDF has only one orientational mode and as such, can not resolve more than one fiber
orientation inside a voxel. This shortcoming of DTI has prompted interest in the development
of more sophisticated models. Tuch et al. (1999) introduced a high angular resolution diffusion
imaging (HARDI) method that ,: :_. I. 1 the apparent diffusion coefficients could be evaluated
along many orientations independently without fitting a s, I Il" function to the data, i.e. using
the Si. ii: iTanner expression (S'. i1. I. 1'"..)
S(u) So ebD(u) (2)
where u is a unit vector specifying the direction of the diffusion sensitizing gradient, S(u) is the
signal value on a sphere in qspace whose radius is related to the diffusion weighting factor b
(where b =47q2t and t = A 6/3 is the effective diffusion time), and So is the signal value at
b = 0. The result is an angular distribution of apparent diffusivities, D(u) herein referred to as
the diffusivity profile. It has been shown that the diffusivity profile has a complicated structure
in voxels with orientational heterogeneity (von dem Hagen and Henkelman, 2002; Tuch et al.,
2002). Several studies proposed to represent the diffusivity profile using a spherical harmonic
expansion (Frank, 2002; Alexander et al., 2002). A schematic description of this approach is
given below:
LS
D(u) aim, (3)
SHT
where SHT and LS stand for spherical harmonic transform and Laplace series respectively.
However, one major difficulty with employing HARDI in studies involving orientation map
ping has been that the peaks of the diffusivity profile do not necessarily yield the orientations of
the distinct fiber populations (Figure 1). Ozarslan and Mareci (2003) have shown that the (SHT)
approach could be seen as a generalization of DTI since the coefficients of the Laplace series
(obtained from the SHT of the diffusivity profile) are related to the components of higherorder
Cartesian tensors. Later, Ozarslan et al. (2004a,b) proposed to use the higherorder Cartesian
tensors to generate signal values (assuming exponential attenuation) on the threedimensional
qspace, and evaluated an FFT to approximate the displacement probabilities. Jansons and
Alexander (2003) proposed a method to calculate a displacement probability map from HARDI
datasets by enforcing the unusual condition that the probabilities are nonzero only on a spher
ical shell. Although the results are encolur i:in:. both of these schemes are computationally
expensive.
Another generalization of DTI that employs higherorder Cartesian tensors was proposed
by Liu et al. (2003). This approach necessitates sampling of qspace in several spherical shells,
undesirably increasing the required number of acquisitions. Also, it is difficult to reliably extract
the phase of the MR signal required by this scheme. Tuch et al. (2003) proposed a method
in which the radial integral of the displacement PDF is approximated by the spherical Radon
transform. However, the end result is the radial projection of the threedimensional displacement
PDF via a line integral and not the real probabilities as would be obtained through the radial
integral with measure r2 dr. More importantly, this scheme provides only an approximation
probability
Figure 1: Apparent diffusivity (left column) and displacement probability (right column) profiles
calculated from simulations of 1, 2 and 3 fiber systems (top to bottom). Black lines depict the
exact orientations of the simulated fibers specified by the azimuthal angles 01 = 300, 02
{200, 100} and 03 {= 200, 550, 100} for the 1, 2 and 3 fiber systems respectively. Polar angles
for all fibers were taken to be 900 so that all fibers lie on the image plane. The peaks of the
diffusivity profile do not necessarily yield the orientations of the distinct fiber populations. This
can sometimes lead to erroneous fiber structure interpretation from HARDI data.
to the true radial integral because the result is a convolution of the probability values with
a 0th order Bessel function (Tuch, 2004) that gives rise to an undesirable ,iii_ Iii, i.''Ii.. "
of the pr'!1 .1.il liv along one direction with probabilities from other directions. Finally, there
have been several studies that have modeled diffusion using multicompartmental models. These
studies assume distinct fiber populations with no exchange. Moreover, the number of such
compartments has to be prespecified (Inglis et al., 2001; Parker and Alexander, 2003; Maier
et al., 2004; Assaf et al., 2004) or the signal from each fiber population is undesirably forced to
have prespecified attributes (such as anisotropy) (Tournier et al., 2004).
In this work, we introduce a new method, called the diffusion orientation transform (DOT),
that describes how the diffusivity profiles can be transformed into pr'l. 1.;i1lilv profiles. Our
method is based on the HARDI acquisition scheme and can be extended to more general acqui
sition strategies. We express Eq. 1 in spherical coordinates, then under the monoexponential
attenuation assumption, evaluate the radial part of the integral analytically. The probability
diffusivity
values on a fixed radius can be reconstructed either directly or parametrically in terms of a
Laplace series. We prove that this expansion converges to the true probability profile. Our
technique can be regarded as a transformation of diffusivity to probability profiles whose peaks
correspond to distinct fiber orientations (Figure 1). Our method is robust and fast. Although we
present results on excised, chemicallyfixed rat nervous tissue, the requirements of our method
make it suitable for the clinical environment. We discuss our assumption that the MR signal
decays monoexponentially and further demonstrate that a trivial extension to multiexponential
attenuation results in improved reconstruction of the probabilities.
2 Theory
In this section we show that starting from the signal attenuations of a HARDI acquisition, it
is possible to calculate the orientation maps without the need to fit a particular model to the
data. We achieve this by two different approaches.
2.1 Parametric Reconstruction
The Fourier transform that relates the signal attenuation to the water displacement probability
(Eq. 1) can be written in spherical coordinates. This is a consequence of the pointwise convergent
expansion of the plane wave in spherical coordinates (Schwabl, 1989) given by
00 I
2iq.R 47r E (i)zj,(27rqr) Yim(u)* Ym(r) (4)
1=0 m=1
where q = qu and R = rr, with q = q and r = IR. Inserting this expression into Eq. 1, we
get
00 1
P(Ror) (i) m(r) j du I m(u)* I(u) (5)
1=0 m=l
where
11(u) 47 dqq2 j(27qRo) exp(472q2tD(u)) (6)
Here r was set to a particular radius Ro and it is assumed that signal attenuates along each
radial line in qspace as described by the Sj' ii. 1Tanner relationship given in Eq. 2. Note
that the function P(Ror) is not the isosurface of the threedimensional displacement probability
function, but it is the probability of finding the particle, initially at the origin, at the point Ror,
i.e. we will be interested in the probability values on a sphere of radius Ro. The integral in Eq.
6 can be evaluated and it is given by
RI F('T3) 1+3 3 R_
I2(u) 2a + I 2 2' 4D(u) (7)
21+3 73/2 (D(u)t)(l+3)/2 ]F( + 3/2) W 2 4D(u)t
Table 1: Az(u) and Bl(u) functions up to I = 8. In this table, /3 stands for /(u).
1 Ai(u) B (u)
0 1 0
2 (1 +6/2) 3
4 1 + 20/32 + 210/34 (1 14/32)
6 (1 + 42/32 + 15754 6) (1 36/32 + 396/34)
8 1 + 72 2 + 9/4 + 45045/36 315 (1 66/32 + 1716/34 17160/36)
+675675/3 
where 1Fi is the confluent hypergeometric function of the first kind (see Appendix A). Using the
recurrence relations of the confluent hypergeometric functions provided in Eqs. 30 iteratively,
these functions can be written as the sum of two terms, one of them being proportional to
2 4D(u)t) 2 4D(u)t"U
1Fi (; ; where the other term will be proportional to iF 2; 4; )). Using
Eqs. 31 and 32, it can be seen that these functions are proportional to the error function and
Gausssian respectively.
Therefore, the resulting expression is given by
exp(/3(u)2/4) erf(/3(u)/2)
I1(u) = Al(u) (4 (u)/2 + B(u) (8)
(4 D(u) t)3/2 4x R3
where
/3(u) = (9)
gDu) t
Az(u) and Bl(u) functions up to I = 8 are given in Table 1. Note that throughout the paper,
only the evenorder terms will be included as a consequence of the antipodal symmetry of the
diffusivity profiles as well as displacement PDFs. The derivation of the particular forms for the
Al(u) and Bl(u) for arbitrary (even) values of I are provided in Appendix B.
In Figure 2a, we plot the I values as a function of Ro calculated with double precision using
Eq. 8 where D = 1.5 x 103mm2/s and t = '"n, Very large values taken by the higherorder
terms near the origin are due to roundoff errors. However, this is not a big concern because we
will be mostly interested in the values of this function in the 10 20pm range. Note that the
contribution from higherorder terms is rapidly collapsing for Ro values in this range. Shown in
Figure 2b are the curves generated by keeping Ro fixed at the value of 15/m and varying the
diffusion coefficients between 3 x 104 and 3 x 103 mm2/s. This plot indicates the nontrivial
manner in which an angular diffusivity profile influences the I and hence the probability values.
In order to estimate the probability on the surface of a sphere of radius Ro, we go back to
Eqs. 5 and 6. Since Ii(u) is a function of orientation, we can expand it in a Laplace series, i.e.
00 1'
I1(u) = = auE m c Yim,(u), (10)
l'=0 m'=l'
7
8xO15  =6 8x10 5 =6 
... ;  l=8
6xiO5 6xJO5
2x105 xl 105 
4 10 0 7 
0 10 20 30 40 0.0 0.5 1.0 1.5 2.0 2.5 3.0
Ro (4m) D (pm /ms)
Figure 2: Dependence of the radial integral I, on Ro (top) and on diffusivity (bottom). The
curves are drawn for 1 values ranging from 0 to 8.
where the coefficients aull/, are given by a SHT,
alrn =/ Ym(u)*1(u)du. (11)
Comparing the integration over u in Eq. 5 with the expression in Eq. 11, it can be seen that
P(Ror) has the Laplace series expansion
00 l
P(Ror)= 1 %., YI,(r) (12)
l=0 m=l
with
/,.. (i)1Caum1=. (1)1/2al (13)
where in the last step we have used the fact that 1 is even. The convergence of the resulting series
in Eq. 5 to the desired probability value is proved in Appendix C. Note that the coefficients of
this Laplace series for a particular value of I come from the lth order Laplace series coefficients
of 17(u).
2.1.1 Implementation Aspects
In summary, given the HARDI data, the estimation of the probability of finding the particle at
the point Ror away from the origin involves the following steps:
1. Compute the diffusivity D(u) along each direction using Eq. 2.
2. Then compute 17(u) using Eq. 7 or Eq. 8 with Table 1.
3. For each 1, compute allm, the 1th order spherical harmonic transform of Ii(u).
4. Then evaluate Eq. 12.
1x104
x10 4
Implementation of the items 1, 2 and 4 above are trivial. Our data acquisition scheme
involves sampling the sphere on the vertices of a tessellated icosahedron. With this method
46 or 81 points are sampled on the unit hemisphere from second or thirdorder tessellations
respectively. Following Ritchie and Kemp (1999), we compute the spherical harmonic transform
given in Eq. 11 by discretizing the integrals on the sphere with integration weights calculated
from the areas of the ..., :_ons specified by the dual tessellation. We also exploit the fact that
the probabilities are real. This condition ensures that the expression
Pl() = (1), *, (14)
holds. Thus, it is unnecessary to evaluate the integrals that generate /,.. coefficients with
negative m values. The calculation of the aiim coefficients takes only 25 to 60s for the entire
dataset, depending on the matrix size and number of angular samples, when using a modest
Athlon XP 1800 processor (AMD, Sunnyvale, CA).
Schematic description of the method is described below. Note that this is our revision of Eq.
3 provided in the introduction.
D(u) E.[8] (u) SHT ) all x(1)2 LS P(Ror) (15)
2.2 Nonparametric Reconstruction
An alternative form to the Rayleigh expansion in Eq. 4 is given by
00
27riq.R (i)' (21 + 1) j(27qr)P(u r) (16)
1=0
which is just a consequence of the addition theorem for spherical harmonics provided in Eq.
41. In Eq. 16, Pl is the 1th order Legendre polynomial. Employing this form of the PR 1. ih
expansion in our formalism does not change the radial integral and the probability values are
given by
P(Ror) 7 (i) (21 + 1) dull(u) Pl(u.r)
I=0
Sdu (1)12 l P (u r) 1 h(u) (17)
l=0
with the definition of Ii as in Eqs. 68.
The above expression provides an alternate estimation of the results that could be obtained
from the parametric reconstruction. The schematic description of the nonparametric reconstruc
tion is given by
D(u) Eq.8 I(u) Eq.17 P(Ror) (18)
The above formulation can be easily expressed in matrix form. Suppose that the HARDI
experiment is performed with diffusion sensitizing gradients applied along NG directions. The
direction describing the jth gradient will be represented with the unit vector uj. Similarly,
let ri denote the unit vector describing the ith direction along which the probability will be
estimated where the total number of such directions is NR. Then Eq. 17 can be expressed
simply by
00
Y M ZZ (19)
l=0
where Y is the NR dimensional vector of probabilities. In Eq. 19 the components of the NG
dimensional vector ZZ are given by
(Zi)j I1(uj) (20)
and the components of the NR x No dimensional matrix M1 are given by
(M)ij (_1)1/2 (21 + 1) Pl(u r) (21)
where wj are the integration weights associated with each of the gradient directions. Note
that the matrices MI can be computed once for each gradient sampling scheme. Therefore,
the only computational burden comes from the pixelbypixel estimation of Ii(uj) (which is a
straightforward operation) and the matrix multiplication in Eq. 19.
2.3 Parametric vs. Nonparametric Reconstruction
We have provided two methods for the reconstruction of probability profiles. The parametric
reconstruction enables one to express the probabilities in terms of a Laplace series, whereas the
nonparametric reconstruction provided the probability values directly. It is simpler to imple
ment the latter scheme as no SHT transform is necessary. However, when the Laplace series is
terminated at = Imax, the parametric reconstruction expresses the probability values in terms
of (lmax+1)(max +2)/2 numbers, which are typically much smaller than the number of directions
along which the probabilities are estimated (Np) when one visualizes the pr'! 1.1dilli , surfaces.
This enables more feasible storage of the probability profiles in computer memory.
Using either one of the schemes, once the pr! .1. 1 .illi , values are evaluated along many points,
the following parametrized surface can be visualized (see Fig. 1):
( sin0 cos
X(O, ) =P(Ror) r =P(0, ) sin 0 sin ( (22)
cos 0
where 0 is the polar and ) is the azimuthal angle associated with the unit vector r.
/////
Figure 3: Probability maps estimated on a sphere of radius 8 to 16/m in equal steps of 2pm
(from left to right.) Top row shows these surfaces when there is only one orientation, where the
bottom row shows them from a voxel with two distinct orientations. As the radius of the sphere
on which the probability values are estimated is increased, the two fiber orientations are better
resolved.
3 Simulations
We have applied the scheme described above to the simulations of single fiber and crossing fiber
systems. The simulations employed the exact form of the MR signal attenuation from particles
diffusing inside cylindirical boundaries (S6derman and J6nsson, 1995) with parameters provided
in Ozarslan et al. (2005). As already shown in Figure 1, we have computed the probability
profiles from fiber configurations whose fiber orientations are specified by the azimuthal angles
01 300, 2 = {200, 100} and 03 = {200, 550, 100} for the 1, 2 and 3 fiber systems respectively.
Polar angles for all fibers were taken to be 900 so that a view from the +zaxis will clearly depict
the individual fiber orientations. Computations with other polar angles yielded similar quality
results. In all computations, the Laplace series were terminated after 1 = 8.
Figure 3 shows the effect of varying Ro on the constructed probability surfaces. Increasing
Ro gives rise to the sharpening of the displacement PDFs. This could be predicted from Fig.
2a that indicates that for small Ro the largest contribution comes from Io, which upon the
spherical harmonic transform forms the isotropic part of the constructed probabilities. When
Ro is greater than the radius of the cylinder confining the water molecules and the characteristic
length /6Dt associated with the diffusion process (which is 15/m for the , . i in Figure 2a),
the distribution of probability on the surface becomes sharper and individual fiber populations
are better resolved.
We also computed the probability surfaces for a simulated image of fiber crossings shown in
Figure 4. The surfaces are consistent with the underlying known fibrous structure. The circular
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
resolved when the different fiber bundles make larger angles with each other. Fig. 4b shows
the probability profiles when there is no noise added to the signal values. Similar to Jansons
arr r ,
Figure 4: (a) Simulated , . i, of two crossing fiber bundles. (b) Probability surfaces computed
using the expansion of the probability on the surface of a sphere. (cf) Surfaces in the framed
area of panel b recomputed under increasing levels of noise added to the signal values. These
panels represent images with signaltonoise ratios (SNRs) between 50:1 and 12.5:1 in the non
diffusionweighted image.
and Alexander (2003), we have added Gaussian noise of increasing standard deviation to the
real and complex parts of the signal. When the signal intensity in the image with no diffusion
weighting is taken to be centered around 1, and Gaussian noise of standard deviations 0.02
through 0.08 is added (in equal steps), the probability profiles shown in Fig. 4cf are obtained.
These panels represent images with signaltonoise ratios (SNRs) between 50:1 and 12.5:1 in
the nondiffusionweighted (So) image. Note that in our standard HARDI protocol, we obtain
SNR values in excess of 30 in diffusionweighted scans and about 100 in nondiffusionweighted
images. Therefore, in real experiments, one can expect to achieve results that will be of similar
or better quality with the image presented in Fig. 4c.
To provide a more quantitative assessment of the DOT method and its sensitivity to increas
ing noise levels, we took the HARDI simulations of 1, 2 and 3 fiber profiles presented in Figure
1 and numerically computed the fiber orientations by finding the maxima of the probability
profiles (see Table 2). In this table, 0 denotes the angle between the computed and the true
fiber orientations in degrees whereas a is the noise level. Note that when no noise was intro
duced (a = 0), there was a small deviation of the computed fiber direction from the true fiber
orientation because of the finite sampling of the hemisphere (at 81 gradient orientations), the
termination of the LS at order 8 and the precision of the numerical procedure used to compute
the maxima of the probability profiles. The simulations of the signal profiles with noise were
"I IN
IN
N N
Ilk IN IN
IN 14 *N
Table 2: The angle between the computed and true fiber orientations (deviation angles) in
degrees. Second column presents the deviation angle of each fiber when the DOT of noiseless
signal profile is taken. Columns 36 show the mean and standard deviation values for the
deviation angle when Gaussian noise of standard deviation 0.02 to 0.08 (from left to right) was
added to the signal profiles. The computations for the DOT of noisy signals were repeated 100
times.
S(a 0) (a = 0.02) 0 (a 0.04) 0 (a 0.06) 0 (a = 0.08)
1 fiber {0.364} 0.77 0.42 1.44 0.79 2.20 1.09 3.08 1.66
2 fibers {1.43, 0.80} 2.33 1.10 3.66 2.01 6.00 5.57 8.07 7.92
3 fibers {2.87, 0.60, 4.57} 5.81 5.84 11.5 10.1 14.7 10.3 17.6 11.9
repeated 100 times for each noise level to provide a distribution of deviation angles. We report
the mean and standard deviations of these distributions in columns 36 of Table 2. As expected,
the 0 values increase with increasing noise and it is more challenging to accurately resolve the
distinct fiber orientations when there are more fiber orientations.
4 Imaging Parameters
To test the performance of the DOT, we calculated the orientation probabilities on HARDI
data from three anatomical regions of excised, perfusionfixed rat nervous tissue (optic chiasm,
brain and spinal cord). These experiments were performed with the approval of the University
of Florida Institutional Animal Care and Use Committee. The images were acquired at 17.6T
(brain) or 14.1T (spinal cord and optic chiasm) using Bruker Avance imaging systems. A
diffusionweighted spin echo pulse sequence was used. Diffusionweighted images were acquired
along 81 (brain) or 46 (spinal cord and optic chiasm) directions with a bvalue of 1500s/mm2
(brain and spinal cord) or 1250s/mm2 (optic chiasm) along with a single image acquired at
b Os/mm2. Echo times were 23, 28, 25ms ; repetition times were 0.5, 2, 1.17s ; A values
were 12.4, 17.8, 17.5ms ; 6 values were 1.2, 2.2, 1.5ms ; bandwidth was set to 35, 32, 39kHz
; signal averages were 10, 6, 7 ; matrix sizes were 128 x 128 x 5, 100 x 100 x 60, 72 x 72 x 40
and resolutions were 33.6 x 33.6 x 200pm3, 150 x 150 x 300pm3, 60 x 60 x 300pm3 for optic
chiasm, brain and spinal cord data respectively. The optic chiasm images were signal averaged
to 67.2 x 67.2 x 200/m3 resolution prior to probability calculations. In Figure 5 we show a
particular axial slice from a HARDI dataset collected from excised rat spinal cord.
5 Scalar Indices
Many clinical studies employ scalar rotationally invariant measures derived from diffusion MRI
data to quantify the changes occurring with development and pathologies. Recently, Ozarslan
Figure 5: Representative HARDI data set from an excised, perfusionfixed rat spinal cord. At
the upper left corner is the image with no diffusion weighting followed by 45 diffusionweighted
images.
et al. (2005) demonstrated that generalized models more accurately quantify anisotropy measures
compared to DTI. In this section, we discuss the estimation of the generalized scalar indices
from the probability values evaluated using DOT, and demonstrate the images constructed by
computing these measures on a pixelbypixel basis. For completeness, formulation of anisotropy
in terms of both variance and entropy are provided.
5.1 Anisotropy from Variance
In Ozarslan et al. (2005) generalized anisotropy indices based on the variance of the values
of an arbitrary integrable positivedefinite function defined on the unit sphere were presented.
When applied to functions represented in terms of spherical harmonics, like the parametrically
reconstructed P(Ror) in this work, the variance takes a particularly simple form given by
00 I
V " 'I 2 (23)
Sl=2 m=l
Then, using the scaling relationship provided in Ozarslan et al. (2005), it is possible to map
the values of the variance to the interval [0, 1) where the resulting index is called generalized
anisotropy (GA).
5.2 Anisotropy from Entropy
The function defined on the sphere can be taken as a PDF simply by normalizing the integral of
the probability profile over the sphere via a multiplication of the /,,. coefficients by 1//47roo.
Figure 6: I1(u) images up to I = 10 from coronal sections of rat brain when u is chosen to point
through the image plane.
Then, it is meaningful to define the entropy associated with this distribution. By using the
expression in Ozarslan et al. (2005) for the entropy of a general function on the unit sphere, it
is possible to show that the entropy in our case is given byl
In (47r(P(Ror)) ) i P I f du P(Ror) lnP(Ror)
47(P(Ror))
S In p0 (24)
where Aim are given by the SHT of In P(Ror). Similar to the transformation of variance values
into the GA index, the entropy values can be transformed into an anisotropy index that was
called scaled entropy (SE) (Ozarslan et al., 2005).
6 Results
The probability maps were calculated by following the procedure described in the theory section.
Terms up to I = 8 were included in all calculations. Representative images of the Iz(u) values,
when u is chosen to point through the image plane, are presented in Figure 6. Note that
the intensity values in the Iio(u) image are very small. The /'.. coefficients that generate the
probability surfaces are shown in Fig. 7 for the same slice. It was not necessary to show the
coefficients with negative m values because of Eq. 14. Note that this relationship also ensures
that pio are real.
The computed /,, components were used in the calculation of the scalar measures described
in the previous section. In Figure 8, we show the variance and GA maps computed from the
optic chiasm and brain datasets. The GA index was calculated both from the ,.. coefficients
1The expression for the entropy is a slightly modified version of the expression in Ozarslan et al. (2005), where
the 3 in the argument of the natural logarithm is replaced by 47r. The reason for this modification stems from the
difference between the normalization conditions imposed on the functions DN(u) in Ozarslan et al. (2005) and
P(Ror) in this work.
lx10I
POD
IRe{p,2,}
IRm{pm,.
IRe{p4m}
m=0
Figure 7: The LS coefficients
rat brain data set.
1x105
L
5x106
to
m=1 m=2 m=3 m=4
of the probability profile up to I = 4 computed from the excised
and a secondorder tensor fit to the data. It is apparent that although the GA values are similar
in the unidirectional section of the optic chiasm, in the region of decussating optic nerve fibers,
GA values implied by DTI were lower than those calculated from the probability surfaces. Also
included are the entropy (a) and the SE maps calculated from both samples. It should be noted
that V, GA, a and SE values depend on the choice of Ro.
GA(DTI)
Figure 8: From left to right: Variance, GA, GA from the rank2 tensor, entropy and SE images
calculated for excised rat optic chiasm (top) and brain (bottom) samples.
L 
GAP
mamxll) maxllXl) m
optic nerve
optic chiasm
optic tract
Figure 9: Probability maps calculated from a rat optic chiasm data set overlaid on an axially
oriented GA map. The decussations of myelinated axons from the two optic nerves at the center
of the optic chiasm are readily apparent using the DOT method. These crossing fibers carry
information from the temporal visual fields to the contralateral hemispheres. In figures 9 through
12, the orientation surfaces are colorcoded such that portions of the surfaces pointing towards
or away from the reader are green and blue respectively (see inset).
Visualization of the probability profiles was done by computing the probabilities along many
directions and displaying the parametrized surface defined in Eq. 22. To increase the sharpness
of the pr'! 1 .1.ilir, profiles, we have subtracted the minimum probability existing in the profile
from all probability values. This process was followed by a normalization of the surface to fill
the cube (voxel) it will be located in so that all visualized surfaces have similar sizes. As a result,
the visualized surfaces emphasize the dire l i.11 I ,il ', but are not intended to provide information
on the true values for the probabilities.
We have overlaid these orientation surfaces on generalized anisotropy (GA) maps (Ozarslan
et al., 2005) computed from the displacement probabilities shown above. The coloring schemes
proposed for DTI orientation visualization (Pajevic and Pierpaoli, 1999) are not readily appli
cable to probability surfaces. The i" .i:' ,.1i  on the image plane is obvious. However, one
may miss the orientations through the image plane. To prevent this, we colorcoded the surfaces
such that as the values of the zcomponent of the parametrized surface vary from the maximum
probability value present in the probability profile to minus this maximum pr'1 1..iliilv value,
parietal
cortex
I corpus
external capsule callosum
~ stratum
stratum f radiatum
lacunosummoleculare
hilum_ Lmolecular
\ layer
Figure 10: P .1. d.1ili , map of a coronallyoriented GA image of the rat brain. Diffusion fiber
orientations in the parietal cortex were collinear with the apical dendrites and axons of cortical
pyramidal neurons found in cortical layers IIIV. In the dorsal hippocampus, the molecular
layer and stratum radiatum fiber orientations paralleled the apical dendrites of granule cells
and pyramidal neurons respectively, whereas in the stratum lacunosummoleculare orientations
paralleled Schaffer collaterals from CA1 neurons and perforant fibers from the entorrhinal cortex.
the color of the surface changes from green to blue. In all calculations Ro was set to 16/m, and
the last term kept in the Laplace series was 1 = 8.
The rat optic chiasm is a distinct white matter structure with both parallel and decussating
optic nerve fibers, thus providing an excellent experimental validation for our approach. Figure
9 shows the displacement probabilities computed from the optic chiasm image where every other
pixel was included for the sake of clarity. This figure demonstrates the distinct fiber orientations
in the central region of the optic chiasm where myelinated axons from the two optic nerves cross
one another to reach their respective contralateral optic tracts.
Figure 10 shows the displacement probabilities calculated from excised coronal rat brain
MRI data. At the top left is a diffusionweighted image that shows the selected ROI. This
region is expanded in the large image and depicts the orientations of the highly anisotropic
and coherent fibers of the external capsule and corpuscallosum bordered inferiorly by the hip
pocampus and superiorly by radial cortical tI ,i. I 1. i, . Note, voxels with crossing orientations
located superiorly to the external capsule represent the interdigitation of fibers entering the
cortex from the external capsule and the corpus callosum into the radial orientations of the cor
Sr, on of th r cuneate nucleus
St ie o and fasciculus
spinal trigeminal tra~l
I hyperglosal nucleus
spinal trigeminal nucleus medial longitudinal
fasciculus
olivocerebellar and dorsal
spinocerebellar tracts reticular nucleus
medial lemniscus
ventral
spinocerebellar tracts
pyramidal tract
Figure 11: The Diffusion Orientation Transform (DOT) described in this paper also character
ized the complex cytoarchitecture of the rat brainstem well as shown in this probability map
from one side of the rat medulla.
tex. Future investigations employing this method should improve our understanding of normal
and pathologicallyaltered Inti lci...ii ,_ in regions of complex fiber architecture such as the rat
brainstem (Figure 11).
Finally, we show the probability maps computed from excised spinal cord data in Figure 12.
Again the ROI is specified on a diffusionweighted image shown on the top left section. The
corresponding orientation maps are depicted in the top right panel. Selected pixels of this image
were enlarged on the bottom panels of the figure. To demonstrate the shapes more clearly, seven
selected surfaces were rotated by 900 about the xaxis so that the upanddown direction in
the individual surfaces shown in blue correspond to inandout direction in red images. The
magnified surfaces may represent locations where ventral root fibers from amotor neurons cross
white matter to enter the gray matter of the spinal cord.
7 Discussion
7.1 Exponential Attenuation Assumption
We have assumed so far that the signal attenuation along each radial line in qspace is char
acterized by an monoexponential decay. Therefore, it was possible to extract orientational
information from data acquired on a single spherical shell and at the origin of the qspace. We
would like to note that this is the very assumption that is intrinsic to DTI, establishing the
correspondence between the diffusion tensor and the assumed Gaussian PDF whose orienta
X
X X
Figure 12: Pi .1.il1i, maps calculated from a diffusionweighted dataset acquired from an
excised rat spinal cord. The surfaces in the bottom row depict the probability profiles selected
from the image matrix and rotated 900 about the xaxis. These crossing fiber orientations may
represent coherent ventral root spinal nerve fibers penetrating along the xaxis perpendicular
to ascending and descending white matter axons in the anterior funiculus to reach the anterior
horn motor neurons.
tional mode is estimated from the principal eigenvector of the diffusion tensor. The satisfactory
performance of DTI in systems with single fiber orientations has prompted us to keep the mo
noexponentiality assumption for the radial behavior while complicating the angular structure;
this results in nonGaussian pr 1 .1 iilli, profiles. This assumption worked with our simulations
and with real datasets.
In Figure 13, we show simulated signal values from a one fiber , I. i1 with a fiber radius of
10pm when the angle between the fiber orientation and the gradient direction varies between
0 and 900 in steps of 300. Also shown are the signal values assumed by the exponentiality as
sumption when the experiment is performed at b = 1500 s/mm2. Note that this is a logarithmic
plot, therefore the true deviations between the real and assumed signal values are much smaller
than what they appear on the right side of the plot when the signal values are small. Because
the qspace is the frequency space for the probabilities, from a signal processing pointofview
the exponentiality assumption can be thought of as a lowpass filtering of the true probability
values. Therefore, the result is a broadened PDF. A correction scheme for this effect is described
.+ 0 ........ 30
0:'x  300
01 60,
10 1 k x 90
S \ .
0 ,. x X X X
103 A
+' A ,
10 4 + +
105
0 2000 4000 6000 8000 10000
b (s/mm2)
Figure 13: Logarithmic plot of the signal attenuation values as a function of the bvalue. The
I, .. .1, indicate the signal attenuations calculated from a simulation of a single fiber , ., i_
while the lines indicate the monoexponential fits when the HARDI experiment is performed at a
bvalue of 1500s/mm2. The curves correspond to different angles between the diffusion gradient
and fiber orientations.
below.
We also have investigated the effect of bvalue on the constructed probability surfaces. To
this end, we simulated HARDI experiments performed at increasing bvalues on 1 and 2 fiber
systems. We also repeated the simulations for fibers of radii 5 and 10 pm. The selected results
are provided in Figure 14. The most reassuring finding is that there has been no realizable
alteration in the peaks of the distributions indicating that the calculated fiber orientations are
robust to the choice of bvalue. However, it is evident that the probability surfaces are sharper
and multiple orientations are better resolved at higher bvalues. As we have demonstrated before,
a bvalue of 1500 s/mm2 seems sufficient to resolve the fiber crossings when the radii of the fibers
are 5 pm. However, when the radii are doubled, it is advantageous to collect the data at higher
bvalues. It should be noted that spurious peaks start to develop at very high bvalues (see the
first two rows of the last column). This may be explained by the crossing of the signal decay
curves in Figure 13, which :_:_ I that at high bvalues the order of signal values from different
orientations may be altered.
7.2 Extension to Multiexponential Attenuation
We have thus far employed the monoexponentiality assumption of the signal attenuation. How
ever, the same formalism provides a surprisingly simple extension to multiexponential attenua
bvalue (s/mm2)
500 1500 3000 5000 8000
5 ym 4r oArO 0 O
10m m 4
Figure 14: Simulations of 1 and 2 fiber systems as a function of bvalue where the radii of the
fibers were taken to be 5 and 10 pm.
tion2, which has been shown in numerous articles to provide a very accurate characterization of
the radial behavior (in qspace) of the MR data collected from tissue (e.g. Niendorf et al., 1996).
To derive the correct generalization, we start by replacing the S' i1: 1iTanner equation (Eq.
2) with the expression:
S(b, u) NE
s(b f(u) bD(u, (25)
where NE is the number of terms exponentialss, transients) in the series, Di(u) is the ith
diffusion coefficient for the gradient direction u, and fi(u) is the \..!iiin . fraction" of the ith
exponential I lif, in:_ the relationship
NE
ft(u)= 1. (26)
Carrying out the same algebra as before, Eqs. 5 and 17 hold with the definition
NE
1(u) = f,(u) 1(u), (27)
i= 1
2We would like to stress that we utilize the multiexponential fit solely to provide an approximation and extrap
olation to the signal attenuation and by no means do we intend to make inferences about the compartmentation
in tissue from this fit.
where
Ili(u) 47/ dqq2 jl(27qRo) exp(472q2tDi(u)) (28)
JoO
which is the same expression when D(u) in Eq. 6 is replaced by Di(u). Therefore, either of the
forms given in Eqs. 7 or 8 can still be used to calculate Ii(u) from Di(u).
The extension to multiexponential attenuation requires the following modifications to the
implementation of the DOT:
1. Fit multiexponential function (Eq. 25) along each radial line in qspace to estimate fi(u)
and Di(u).
2. For each diffusion coefficient Di(u) corresponding to each term in the series, calculate
Ii(u) from Eq. 28.
3. Calculate I(u) from Eq. 27.
4. Apply either the parametric or nonparametric reconstruction as before.
Figure 15 shows the biexponential fits to the data points already presented in Figure 13.
The improvement in the functional fits is evident.
We have tested the proposed extension scheme on our simulated data from 1, 2 and 3
fiber systems. The results are shown in Figure 16. It is clear that the monoexponential and
multiexponential fits provide the same orientational information, yet the constructed probability
100
+  0
0 ........ 300
,10 1 A  600
101 nx ..o
103 .
104 90
xx%
b (s rrn )
Figure 15: The ,i.1 indicate the signal attenuations calculated for the same ,i as
in Figure 13 whereas the lines indicate the curves obtained from a biexponential fit to these
data points. The curves correspond to different angles between the diffusion gradient and fiber
orientations.
monoexponential
(b=2500s/mm2) biexponential triexponential
(b=2500s/mm )
Figure 16: Simulations of 1, 2 and 3 fiber systems with mono (b 2500s/mm2), bi and
triexponential fits (from data up to b 9000s/mm2). Similar to Figure 1, the orientations
of the simulated fibers are specified by the azimuthal angles i1 = 300, 2 = {200, 100} and
3 = {200, 550, 100} for the 1, 2 and 3 fiber systems respectively. All fibers lie on the image
plane.
surfaces in the latter case resolve the distinct fiber orientations better, most notably in the
3 fiber system. However, the results indicate that transition from bi to triexponential fits
does not result in a significant improvement. This demonstrates the sufficient accuracy of the
biexponential fits to the signal attenuation values.
Unfortunately, using a biexponential attenuation fit in our formalism would necessitate col
lecting about three times the number of data points when compared with the case in which the
monoexponentiality assumption is made. This is because there are 2 NE unknowns in the fit,
and if one choses to collect data at b = 0, then 2 NE 1 spherical shells have to be sampled
for the NEexponential fits.
8 Conclusion
The DOT technique provides a direct estimation of displacement probability surfaces within
each voxel from multiorientational diffusionweighted MRI data. The method is robust and
fast. DOT can be implemented nonparametrically for direct estimation of probability values
along desired directions, or by using an SHT that gives the Laplace series coefficients of the
probability profile. In either case, high resolution probability surfaces can be reconstructed
easily from the signal values. Further, the behavior of the MR signal intensities with increasing
bvalues can be characterized by mono or multiexponential fits. Our findings indicate that
multiexponential fits result in higher quality reconstructions. However, when the acquisition
time or the available gradient strength is limited, the monoexponentiality assumption can be
employed that results in some broadening of the PDF whose angular structure is smoother. As
demonstrated in excised rat nervous tissue, the potential applications of our approach include
more accurate estimates of fiber orientations that will improve the existing fiber tractography
schemes and enable the reliable mapping of more neural connections between different parts of
fibrous tissues.
Acknowledgments
All MRI data was obtained at the Advanced Magnetic Resonance Imaging and Spectroscopy
(AMRIS) f 1 ili, in the McKnight Brain Institute at the University of Florida. This research
was supported by the National Institutes of Health Grants R01NS42075, R01NS :1,.','2 and
P41RR16105, and the National High Magnetic Field Laboratory (NHMFL), Tallahassee.
Appendix A. Confluent Hypergeometric Functions of the First
Kind
The confluent hypergeometric function of the first kind (also known as Kummer's function of
the first kind or Kummer's function) 1Fi(a; b; x) is given by the series (Abramowitz and Stegun,
1977)
IFi(a; b; x) (0 (29)
S( (b)k k! (29)
k=0
where (a)k = a(a + 1)(a + 2) ... (a + k 1) with (a)o = 1.
Among others, the confluent hypergeometric function of the first kind satisfies the recurrence
relations
(ba)iFi(a 1;b;)+(2a b +)Fi(a;b;) aiFi(a +l;b;) = 0
b(b1)iFi(a;b1;x)+b(1bx)iFi(a;b;x)+x(ba)iFi(a;b+1;x) = 0
(1+ab)iFi(a;b;x)aiFi(a+1;b;x)+(b1)iFi(a;b1;x) = 0.
(30)
Many of the commonly used functions are special instances of the confluent hypergeometric
function of the first kind. For example,
IFi(a;a;x) ex (31)
F1 ; 2) 2 erf(x) (32)
(2' 2 ) 2x
Finally, the asymptotic behavior of the function iFi(a; b; x) as Ix x oo when x is real, is
given by
1Fl(a;b;x) eCxa (a),(l+ab), ]
F(b) F(b a) n!
+ Zab S (b a) (1 a) (33)
F(a) + o( (33)
Appendix B. At and B1 Coefficients
Although the recursion relations of the confluent hypergeometric functions provided in Eqs. 30
are useful in seeing that the Ii(u) functions can be expressed as the sum of two terms (one
involving exponential and the other involving error functions as in Eq. 8), the derivation of the
analytical form of the Al and Bi coefficients using these recursion relations is a formidable task.
Therefore, in order to find analytical expressions for the Al and Bi coefficients for a general 1
value, we make a termbyterm comparison of the asymptotic form of the 4I function evaluated
from Eq. 33 with the asymptotic form of Eq. 8. After some tedious algebra we have found that,
if Aln and Bin are defined such that
1/2 1/21
Al(u) = Al/3(u) 2 and Bl(u) = Bin/3(u)2 (34)
n=0 n=0
the following expressions hold:
SA ,if n < 2
Am= A o ynl (1)t (2t3)!!( )n 1 )! if > 2 (35)
St=i (1/2)(ntI)! 21/22n +t
where (1 + 1)!! (I + 1)(1 1)... 1 and
A (n1)!/2" an (36)
and
S(1 ). 1 +1)!!
B1 22 (37)
lr F ( 2 ) n! 21/212n
We have verified using Mathematica that these expressions indeed yield the correct coeffi
cients for the Ii(u) functions.
Appendix C. Convergence of the Laplace Series for the Probabil
ity Profile
Theorem: The series given by (see Eq. 5)
l I
P(Ror) E (i) Ytm(r) / du Yl(u)* I1(u) (38)
1=0 m=l
is convergent.
Proof: We start by inserting the upper bound for the spherical Bessel functions of order 1
(Abramowitz and Stegun, 1977)
ir(irqRo)' (2irqRo)'
Ij(2qRo) < (39(q) (27qRo)' (39)
jl(2rq ) 2(1 + ) (21 + 1)!!(39)
into Eq. 6. This yields the upper bound for the functions 11(u)
(1 + 1)!! R2 1/2 (40)
( (21 + 1)!!21/2(47rD(u)t)3/2 \D(u)tJ
Note that using the addition theorem for spherical harmonics (Arfken and Weber, 2001),
21+ 1
E m(r)Y m(u)* 21 Pl(u. r), (41)
47
m=l
where Pl(x) is the lth order Legendre polynomial, it is possible to express Eq. 5 in the following
form:
P(Ror) = 9 (42)
l=0
where
S i)(21 1) dul1(u)Pl(r.u). (43)
47
Using the generating function for the Legendre polynomials, it is possible to prove that (Arfken
and Weber, 2001)
P (cos)l < 1. (44)
Using Eqs. 40 and 44, it is easy to see that
(1 + 1)!! R2 1/2
S(21 1)!! 2/2 (4Dmint)3/2 Dmint (45)
where Dmin = min D(u). Note that
lim 1+2 lim 1+3 R =0 (46)
1>o0 (I 100 (21 + 1)(21 + 3) 2Dmint
Therefore, using the d'Alembert (Cauchy) ratio test, the series Zo0 l converges. Using the
comparison test, it is straightforward to see that the series Z0o I cll converges since 0 < QiI < 1.
Therefore, the series in Eq. 38 is absolutely convergent.
27
References
Abramowitz, M., Stegun, I. A., 1977. Handbook of Mathematical Functions: With Formulas,
Graphs, and Mathmetical Tables. Dover Publications, New York.
Alexander, D. C., Barker, G. J., Arridge, S. R., 2002. Detection and modeling of nonGaussian
apparent diffusion coefficient profiles in human brain data. Magn Reson Med 48 (2), 331 340.
Arfken, G. B., Weber, H. J., 2001. Mathematical Methods for Physicists. Academic Press, San
Diego.
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), ,1',978.
Basser, P. J., 2002. Relationships between diffusion tensor and qspace MRI. Magn Reson Med
47, 392397.
Basser, P. J., Mattiello, J., LeBihan, D., 1994a. Estimation of the effective selfdiffusion tensor
from the NMR spin echo. J Magn Reson B 103 (3), 247254.
Basser, P. J., Mattiello, J., LeBihan, D., 1994b. MR diffusion tensor spectroscopy and imaging.
Bi., .i. J 66 (1), 259267.
Basser, P. J., Pajevic, S., Pierpaoli, C., Duda, J., Aldroubi, A., 2000. In vivo fiber tractography
using DTMRI data. Magn Reson Med 44, 625632.
Callaghan, P. T., 1991. Principles of Nuclear Magnetic Resonance Microscopy. Clarendon Press,
Oxford.
C'I !:'. i.t, T. L., Bin l,!l:_. J. A., Pipe, J. G., 1990. Anisotropic diffusion in human white
matter: demonstration with MR techniques in vivo. Radiology 177 (2), 328329.
C'!, I., !.11, G. G., C'l! .!:_. D. C., Hazlewood, C. F., Rorschach, H. E., 1976. Nuclear magnetic
resonance measurement of skeletal muscle: anisotropy of the diffusion coefficient of the intra
cellular water. Biophys J 16 (9), 1043 1053.
Conturo, T. E., Lori, N. F., Cull, T. S., Akbudak, E., Snyder, A. Z., Shimony, J. S., McKinstry,
R. C., Burton, H., Raichle, M. E., 1999. Tracking neuronal fiber pathways in the living human
brain. Proc Natl Acad Sci 96, 10422 10427.
Frank, L. R., 2002. Characterization of anisotropy in high angular resolution diffusionweighted
MRI. Magn Reson Med 47 (6), 10831099.
Inglis, B. A., Bossart, E. L., Buckley, D. L., Wirth, E. D., Mareci, T. H., 2001. Visualization
of neural tissue water compartments using biexponential diffusion tensor MRI. Magn Reson
Med 45 (4), 580587.
Jansons, K. M., Alexander, D. C., 2003. Persistent angular structure: new insights from diffusion
magnetic resonance imaging data. Inverse Problems 19, 1031 1046.
Liu, C. L., Bammer, R., Moseley, M. E., 2003. Generalized diffusion tensor imaging (GDTI): A
method for characterizing and imaging diffusion anisotropy caused by nonGaussian diffusion.
Isr J Chem 43 (1 2), 145 154.
Maier, S. E., Vrli I. , .ii. S., Mamata, H., Westin, C. F., Jolesz, F. A., Mulkern, R. V., 2004.
Biexponential diffusion tensor .111 i, i, of human brain diffusion data. Magn Reson Med 51 (2),
321 330.
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, 265269.
Moseley, M. E., Cohen, Y., Kucharczyk, J., Mintorovitch, J., Asgari, H. S., Wendland, M. F.,
Tsuruda, J., Norman, D., 1990. Diffusionweighted MR imaging of anisotropic water diffusion
in cat central nervous system. Radiology 176 (2), 439445.
Niendorf, T., Dijkhuizen, R. M., Norris, D. G., van Lookeren Campagne, M., Nicolay, K., 1996.
Biexponential diffusion attenuation in various states of brain tissue: implications for diffusion
weighted imaging. Magn Reson Med 36 (6), 847."7.
Ozarslan, E., Mareci, T. H., 2003. Generalized diffusion tensor imaging and analytical relation
ships between diffusion tensor imaging and high angular resolution diffusion imaging. Magn
Reson Med 50, 955 ', 1'
Ozarslan, E., Vemuri, B. C., Mareci, T. H., 2004a. Generalized diffusion tensor imaging of
excised rat brain. In: Proc. of the 12th Scientific Meeting of ISill;\!. p. 92.
Ozarslan, E., Vemuri, B. C., Mareci, T. H., 2004b. Multiple fiber orientations resolved by gen
eralized diffusion tensor imaging. In: Proc. of the 12th Scientific Meeting of ISi I\!. p. 89.
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), ,N1,.876.
P,,i,. 1ii S., Pierpaoli, C., 1999. Color schemes to represent the orientation of anisotropic tissues
from diffusion tensor data: application to white matter fiber tract mapping in the human
brain. Magn Reson Med 42 (3), 526540.
Parker, G. J. M., Alexander, D. C., 2003. Probabilistic Monte Carlo based mapping of cerebral
connections utilising wholebrain crossing fibre information. In: Taylor, C. J., Noble, A. (Eds.),
Information Processing in Medical Imaging (IPMI'03). Vol. 2737 of Lecture Notes in Computer
Science. Springer, pp. 242254.
Ritchie, D. W., Kemp, G. J. L., 1999. Fast computation, rotation, and comparison of low
resolution spherical harmonic molecular surfaces. J Comput Chem 20 (4), 383395.
Schwabl, F., 1989. Quantum Mechanics. SpringerV. i1 :_. Berlin.
Siderman, O., J6nsson, B., 1995. Restricted diffusion in cylindirical geometry. J Magn Reson
A (117), 9497.
Sj ii. i1.. E. O., 1 1".' Use of spin echoes in a pulsed magneticfield gradient to study anisotropic,
restricted diffusion and flow. J Chem Phys 43 (10), 35973603.
Sj' ii. 1. E. O., Tanner, J. E., 1i.". Spin diffusion measurements: Spin echoes in the presence
of a timedependent field gradient. J Chem Phys 42 (1), 288292.
Tournier, J. D., Calamante, F., Gadian, D. G., Connelly, A., 2004. Direct estimation of the fiber
orientation density function from diffusionweighted MRI data using spherical deconvolution.
Neurolmage 23, 11761185.
Tuch, D. S., 2004. Qball imaging. Magn Reson Med 52, 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., 2003. Diffusion MRI of complex neural
architecture. Neuron 40, '.895.
Tuch, D. S., Weisskoff, R. M., Belliveau, J. W., Wedeen, V. J., 1999. High angular resolu
tion diffusion imaging of the human brain. In: Proc. of the 7th Annual Meeting of IS\lI;\!,
Philadelphia. p. 321.
von dem Hagen, E. A. H., Henkelman, R. M., 2002. Orientational diffusion reflects fiber structure
within a voxel. Magn Reson Med 48 (3), 454459.
Wedeen, V. J., Reese, T. G., Tuch, D. S., Weigel, M. R., Dou, J. G., Weiskoff, R. M., C'i. I. i,
D., 2000. Mapping fiber orientation spectra in cerebral white matter with Fourier transform
diffusion MRI. In: Proc. of the 8th Annual Meeting of ISMRM, Denver. p. 82.
