Neuronal Fiber Tracking in DTMRI1
T. E. McGraw, B. C. Vemuria Y. Chen, M. Rao b T. Mareci
aDept. of Computer and Information Sciences and i.,'l.' ring
University of Florida, Gainesville FL, 32611
bDept. of Mathematics, Uil,, I, 't, of Florida, Gainesville FL, 32611
CDept. of Biochemistry, Uil.,, I ~,t, of Florida, Gainesville FL, 32610
Abstract
Diffusion tensor imaging (DTI) can provide the fundamental information required
for viewing structural connectivity. However, robust and accurate acquisition and
processing algorithms are needed to accurately map the nerve connectivity. In this
paper, we present a novel algorithm for extracting and visualizing the fiber tracts in
the ('XS specifically the spinal cord and brain. The automatic fiber tract mapping
problem will be solved in two phases, namely a data smoothing phase and a fiber
tract mapping phase. In the former, smoothing is achieved via a weighted TVnorm
minimization, which strives to smooth while retaining all relevant detail. For the
fiber tract mapping, a smooth 3D vector field indicating the dominant anisotropic
direction at each spatial location is computed from the smoothed data. Neuronal
fibers are then traced by calculating the integral curves of this vector field.
Results are expressed using three modes of visualization. (1) Line integral convo
lution (LIC) produces an oriented texture which shows fiber pathways in a planar
slice of the data. (2) A streamtube map is generated to present a threedimensional
view of fiber tracts. Additional information, such as degree of anisotropy, can be
encoded in the tube radius, or by using color. (3) A particle system form of visual
ization is also presented. This mode of 1ibpiv allows for interactive exploration of
fiber connectivity with no additional preprocessing.
Key words: DTI, DTMRI, Diffusion Tensor, Total Variation
PACS:
Email address: {tmcgraw, vemuri}@cise.ufl.edu (T. E. McGraw, B. C.
Vemuri).
1 This research was funded in part by NIH grant RO1NS 121i7",.
Preprint submitted to Elsevier Science
10 D~ecemnber 2002
1 Introduction
The understanding of neurological function, and medical diagnosis of disease
and injury require knowledge of the structural organization of the brain. Mag
netic Resonance Imaging provides a noninvasive means of studying anatomi
cal connectivity. This connectivity, in the form of axonal nerve fiber bundles,
can only be measured indirectly. The presence of these fibers must be inferred
from the behavior of water molecules near the tissue.
Diffusion tensor magnetic resonance imaging (DT '.l I) is a method of mea
suring the rate of water diffusion in biological structures. A diffusion tensor
at each location on a regular lattice describes a volumetric average of the
directional properties of diffusion within each voxel. The observation that dif
fusion is anisotropic in areas of whitematter fiber bundles allows tracking of
fibers through the lattice. However, acquisition noise corrupts the data mea
surements. Voltage variations in the receiving coil of the '.11 I machine due to
thermal noise are a major source of signal degradation. The noise is assumed
additive and is modeled by a zeromean Gaussian distribution. Extracting
fiber paths from a noisy field is an illposed problem. The goals of the research
reported here are to extract fibers from noisy DT'.II1l data and to visualize
the same using techniques which present the fibers in the most biologically
informative way possible.
1.1 Overview
The fiber tracts are estimated in a two stage process, (a) smoothing of the raw
directional images acquired for varying magnetic field strengths and estimating
the diffusion tensor, D, field over the 3D image lattice, (b) computing the
dominant eigenvector field from this regularized D and estimating regularized
streamlines/integral curves as the desired fiber tracts.
Computed fibers must be visualized in a way that is useful for anatomical study
or medical diagnosis. Existing vector field visualization techniques, such as
streamlines/streamtubes can be used to visualize fibers. These techniques can
also be modified to convey more information about diffusion by incorporating
quantities derived from the tensor. Here, we adapt line integral convolution
(LIC), streamtubes, and particles to the task of tensor field visualization.
We will now review the literature on the raw data smoothing, streamline
generation, and visualization techniques where appropriate.
1.2 Contributions
Although it seems natural for one to perform denoising of DT'.i I1 data prior
to computing any meaningful quantitative measures (e.g., dominant diffusion
directions, etc), most of the reported research in literature to date does not
focus on this important issue. Instead, the diffusion tensor is computed di
rectly from the raw noisy measurements. Dominant direction and magnitude
of diffusion (Eigen vectors and values) are computed from this noisy D and
then smoothed. It should be noted that Eigen values and vectors computed
from a noisy estimated tensor may not be meaningful (unless the data is al
ways assumed to be of a high signaltonoise ratio.) We propose to smooth
the data before these calculations to provide for better accuracy in estimated
nerve fiber maps.
A new application of particlebased visualization is presented. This technique
has previously been used in fluid mechanics for visualization of velocity fields.
Here we will adapt the technique to tensor field visualization. This technique
has no associated preprocessing time, and allows realtime interactive visual
ization of the nerve fibers. '1 1i1, other techniques in literature require extensive
computing time for streamline integration prior to visualization.
The rest of this paper is organized as follows. In Section 2 we will discuss the
diffusion process and give a short introduction to the DT'.II I data acqui
sition process. An overview of PDE based image denoising will follow. Some
linear and nonlinear techniques for scalar and vectorvalued images will be
reviewed. Variational techniques for formulating image denoising, and the TV
norm will be presented. Section 3 contains a description of our vectorvalued
image restoration technique, a weighted TV norm minimization. The noise
model used to formulate this technique is discussed, as well as the motivation
for modifying the TV norm usually used for vectorvalued images. In Section
4 we present the process of tracking neuronal fiber bundles through the re
stored DT'.II I data. Previous methods of fiber tracking will be discussed.
A quantification of diffusion tensor anisotropy will be given. In this section
streamline regularization will be presented as a way of enforcing a smoothness
constraint on the fibers. Section 5 contains the linearization and discretization
of the PDEs involved in data denoising and fiber regularization. The numerical
methods used to solve linear equations will be described. In Section 6 we will
present visualization techniques for DT 'IIII data. These include rendered el
lipsoids, streamtubes, LIC and particles. In Section 7 experimental results will
be presented. The data denoising results for rat spinal cord and brain data
sets will also be presented. Fiber tracking results using several visualization
techniques will be presented. Conclusions and possible directions for future
work will be discussed in Section 8.
2 Literature Review
Fundamental advances in understanding living biological systems require de
tailed knowledge of structural and functional organization. This is particu
larly important in the nervous system where anatomical connections deter
mine the information pathways and how this information is processed. Our
current understanding of the nervous system is incomplete because of a lack
of fundamental structural information [21] necessary to understand function.
In addition, understanding fundamental structural relationships is essential to
the development and application of therapies to treat pathological conditions
(e.g. disease or injury).
However, most imaging methods give only an anatomically isolated represen
tation of living tissue because the images do not contain connectivity informa
tion. Such information would allow the identification and correlation of system
elements responding during function. For example in brain trauma, the rela
tionship between anatomy and behavior will only become apparent when we
are able to discriminate the afferent nerve fiber pathways transmitting the
sensation from a stimulus to the brain or the efferent pathways transmitting
impulses from the brain region controlling behavior.
Research during the last few decades has shown that the central nervous sys
tem is able to adapt to challenges and recover some function. However, the
structural basis for this adaptive ability is not well understood. For the en
tire central nervous system, understanding and treating evolving pathology,
such as brain trauma, depends on a detailed understanding of the anatomical
connectivity changes and how they relate to function.
Recently :'.l measurements have been developed to measure the tensor of
diffusion. This provides a complete characterization of the restricted motion
of water through the tissue that can be used to infer tissue structure and
hence the fiber tracts. In a series of papers, Basser and colleagues [26,36]
have discussed in detail general methods of acquiring and processing the com
plete apparentdiffusiontensor of .Mli measured translational selfdiffusion.
They showed that directly measured diffusion tensors could be recast in a
rotationally invariant form and reduced to parametric images that represent
the average rate of diffusion (tensor trace), diffusion anisotropy (relationship
of eigenvalues), and how the diffusion ellipsoid (eigenvalues and eigenvectors)
can be related to the laboratory reference frame. The parametric images [5,35]
of volume ratio, fractional anisotropy, and lattice anisotropy index represent
scalar measures of diffusion that are independent of the lab reference frame and
subject orientation. Therefore, these measures can be used to characterize the
tissue pathology, e.g., ischemia, independent of the specific frame of reference
used to acquire the images. The development of diffusion tensor acquisition,
processing, and analysis methods provides the framework for creating fiber
tract maps based on this complete diffusion tensor analysis [19,24,29,30]. This
has been used to produce fiber tract maps in rat brains [30,47] and to map
fiber tracts in the human brain [24] then the first steps were taken to relate
this structural connectivity to function [19].
The directional properties of diffusion can be characterized by a diffusion
tensor, a 3 x 3 symmetric matrix of real values. In order to calculate the 6
independent components of the tensor, the subject is imaged in seven different
directions with several magnetic field strengths. The StejskalTanner equation
[23]
S = So Ej bjD (1)
allows the diffusion tensor, D, and the T2 weighted image So to be calcu
lated given the samples S and field gradient strength, b. Previous work has
concentrated on smoothing the field of eigenvectors of D. More recent work
[14] has formulated regularization techniques for the tensor field itself, even
constraining the resulting tensors to be positivedefinite. Westin et al. [46]
reported techniques for processing DT. 11 I data using tensor averaging. Dif
fusion tensor averaging is an interesting concept but does not address the issue
of estimating a smooth tensor from the given noisy vectorvalued data. More
recently, Poupon et al. [38] developed a Bayesian formulation of the fiber tract
mapping problem. Prior to mapping the fibers, they use robust regression to
estimate the diffusion tensor from the raw vector valued image data. No im
age selective smoothing is performed in their work prior to application of the
robust regression for estimating the diffusion tensors. Here we will take the
approach of smoothing the observed vectorvalued image S prior to calculating
D. The decision to smooth in this manner will be justified in Section 3.
In summary, the anisotropy of water translational diffusion can be used to
visualize structure in the brain and provides the basis for a new method of
visualizing nerve fiber tracts. Initial results have been very encouraging and
suggest that this approach to fiber mapping may be applied to a wide range
of studies in living subjects.
Random molecular motion (Brownian motion) can cause transport of matter
within a system. Within a volume of water, molecules freely diffuse in all
directions. The water abundant in biological systems is also subject to such
stochastic motion. The properties of the surrounding tissue can affect the
magnitude of diffusion, and the directional properties as well. Tissue can form
a barrier to diffusion, restricting molecular motion.
Within an oriented structure, such as a bundle of axonal fibers, diffusion can
be highly anisotropic. The white matter of the brain and spinal cord is char
Fig. 1. Diffusion Ellipsoid
acterized by many such bundles.
The directional properties of diffusion can be characterized by a tensor. The
diffusion tensor, D, is a symmetric, positivedefinite 3 x 3 matrix. We will
make use of the eigenvalues and eigenvectors of this tensor, sorting the eigen
values (A, A2, A3) from largest to smallest, and labelling the corresponding
unit eigenvectors (el, e2, e3). The eigenvalues represent the magnitude of diffu
sion in the direction of their corresponding eigenvector. For isotropic diffusion
A, = A2 = A3. A popular representation for describing anisotropic diffusion is
the diffusion ellipsoid. This ellipsoid is the image of the unit sphere under the
transformation defined by the tensor, D. The eigenvectors of D form an or
thogonal basis, representing the orientation of the ellipsoid. The length of each
axis of the ellipsoid is the corresponding eigenvalue. For isotropic diffusion, the
diffusion ellipsoid is a sphere.
2.1 DTMRI Acquisition
By carefully designing gradient pulse sequences, the measured signal from
protons in water molecules undergoing diffusion can be attenuated. The first
gradient pulse induces a known phase shift in proton precession. After some
delay, a second gradient pulse is applied, inducing the opposite phase shift.
Protons which have not moved between the two gradient pulses are returned
to their previous phase. Protons belonging to molecules which have changed
location have some net change in phase, changing their T2 relaxation time.
Conventional '.11 1 images of white matter in the brain suggest a material of
homogeneous composition. The fibrous nature of white matter is apparent in
DT'. 11I however.
2.2 Stejskal Tanner Equation
The intensity of the received signal, S, at each voxel depends on the properties
of the diffusionencoding gradient, and the apparent diffusion tensor, D at that
location. The StejskalTanner equation relates all of these quantities
S E (2)
SO = e' tlbi'jDiJ (2)
So
In (2), So is the intensity with no diffusionencoding gradient present, and b is
a matrix characterizing the gradient pulse sequence. The physics behind (2) is
beyond the scope of this paper. The motivated reader may refer to the work
of Haacke et al. [23] for more details.
This imaging process must be performed with 7 noncoplanar gradient direc
tions in order to fully generate a diffusion tensor image. Multiple samples,
usually 3 or 4, with varying gradient strengths are taken for each gradient
direction.
In So
Dxx
In S 1 bx bl b 2b' 2b 2b' D,
D,, (3)
In S28 1  1 ?" 2b28 2b 2 2b D
Dyz
DX
The overconstrained linear system, (3) is solved for So and the elements of the
symmetric tensor D by a least squares linear regression [?]. Noisy Si values
can cause the estimated D to be erroneous. In the following sections we will
describe possible methods of denoising vectorvalued data.
2.3 PDE Based Image Restoration
Image data smoothing or denoising is a fundamental problem in image pro
cessing. Image denoising (noise removal) is a technique that enhances images
by attempting to reverse the effects of degradations occurring during acquisi
tion or transmission. Image noise makes it difficult to perform other processing
tasks such as edgedetection, segmentation, or in our case, fiber tracking. For
this reason, denoising is the first step in most image analyses.
The goal of image denoising is to recover an unknown, ideal image given an
observed image. In our case, we wish to recover smooth S values from which
the DT image is calculated. The most common degradation source is the
noise from the image acquisition system and is commonly modeled by additive
Gaussian random noise. In the following, we will briefly review representative
schemes, specifically, partial differential equation (PDE) based methods that
lend themselves to fast numerical implementations.
2.3.1 ScalarValued Images
Image denoising can be formulated using variational principles which in turn
require solutions to PDEs. For scalar valued image smoothing using non
linear diffusion filters with scalar diffusivity coefficient, we refer the reader
to the following articles and references therein [1,11,12,15,28,31,32,34,43,45].
Anisotropic diffusion filters that use a tensorial diffusivity parameter were
introduced in Weickert [45]. These filters can be tailored to enhance image
structures (edges, parallel lines, curves etc.) that occur in preferred directions.
The total variation (TV) norm introduced in Rudin et.al.[40] is a popular
norm used for image restoration. In the case of scalar valued functions, min
imizing the TVnorm amounts to minimization of the L1 norm of the image
gradient. Minimizing the TV norm produces very smooth images while per
mitting sharp discontinuities between regions [40]. Since most images consist
of smooth regions separated by discontinuities (edges), this is a useful model
for image denoising. The scalarvalued functions, the TVnorm is defined as
TVf, (j(x)) = VI(x)dx, C (4)
The EulerLagrange equation for the minimization of 4 written in gradient
descent form is
9I(x,t) VI(x, t)
at i I(x, t)
I(x,0) Iox) (5)
An example of image denoising by TV norm minimization is shown in 2. The
restored image on the right is characterized by a high degree of intraregion
smoothing, and preservation of edges.
Fig. 2. Noisy image (left) and restored (right) by TV norm minimization.
2.3.2 Vector Valued Images
The main application of vectorvalued image restoration thus far in litera
ture has been the restoration of color images. The simplest implementation
of vectorvalued smoothing is uncoupled, also known as channelbychannel,
smoothing. By treating each component of the vector field as an independent
scalar field, we can proceed by smoothing channelbychannel using a scalar
image smoothing technique. This can, however, result in a loss of correlation
between channels as edges in each channel may move independently due to
diffusion. To prevent this, there must be coupling between the channels.
Whitaker and Gerig introduced anisotropic vectorvalued diffusion which was
a direct extension of the work of Perona and Malik [34]. Sapiro et al. [41]
introduced a selective smoothing technique where the selection term is not
simply based on the gradient of the vector valued image but based on the eigen
values of the Riemannian metric of the underlying manifold. In 2 dimensions,
the underlying manifold can be thought of as the parametric surface described
by the image I(x). By geometrically smoothing this surface, we also smooth
the vectorvalued image.
More recently, Kimmel et al. [25] presented a very general flow called the
Beltrami flow as a general framework for image smoothing and show that
most flowbased smoothing schemes may be viewed as special cases in their
framework. There are many other PDEbased image smoothing techniques
which we will not cover here, but will refer the reader to Weickert [45] and
Caselles et al.[10].
Blomgren and Chan [8] introduced the TVn,m norm for vector valued images.
TV,,m(I(x)) = E[TVn,l([i)]2 (6)
For m = 1, equation 6 reduces to the scalar TV norm given in equation 4.
The EulerLagrange condition for the minimization of 6, written in gradient
descent form is
9I (x, t) TVV,(J) V i
(X, t) TV,(I) V ( v I(x, 0) = I(x) (7)
Ot TVn,,m(I) VIi
This was shown to be quite effective for color images, preserving edges in
the color space while attenuating noise. However, for much larger dimensional
data sets (m=7) as in the work proposed here, the Color TV method becomes
computationally very intensive and thus may not be the preferred method in
such applications.
3 DTMRI Image Restoration
The degradation associated with the measurement of the S values is modelled
by an additive gaussian process. Let S(X) be the vector valued image that
needs to be smoothed where, X = (x, y, z) and let S(X) be the unknown
smooth approximation of the data to be estimated. We have S(X) = S(X)+q.
Although we can consider the noise to be on the components of the tensor,
D, the form of this noise is no longer so simple. In fact, even computing D
from S, using the StejskalTanner relation may be meaningless. Substituting
the noise model into (2).
3 3
n S = in So exp( E bijDi,j + q (8)
i=1 j=
The physical principles governing diffusion do not allow for negative values of
S, as evidenced by the dependence of (8) on the natural logarithms of S and
So.
Lowintensity S measurements may be overwhelmed by noise. However, since q
is a zeromean distribution, these measurements may become negative. Clearly,
we cannot calculate D in this case. Nor should we propagate the noise model
through the SkeskalTanner equation in order to write D as a function of q. We
propose to instead smooth the vector of S measurements before any further
analysis.
Smoothing the raw vector valued image data is posed as a variational principle
involving a first order smoothness constraint on the solution to the smoothing
problem. We propose a weighted TVnorm minimization for smoothing the
vector valued image S.
The variational principle for estimating a smooth S(X) is given by
7 7
minE(S)= [g(A+,A _) VSi + Si S]d (9)
s i=i i=1i
where, Q is the image domain and / is a regularization factor. The first term
here is the regularization constraint on the solution to have a certain degree of
smoothness. The second term in the variational principle makes the solution
faithful to the data to a certain degree. The weighting term in this case g(s) =
1/(1 + s) where s = FA is the fractional anisotropy defined by Basser et al.
[2]. This selection criteria preserves the dominant anisotropic direction while
smoothing the rest of the data. Note that since we are only interested in
the fiber tracts corresponding to the streamlines of the dominant anisotropic
direction, it is apt to choose such a selective term.
Here we have used a different TV norm than the one used by Blomgren and
Chan [8]. The TVn,m norm is an L2 norm of the vector of TVn,1 norms for each
channel. We use the L1 norm instead.
The gradient descent of the above minimization is given by
s div (g(A _)VS Si Si) i 1,...,7
at \ Ivs I )
OS
i anxa+ = 0 and S(x,t = 0) = S(x) (10)
On
The use of a modified TV norm in (9) results in a looser coupling between
channels than the use of the true TVn,m norm would have. This reduces the
numerical complexity of (10) and makes solution for large data sets computa
tionally feasible.
Note that the TVn,m norm appears in the gradient descent solution 7 of the
minimization problem. Our data sets consist of 7 images, corresponding to
gradient directions. Each of these images consists of several samples (usually
3 or 4) corresponding to different gradient strengths. Calculating the TVn,m
norm by numerically integrating over the 3dimensional data set at each step
of an iterative process would have been prohibitively expensive.
4 Neuronal Fiber Tracking
Water in the brain preferentially diffuses along white matter fibers. By track
ing the direction of fastest diffusion, as measured by 'I111I noninvasive fiber
tracking of the brain can be accomplished. Fibers tracks maybe constructed
by repeatedly stepping in the direction of fastest diffusion. The direction along
which the diffusion is dominant corresponds to the direction of eigen vector
corresponding to the largest eigenvalue of the tensor D. In Conturo et al.,
[18], fiber tracks were constructed by following the dominant eigenvector in
0.5 mm steps until a predefined measure of anisotropy fell below some thresh
old. This usually occurred in grey matter. The tensor, D, was calculated at
each step from interpolated DT'.II1l data. This tracking scheme is primarily
based on heuristics and is not grounded in well founded mathematical princi
ples. Mori et al. [30] achieved fiber tracking by using several heuristics. The
tracking algorithm starts from a voxel center and proceeds in the direction
of the major axis of the diffusion ellipsoid. When the edge of the voxel is
reached, the direction is changed to that of the neighboring voxel. Tracking
stops when a measure of adjacent fiber alignment crosses a given threshold.
The scheme however is resolution dependent since the '.11 I data only reflects
average axonal orientation within a voxel. Small fibers adjacent to each other
may not be distinguished. Another problem occurs with branching fibers. This
method will only track one path. In this situation, multiple points within a
bundle may be independently tracked. Coulon et al. [20] determined fibers
tracts after smoothing the eigen values and vectors. Once again, this scheme
is faced with the same problem i.e., the eigen vector and the eigen values are
computed from a noisy tensor field and hence may not be meaningful at sev
eral locations in the field. Parker et al. [33] also presented a follow up article
on fiber tract mapping wherein, they use the idea of the fast marching method
from Sethian et al. [42] for growing seeds initialized in the smoothed domi
nant eigen vector field. Batchelor et al. [7] reported an interesting fiber tract
mapping scheme wherein, they produce a map indicating the probability of a
fiber passing through each location in the field. However, they do not address
the issue of computing the diffusion tensor from a noisy set of vector valued
images.
In the context of the fiber tract mapping, two subtasks are involved namely,
(1) estimating a denoised diffusion tensor from noisy vectorvalued image mea
surements and (2) estimating the streamlines of the dominant eigen vectors
of the diffusion tensor. The denoising involves a weighted TVnorm minimiza
tion of the vector valued data S followed by a linear least squares estimation
of D from the smoothed S and then an estimation of the stream lines of the
dominant eigen vectors of D.
Diffusion at each point can be characterized by a diffusion ellipsoid. The el
lipsoid axes are the eigenvectors of the diffusion tensor. The radii are the
corresponding eigenvalues. The shape of the ellipsoid reflects the isotropy of
diffusion. Nearly spherical diffusion ellipsoid represent areas of free water,
where diffusion is unimpeded. Areas of whitematter fiber bundles have elon
gated ellipsoids, since water diffusion is restricted in directions perpendicular
to fiber direction. The phenomenon was quantified by Basser and Pierpaoli [2]
using a scalar measure called the fractional anisotropy, given by
FA 3 (A )2 + (A2 A)2 + (A3 )2 (11)
FA = ++ (11)
2 Al +A +A
This quantification of anisotropy is used as a stopping criterion for streamline
integration, since we are interested in white matter fibers.
Given the dominant eigen vector field of the diffusion tensor in 3D, tracking
the fibers (space curves) tangential to this vector field is equivalent to finding
the stream lines/integral curves of this vector field. Finding integral curves
of vector fields is a well researched problem in the field of Fluid '.1., liiil
[17]. We numerically integrate the given vector field using fourth order Runge
Kutta integration [39].
5 Numerical Methods
In this section we present the details of the numerical solution method em
ployed to solve equation (10). The nonlinear PDE for denoising the raw data
is
VS
div(g ) (S S) = 0 (12)
VS
where g : R3 + IR is the selective smoothing and p is the constant regulariza
tion parameter.
Equation 12 is nonlinear due to the presence of Vs in the denominator of
the first term. We linearize 12 by using the method of l;i,,.I diffusivity"
presented by Chan and '..ilIt [13]. By considering Vs to be a constant for
each iteration, and using the value from the previous iteration we can instead
1
V (Vg Vst + gAst+l) + /(st+l  s) = 0
Vsi
(13)
Here the superscript denotes iteration number. FiNI, rewrite 13 with all of
the unknown st+1 terms on the righthand side
Ast+ + P
V st
9
/i Vstlso + Vg Vst
9
(14)
To write 14 as a linear system discretize the laplacian and gradient terms.
Using central differences for the laplacian we have
Ast+1 =t+l ,1 + l +, t+
X al z xyl z xyz1
6 + St+l + t+l + St+l (15)
Define the standard forward, backward and central differences to be
5XY1,z SXNyz
5x~yz1 sx~yz
5J~yz J 1,ly~z
8 J'YZ SJ;,Y 1,Z
1,
2 (sx+i,y,z
DYS 1.
2
1,
Dzs = 2(sxyz+1
2_ l
We can rewrite 14 in discrete form using the definitions in 16
Sdl,y,z Sd,yl,z
Sx,,z1 + (6 +
SV(Dxst)2 + (DySt)2
9
+ (D ..x)2
8x+l,y,z Sx y+1,z Sx,y,+l
1(ls (Dxst)2 + (Dyst)2 + (D .)2 + DxgDxst + DngDyst + DgD())
9
solve
(16)
sxlyz)
sx_1,,yz)
This results in a bandeddiagonal linear system with 7 nonzero coefficients per
row.
6+ 1 1 1 s+1 fot
1 6+ i1viS1 1 1 1 St+ ft
91 = (18)
0 1 6+ Mvs3 1 1 1
93
S. ". t+t
S_ n3 Af3.
where the righthand side of 17 has been replaced with ft. The matrix in
equation 18 is symmetric and diagonally dominant. We have successfully used
conjugate gradient to solve this system.
The solution of 18 represents one fixedpoint iteration. This iteration is con
tinued until st st+l < c, where c is a small constant, typically of value
1 x 103.
6 Neu ronal Fiber Visualization
A very simple visualization strategy is to simply render the diffusion ellipsoid
at a subset of data points. Since a 3D field of ellipsoids would occlude each
other, this visualization is usually done for 2D slices of data. Additionally,
only ellipsoids on a sparse grid can be rendered in order for each ellipsoid to
be discerned. This type of visualization can easily become visually cluttered
and convey so little information as to be useless. Laidlaw[26] however, suc
cessfully incorporated ellipsoids in a layered visualization approach. In this
paper we present 3 possible methods of visualizing the fibers. In descending
order of computational complexity, they are (1) line integral convolution, (2)
streamtubes, and (3) particles.
Streamtubes are a threedimensional alternative to streamlines. In fact, the
streamtube is computed by using a streamline as the centerline of the tube. We
can use the streamtube diameter to encode some additional information about
the tensor field being visualized, such as the FA value. Previously, Laidlaw [27]
has applied the streamtube visualization approach to DT. 11I
6.1 Line Integral Convolution
It is also possible to visualize the 3D vector field corresponding to the dominant
eigenvalues of the diffusion tensor using other visualization methods such as
the line integral convolution technique introduced by Cabral et al. [9] a
concept explored in this work as well. The advantage of this visualization
technique is that it is well suited for visualizing high density vector fields
and does not depend on the resolution of the vector field moreover, it also
has the advantage of being able to deal with branching structures that cause
singularities in the vector field. Since the fiber direction is parallel to the
Fig. 3. LIC visualization of synthetic field.
dominant eigenvector of the diffusion tensor, we can calculate fiber paths as
integral curves of the dominant eigenvector field. The stopping criterion is
based on FA value. When FA falls below 0.17 we consider the diffusion to be
nearly isotropic and stop tracking the fiber at this point.
Once the diffusion tensor has been robustly estimated, the principal diffu
sion direction can be calculated by finding the eigen vector corresponding to
the dominant eigen value of this tensor. The fiber tracts may be mapped by
visualizing the streamlines through the field of eigen vectors.
LIC is a texturebased vector field visualization method. The technique gen
erates intensity values by convolving a noise texture with a curvilinear kernel
aligned with the streamline through each pixel, such as by
so+L
I(xo) = T(a(s))k(so s)ds (19)
soL
where I(xo) is the intensity of the LIC texture at pixel x0, k is a filter kernel
of width 2L, T is the input noise texture, and a is the streamline through
point x0. The streamline, a can be found by numerical integration, given the
discrete field of eigen vectors.
The result is a texture with highly correlated values between nearby pixels on
the same streamline, and contrasting values for pixels not sharing a streamline.
In our case, an FA value below a certain threshold can be a stopping criterion
for the integration since the diffusion field ceases to have a principal direction
for low FA values. Stalling and Hege [44] achieve significant computational
savings by leveraging the correlation between adjacent points on the same
streamline. For a constant valued kernel, k, the intensity value at I(a(s + ds))
can be quickly estimated by I(a(s))+e, where e is a small error term which can
be quickly computed. Previously, Cliiiiig et al. [16] have used LIC to visualize
fibers from diffusion tensor images of the myocardium.
We adapt this technique to tensor field visualization by incorporating the FA
value at each field location in to the LIC texture. By modulating the image
intensity with an increasing function of FA, we highlight the areas of white
matter, and are able to resolve where the fibers eventually track into grey
matter.
6.2 Particles
The LIC and streamtube techniques presented in the previous sections are
computationally expensive methods. All of the streamlines need to be com
pletely traced before an image can be displayed. For interactive display of
fibers, we developed a novel variant of a particle based visualization technique.
Rather than simulating the actual diffusion process, the particles simply ad
vect through a velocity field described by the dominant eigenvector of the
diffusion tensor at each point. By seeding a few streamlines within a region
of interest, and performing a single step of numerical integration at a time,
interactive framerates can be achieved.
Each particle is implemented as a small textured quadrilateral which is always
oriented to face the viewer. We vary the size and color of this quad as a means
of visualizing the FA value at the particle's location in the tensor field.
Unlike LIC and streamtube tracing, the particle visualization requires no pre
processing time. The other techniques require completely integrating many,
perhaps thousands, of streamlines prior to visualization. On the contrary,
particlebased techniques allow immediate visualization.
7 Experimental Results
Our data were acquired at the McKnight Brain Institute on a 17.6 Tesla (750
MHz) Bruker Avance Imaging spectrometer system with a diffusion weighted
spin echo pulse sequence. Imaging parameters were : effective TR = 3000
ms, TE = 27.7 ms, A = 17.8 ms, 6 = 2.4 ms, ,, ,, = 2250 s/mm2, blow =
750s/mm2, Diffusion gradient strengths = 335, 475, . , mT/m, 6 averages
were performed. The image field of view was 15 x 15 x 19.5 mm3, acquisition
matrix was 128 x 128 x 78. Total time of acquisition was approximately 14
hours.
7.1 Data Denoising
In all of the experiments, we first smooth the seven 3D directional images
using the novel selective smoothing technique outlined in section 3. Following
this, the diffusion tensor is estimated from the smoothed data using a standard
least squares technique. The results of FA calculation from the smoothed data,
and from raw data are presented in Figures 4 and 5
Fig. 4. FA image of coronal slice of raw rat brain data.
Fig. 5. FA results for smoothed data.
As expected from a TVnorm minimization technique, the smoothed images
show a reduction in noise level while retaining edges.
7.2 Fiber Tracking
Figure 6 depicts the computed fiber tracts for the reconstructed rat brain
data. The intensity of the LIC texture has be modulated with the FA image
to emphasize the most anisotropic region of each image. A constantvalued
kernel 188 pixels wide was used to calculate the LIC value at each pixel in the
images.
Fig. 6. LIC fiber tracts in coronal (1) and axial (r) slices of raw rat brain data.
The LIC images show the expected fiber paths in the corpus callosum. Addi
tionally, the LIC image computed from the smoothed data shows fiber paths
computed in low FA regions of the brain.
The streamtube maps 8, 9 were generated by starting streamline integration
at each point of a sparse grid within the data satisfying FA > 0.3. Tracking
stopped at FA < 0.17, the presumed criteria for differentiating white matter
from grey matter. The streamtube radius is a function of FA.
The streamtube images of the corpus callosum agree with atlas images of
corpus callosum dissections when inspected visually.
Figure 10 shows a fiber tracing sequence obtained by seeding fiber in the
corpus callosum region of the rat brain data. The brighter particles signify
high FA value. Dimmer particles can be seen tracing into grey matter. The
Fig. 7. LIC fiber tracts in coronal (1) and axial (r) slices of smoothed rat brain data.
Fig. 8. Axial view of streamtubes in corpus callosum.
particlebased visualization shows the expected fiber paths through the cor
pus callosum, and required the least preprocessing of the other visualization
techniques. This visualization technique allows for interactive selection of the
Fig. 9. Detail of corpus callosum.
Fig. 10. Sequence from Fiber visualization (2 seconds between images).
region of interest for fiber tracking as well as selection of the point of view these
fibers are visualized from. This degree of realtime interactivity is impossible
with the other visualization techniques.
8 Conclusions
In this paper, we presented a new weighted TVnorm minimization formu
lation for smoothing vectorvalued data specifically tuned to computation of
smooth diffusion tensor MR images. The smoothed vector valued data was
then used to compute a diffusion tensor image using standard least squares
technique. Fiber tracts were estimated using the dominant eigen vector field
obtained from the diffusion tensor image. Finally, results of fiber tract map
ping of a rat spinal cord, and a rat brain were depicted using LIC, streamtube
and particles. The fiber tracts are quite accurate when inspected visually, and
correspond well with known anatomical structures, such as the corpus callo
sum.
However, quantitative validation of the computed fiber tracts is essential.
This can be achieved by registration of our threedimensional fibers with two
dimensional fluoroscopy images. This will be the focus of our future efforts.
References
[1] L. Alvarez, P. L. Lions, and J. M. Morel, "Image selective smoothing and edge
detection by nonlinear diffusion. ii," SIAM J. Numer. Anal., vol. 29, no. 3, pp.
845866, June 1992.
[2] P. J. Basser and C. Pierpaoli "Microstructural and physiological features of
tissue elucidated by quantitativediffusiontensor MRI," J. .1n ,,i Reson. B vol.
110, 209219, 1996.
[3] P. J. Basser, J. Mattiello and D. LeBihan M diffusion tensor spectroscopy
and imaging" B,,',li,, J. vol. 66, 259267, 1994a.
[4] P. J. Basser, J. Mattiello and D. LeBihan "Estimate of the effective selfdiffusion
tensor from the NMR spin echo," J. .1?,t.i, Reson. B vol. 103, 247254, 1994b.
[5] P. J. Basser "Inferring microstructural features and the physiological state of
tissues from diffusion weighted images" NMR in Biomed. vol. 8, 333344, 1995.
[6] P. J. Basser, S. Pajevic, C. Pierpaoli, J. Duda and A. Aldroubi, "In vivo fiber
tractography using DTMRI data," .1ni,,t. t,, Resonance in Medicine, vol. 44,
2000, pp. 625632.
[7] P. G. Batchelor, D. L. G. Hill, F. Calamante and D. Atkinson, "Study
of connectivity in the brain using the full diffusion tensor from MRI," in
Proceedings of the 17th Intl. Conf. on Information Processing in Medical
Imaging, 2001, Davis, CA, pp. 121133.
[8] P. Blomgren and T. F. Chan,"Color TV: total variation methods for restoration
of vectorvalued images," IEEE I,,,.u,'i i,,., on Image Processing, vol. 7, no. 3,
pp. 304309, March 1998.
[9] B. Cabral, "Imaging vector fields using line integral convolution," Computer
Graphics Proceedings, pp. 263270, 1993.
[10] V. Caselles, J. M. Morel, G. Sapiro and A. Tannenbaum,IEEE 1 ',,. on
Image Processing, special issue on PDEs and geometrydriven diffusion in image
processing and analysis, vol. 7, No. 3, 1998.
[11] F. Catte, P. L. Lions, J. M. Morel, and T. Coll, "Image selective smoothing and
edge detection by nonlinear diffusion," SIAM Journal on Numerical I.1. ,,t,.
vol. 29, pp. 182193, 1992.
[12] P. Charbonnier, L. BlancFeraud, G. Aubert, and M. Barlaud, "Two
deterministic halfquadratic regularization algorithms for computed iii;,ij..
in in Proc. of the IEEE Intl. Conf. on Image Processing (ICIP), 1994, vol. 2,
pp. 168172, IEEE Computer Society Press.
[13] T. Chan and P. Mulet, "On the convergence of the lagged diffusivity fixed
point method in total variation image restoration," SIAM Journal on Numerical
Analysis, vol. 36, no. 2, pp. 354367, 1999.
[14] C. Chefd'hotel, D. Tschumperl, R. Deriche, O. Faugeras, "Constrained flows
of matrixvalued functions : Application to diffusion tensor regularization,"
in Proc. of the European Conference on Computer Vision 2002 (ECCV'02),
Copenhagen, Denmark, May 2002, pp. 251265.
[15] Y. Chen, B. C. Vemuri and L. Wang, "Image denoising and segmentation via
nonlinear diffusion," Computers and Mathematics with Applications, vol. 39,
No. 5/6, pp. 131149, 2000.
[16] P.J. Chiang, B. Davis and E. Hsu, "Lineintegral convolution reconstruction
of tissue fiber architecture obtained by MR diffusion tensor imaging," Bn [..q
Annual .If,, .,.,t. 2000.
[17] A. Chorin, Computational Fluid Mechanics, Selected papers, Academic Press,
NY, 1989.
[18] T. E. Conturo, R. C. McKinstry, E. Akbudak, and B. H. Robinson "Encoding
of anisotropic diffusion with tetrahedral gradients: A general mathematical
diffusion formalism and experimental results" .1 f,,t Reson. Med. vol. 35, 399
412, 1996
[19] T. E. Conturo, N. F. Lori, T. S. Cull, E. Akbudak, A. Z. Snyder, J. S.
ShJinj .Iy R. C. McKinstry, H. Burton and M. E. Raichle "Tracking neuronal
fiber pathways in the living human ',r;IiJ" Proc. Natl. Acad. Sci. USA vol. 96,
1042210427, 1999.
[20] 0. Coulon, D. C. Alexander and S. R. Arridge, "A regularization scheme for
diffusion tensor magnetic resonance images," in Proceedings of the 17th Intl.
Conf. on Information Processing in Medical Imaging, Davis, CA, pp. 92105,
2001.
[21] F. Crick and E. Jones "Backwardness of human neuroanatomy," Nature vol.
361, 109110, 1993.
[22] L.C.Evans, Partial Differential Equations, American Mathematical Society,
Providence RI, 1998.
[23] E. Haacke, R. Brown, M. Thompson, R. Venkatesan. .,I1,II. I, Resonance
ImagingPhysical Principals and Sequence Design, John Wiley and Sons, NY,
1999.
[24] D. K. Jones, A. Simmons, S. C. R. Williams and M. A. Horsfield, "\ N,.ji_ ov;iive
assessment of axonal fiber connectivity in the human brain via diffusion tensor
MRI" .v1,, Reson. Med., vol. 42, 3741, 1999.
[25] R. Kimmel, N. Sochen, and R. Malladi, "Images as embedding maps and
minimal surfaces:movies, color and volumetric medical images," in Proc. of the
IEEE Conf. on Computer Vision and Pattern Recognition, pp. 350355, June
1997.
[26] D. H. Laidlaw, E. T. Ahrens, D. Kremers, M. J. Avalos, `%i11,.1. spinal cord
diffusion tensor visualization using concepts from painting," Siggraph Technical
Slide Set, 1998.
[27] D. H. Laidlaw, S. Zhang, C. Curry, D. Morris, "Streamtubes and streamsurfaces
for visualizing diffusion tensor MRI volume images," in Proc. of the IEEE
Visualization, October 2000.
[28] R. Malladi and J. A. Sethian, "A unified approach to noise removal, image
enhancement and shape recovery," IEEE l1,~. on Image Processing, vol. 5,
no. 11, pp. 15541568, 1996.
[29] N. Makris, A. J. Worth, A. G. Sorensen, G. M. Papadimitriou, O. Wu, T. G.
Reese, V. J. Wedeen, T. L. Davis, J. W. Stages, V. S. Caviness, E. Kaplan, B.
R. Rosen, D. N. Pandya, and D. N. Kennedy .lph.i jct 'ry of in vivo human
white matter association pathways with diffusionweighted magnetic resonance
imaging," Ann. Neurol., vol. 42, pp. 951962, 1999.
[30] S. Mori, B. J. Crain, V. P. Chacko and P. C. M. van Zijl, I In cdimensional
tracking of axonal projections in the brain by magnetic resonance imaging"
Ann. Neurol., vol. 45, pp. 265269, 1999.
[31] M. Nitzberg and T. Shiota, Ni.lii.;ir image filtering with edge and
corner enhancement," IEEE Transactions on Pattern _i.,'iJi,, and Machine
[I. t, 1 ,,., V,,. vol. 14, no. 8, pp. 826832, 1992.
[32] P. J. Olver, G. Sapiro, and A. Tannenbaum, "Invariant geometric evolutions of
surfaces and volumetric smoothing," SIAM J. Appl. Math., vol. 57, pp. 176194,
1997.
[33] G. J. M. Parker, C. A. M. WheelerKinghtshott and G. J. Baker, "Distributed
anatomical brain connectivity derived from diffusion tensor imaging," in
Proceedings of the 17th Intl. Conf. on Information Processing in Medical
Imaging, Davis, CA, pp. 106120, 2001.
[34] P. Perona and J. Malik, "Scalespace and edge detection using anisotropic
diffusion," IEEE 1,'.' PAMI, vol. 12, no. 7, pp. 629639, 1990.
[35] C. Pierpaoli, and P. J. Basser "Toward a quantitative assessment of diffusion
anisotropy," .Vf,,' Reson. Med. vol. 36, 893906, 1996.
[36] C. Pierpaoli, P. Jezzard, P. J. Basser, A. Barnett and G. Di Chiro "Diffusion
tensor MR imaging of the human brain," Radiology vol. 201, pp. 63'7"; is. 1996.
[37] L. Lapidus and G. F. Pinder, Numerical solution of partial differential equations
in science and engineering, John Wiley and Sons, NY, 1982.
[38] C. Poupon, C. A. Clark, V. Frouin, J. Regis, I. Bloch, D. LeBihan, and J.
Mangin, "Regularization of diffusionbased direction maps for the tracking of
brain white matter fascicles," Neurolmage, vol. 12, 184195, 2000.
[39] W. H. Press, B. P. Flannery, S. A. Teukolsky and W. T. Vetterling, Numerical
recipes in C: The art of scientific computing. Cambridge University Press,
Cambridge, England, second edition, 1992.
[40] L. Rudin, S. Osher, E. Fatemi, N, .!liii;ir total variation based noise reduction
algorithms," Physica D, vol. 60, pp. 2','2(l., 1992.
[41] G. Sapiro and D. L. Ringach, "Anisotropic diffusion of multivalued images with
applications to color filtering," IEEE I1,,. on Image Processing, vol. 5, pp.
15821586, 1996.
[42] J. Sethian, Levelset methods, Cambridge University Press, Cambridge, 1998.
[43] J. Slj;ili, "A common framework for curve evolution, segmentation and
anisotropic diffusion," in IEEE Conf. on Computer Vision and Pattern
Recognition, pp. 136142, 1996.
[44] D. Stalling and H. Hege, "Fast and independent line integral convolution," in
ACM SIGGRAPH, pp. 249256, 1995.
[45] J. Weickert, "A review of nonlinear diffusion filtering,," in Scalespace theory
in computer vision,, (Eds.) B. ter Haar Romney, L.Florack, J. Koenderink, and
M. Viergever, Eds., vol. 1252, of Lecture Notes in Computer Science,, pp. 328,
SpringerVerlag, 1997.
[46] C. F. Westin, S. E. Maier, B. Khidir, P. Everett, F. A. Jolesz and R. Kikinis,
"Image processing for diffusion tensor magnetic resonance imaging," in Proc.
of the Second Intl. Conf. on Medical Image Computing and Computer Assisted
Interventions( IICCAI), pp. .1 111 12, 1999.
[47] R. Xue, P. C. M. van Zijl, B. J. Crain, H. Solaiyappan, and S. Mori "In vivo
threedimensional reconstruction of rat brain axonal projections by diffusion
tensor imaging," ift,,,. Reson. Med. vol. 42, pp. 11231127, 1999.
