Fiber Tract Mapping from Diffusion Tensor MRI *
B. C. Vemuri'
Y Chen2
1Dept. of CISE
M. Rao2
2Dept. of Mathematics
University of Florida
Gainesville, Fl. 32611
UFCISE TR01004
T. McGraw', Z. Wang'
3Dept. of Biochemistry
Abstract
To understand evolving 1,ari. '1i. in the central nervous
system (CNS) and develop cet rti've treatments, it is essen
tial to correlate the nerve fiber connectivity with the vi
sualization of function. Such information is fundamental
in CNS processes since anatomical connections determine
where information is passed and processed. Dirtu.\ion ten
sor imaging (DTI) can provide the fundamental informa
tion required for viewing structural connectivity. However,
robust and accurate acquisition and processing (l1i .0. di1 ,
are needed to accurately map the nerve connectivity. In
this paper we present a novel, ul ... aihl for automatic
fiber tract mapping in the CNS specifically, the spinal cord.
The automatic fiber tract mapping problem will be solved
in two phases, namely a data ',' .. 'rii,, phase and a fiber
tract mapping phase. In the former ,a.'. 'hri,,g; is achieved
via a new weighted TVnorm minimization which strives to
smooth while retaining all relevant detail. Existence and
uniqueness results for this minimization are presented in
brief For the fiber tract mapping, a smooth 3D vector field
.,1. i ... 1ri the dominant anisotropic direction at each spa
tial location is computed from the smoothed data. Fiber
tracts are then determined as the smooth integral curves of
this vector field in a variational framework. To facilitate vi
sualization of the computedfiber tracts, we overlay the 3D
fibers on a volume rendering of the DTMR scan. Exam
ples are presented for DTMR data sets from a normal and
injured rat spinal cords respectively.
1 Introduction
Fundamental advances in understanding living biological
systems require detailed knowledge of structural and func
tional organization. This is particularly important in the
nervous system where anatomical connections determine
*This research was in part funded by the NSF grant IIS9811042 and
NIH RO1RR13197
the information pathways and how this information is pro
cessed. Our current understanding of the nervous system
is incomplete because of a lack of fundamental structural
information [13] necessary to understand function. For the
entire central nervous system, understanding and treating
evolving pathology, such as spinal cord injury, depends
on a detailed understanding of the anatomical connectiv
ity changes and how they relate to function. For example,
an evolving spinal cord lesion undergoes an initial response
to an insult which is followed by a sequence of secondary
events leading to further tissue degradation [37]. A method
for defining the structurefunction relationship is needed
that can be used in the whole living organism that facilitates
the study of such an evolving dynamics process.
Recently MR imaging has been used to study the struc
tural connectivity within whole living organisms. The MR
measurement of water translational selfdiffusion provides
a method that can be used to study structural connectivity
with ubiquitous indigenous material, water. In highly orga
nized nervous tissue, like white matter, diffusion anisotropy
can be used to visualize fiber tracts. Douek, et al. [14] have
used diffusion measurements along three orthogonal axes to
estimate diffusion anisotropy in human brain white matter.
From this data, they produced a color map of the fiber ori
entations reflective of white matter organization. Recently
MR measurements have been developed to measure the ten
sor of diffusion. This provides a complete characteriza
tion of the restricted motion of water through the tissue that
can be used to infer tissue structure and hence fiber tracts.
The development of diffusion tensor acquisition, process
ing, and analysis methods provides the framework for creat
ing fiber tract maps based on this complete diffusion tensor
analysis [12, 17, 19, 21].
For automated fiber tract mapping, prior to estimating
the diffusion tensor, the raw data must be smoothed while
preserving relevant detail. The raw data in this context con
sists of seven directional images acquired for varying mag
netic field strengths. Note that atleast seven values at each
3D grid point in the data domain are required to estimate
T. Mareci3
the six unknowns in the symmetric 2tensor and one scale
parameter. The data smoothing or denoising can be formu
lated using variational principles which in turn require solu
tions to PDEs. Recently, there has been a flurry of activity
on the PDEbased smoothing schemes. In [25], Perona and
Malik developed an anisotropic di/ttu.\in scheme for im
age smoothing. The basic idea of this nonlinear smoothing
scheme was to smooth the image while preserving the edges
in it. This was done by using the following equation Is =
div(c(VI)VI), where I is the image to be smoothed and
It describes its evolution over time, and c(VI) is a decreas
ing function of VI. Catte et al., [5], Nitzberg and Mumford
[22] and Alvarez et al. [1] recognized the illposedness of
the PeronaMalik diffusion and proposed modifications to
overcome the same. Since then, several nonlinear diffusion
methods have been developed and a good account of these
can be found in [5, 1, 8, 23, 31, 33]. All of these are non
linear models and differ in the diffusivity coefficient and/or
the diffusion term. Some of them are also supplemented
with a reactive term. Another popular framework for im
age smoothing is the total variation or TV norm framework
pioneered by Rudin et.al., [28] and further developed by
Chan et.al., [6] and Strong and Chan [32]. The total vari
ation methods yield nonlinear diffusion equations that are
always derived from variational principles using the TV
norm. In [20], Malladi and Sethian propose a unified ap
proach to noise removal and image segmentation using the
concept of minmax curvature flow. Based on the image
data, a min/max switch was designed to select min(n, 0.0)
or max(n, 0.0) so that the curvature based curve evolution
smoothes out small oscillations, but maintains the essen
tial properties of the shape. Results of implementation were
shown on a variety of images yielding quality noise removal
and image segmentation. In [31], Shah developed a com
mon framework for curve evolution, image denoising and
segmentation, and anisotropic diffusion. In this work, a
new segmentation functional was developed which lead to
a coupled system of PDEs, one of them performed nonlin
ear smoothing of the input image and the other smoothed
an "edge strength" function. Shah [31] demonstrated that
all the existing curve evolution and anisotropic diffusion
schemes reported in literature can be viewed as special
cases of his method. In [9], Chen et. al., a nonlinear diffu
sion equation supplemented with reactive terms for achiev
ing edge preserving smoothing was presented. All of the
methods discussed thus far are primarily applicable to the
selective 'i, ,. hri,,i ofscalar valued images.
Smoothing vector valued images has been less popular
than the scalar valued image data sets. In this context,
Whitaker and Gerig introduced anisotropic vectorvalued
diffusion which was a direct extension of the work by Per
ona and Malik [25] to vectorvalued images. The selec
tive term in their work was based on the the gradient of
the vector valued image, which is the Jacobian matrix. In
[30] Sapiro et.al., introduced a selective smoothing tech
nique where the selection term is not simply based on the
gradient of the vector valued image. Instead, he showed
that the stopping term should be a quantity related to the
eigen values of the Riemanian metric tensor computed from
the underlying surface defined by the vector valued image.
They applied their selective smoothing technique to smooth
noisy color images leading to impressive results. A very
general flow called the Beltrami flow as a general frame
work for scalar and vector valued image smoothing was in
troduced in Kimmel et. al., [18] and it was shown that most
flowbased smoothing schemes may be viewed as special
cases in their framework. A generalization of the TV norm
to handle vectorvalued image smoothing was presented in
Blomgren and Chan [3]. They showed that their general
ization was natural and had desirable properties such as the
rotational invariance in the image space etc. Existence and
uniqueness of a solution to their evolution equation is yet to
be explored but is not difficult to establish. There are many
other PDEbased image smoothing techniques that we have
not covered here but will refer the interested reader to a re
cent survey by Weickert [35] and also the special issue of
the IEEE Transactions on Image Processing on PDEbased
image processing [4].
Very briefly, we propose a novel and efficient weighted
total variation (TV) norm based image smoothing scheme
where in the raw image data (one image for each of the 7
directions) S is smoothed using a PDE which is obtained
as a consequence of a weighted TV norm minimization de
fined for vector valued functions. The selective term in
our work, is based on the eigen values of a diffusion ten
sor D that can be computed initially from the raw im
age data using the relationship S = S,,. ifl Eij bijDij,
where, S is the vector of signal/image measurements taken
along seven directions X, Y, Z, XY, YZ, XZ, XYZ, So
is a constant, bij is the magnetic field strength (which is
a constant for a given direction) and Dij are the entries
of the (3, 3) matrix representing the diffusion tensor mea
suring the diffusion of water inside the body being im
aged. The selective term in this case g(s) = 1/(1 + s)
where s = FA is the fractional anisotropy defined as [2]
FA = V+(A ) A) ) where, A1, A3 and
are the largest, smallest and average eigen values of the
diffusion tensor D respectively. This selection criteria pre
serves the dominant anisotropic direction while smoothing
the rest of the data. Another measure that works quite well
is s = (A1 A3)/A3. 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 which correspond to the streamlines of the dominant
anisotropic direction, it is apt to choose such a selective
term as opposed to one that preserves edges in signal in
tensity as was done in [24].
1.1 Finding Stream Lines
Water in the brain preferentially diffuses along white mat
ter fibers. By tracking the direction of fastest diffusion, as
measured by MRI, noninvasive fiber tracking of the brain
can be accomplished. 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
eigen value. In Conturo et. al., [11], fiber tracks were
constructed by following the dominant eigenvector in 0.5
mm steps until a predefined measure of anisotropy fell be
low some threshold. This usually occurred in grey matter.
The tensor, D, was calculated at each step from interpolated
DTMRI data. This tracking scheme is primarily based on
heuristics and is not grounded in well founded mathemati
cal principles. Using methods well grounded in mathemati
cal principles will allow us to better understand/quantify the
strengths and weakness of the method/algorithm leading to
a better overall performance.
In Mori et.al., [21] fiber tracking was achieved using sev
eral 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. One possible measure is the
sum of inner products of nearby data points. This method
was able to reconstruct wellknown pathways through a rat
brain. This method is also based on several data dependent
heuristics for achieving the fiber tract mapping.
Given the dominant eigen vector field of the diffusion
tensor in 3D, tracking the fibers (space curves) along this
dominant eigen vector field is basically equivalent to find
ing the stream lines/integral curves in 3D of this vector field.
Finding integral curves of vector fields is a well researched
problem in the field of Fluid Mechanics [10]. The sim
plest solution would be to numerically integrate the given
vector field using a stable numerical integration scheme
such as a fourth order RungeKutta integrator [27]. How
ever, this may not yield a regularized integral curve. In
this paper, we pose the problem of finding stream lines
of the dominant eigen vector field of the diffusion tensor
in a variational framework incorporating smoothness con
straints which regularize the integral curve. The variational
principle formulation leads to a PDE which can be solved
using efficient numerical techniques. We present the com
puted 3D fiber tracts produced for normal and injured spinal
cords of rats using volume rendering techniques.
2 Image Denoising and Diffusion
Tensor Computation
In this paper, we propose a novel technique for smoothing
vector valued data that will be used in computing the dif
fusion tensor field and mapping out the fiber tracts. The
novelty lies in the formulation that leads to a PDE which is
different from the traditionally used PDEs in literature for
vector valued image selective smoothing. The difference
lies in both the selective term used as well as the fact that
the PDE is derived from a minimization principle which
does not involve arc length minimization as is used tradi
tionally in most selective image smoothing schemes that are
based on minimization principles. Note that there are sev
eral PDEbased schemes in literature that are not based on
minimization principles [31, 34, 29, 20, 36].
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. Let
S(X) be the vector valued image that we want to smooth
where, X = (x, y, z) and let S(X) be the unknown smooth
approximation of the data that we want to estimate. 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
min (X) = gj(AX+, A) IVS(X)I +
s .,i
p/2 5 'dX
i=1
where, Q is the image domain and p is a regularization
factor. The first term here is the regularization constraint
on the solution to have a certain degree of smoothness
and selective smoothing is achieved by the term g(A) =
1/[1 + {(A+ A_)/A+}2], where A is the eigenvalue of the
diffusion tensor computed from the initial data. This func
tion has very small value (approaching zero) as the relative
difference in the largest and smallest eigen values becomes
large stopping the smoothing at such locations and vice
versa. Since, the anisotropy in the image is well captured by
the direction of the dominant eigen value, it is apt for us to
preserve any discontinuities in the anisotropy while smooth
ing the data. Note that it is not the edges (local maxima in
the gradient) in the DT image that we are interested in but
its the anisotropy or the lack thereof that is crucial for the
fiber tract mapping. The aforementioned selective smooth
ing criteria is therefore well justified and is also supported
by the superior quality of the preliminary results (in com
parison to the competing method described in [24]). The
second term in the variational principle makes the solution
faithful to the data to a certain degree. The parameter p con
trols how close the smooth approximation should be to the
given data. The gradient descent of the above minimization
is given by
as., /g(A+, A_)VSj
O =div ft) p(Si S) i = 1,...,7
O anx+ = 0 and S(X,t = 0)= S(X)
(2)
Note that this nonlinear PDE is ditirent from the tradi
tional nonlinear dtirti,.\ ion equations that are found in litera
ture [31, 29, 20] which are curve evolution based schemes.
The main difference is that there is NO V S I multiplica
tive factor in the first term on the right hand side. What dif
ference does this make? Firstly, its asymptotic solution con
verges to the correct minimizer of 1 and the proof of con
vergence does not require the use of the viscosity methods
as in [1]. Moreover, as in the traditional curve/surface evo
lution based nonlinear diffusion equations, if we include the
aforementioned multiplicative factor I S I it is not clear
if the asymptotic solution of the gradient descent equation
2 converges to the true minimizer of 1. The above vari
ational principle is a generalization of the traditional TV
norm for the scalar valued functions but differs in obvious
ways from the generalization presented in Blomgren and
Chan [3]. The main difference being that we use a weighted
TVnorm and our generalization does not have a square root
and the square under the summation as in [3]. The exis
tence and uniqueness questions for such a weighted TV
norm minimization have not been answered in literature to
date. We present a sketch of this proof in the next section.
The above nonlinear PDE can be solved using efficient
and stable numerical schemes. In this paper, we used an im
plicit method namely the CrankNicholson scheme [26]. It
can also be solved using the laggeddiffusivity scheme dis
cussed in [7] but we found the former to be more effective
and stable.
2.1 Estimating the Stream Lines/Integral
Curves
Once the diffusion tensor has been robustly estimated, the
fiber tracts may be mapped by choosing seed points in the
image lattice and using numerical integration techniques to
determine the integral curves of the eigen vector field cor
responding to the dominant eigen values. Several numer
ical integration schemes exist in literature [27]. The most
widely used and stable numerical integration scheme for or
dinary differential equations is the RungeKutta scheme of
order four (RK4) [27]. The solution obtained by directly
using the RK4 may not be at times desirable since there
are no regularization constraints on the resulting integral
curves/fibers which are space curves in this case. In order
to have these space curves not exhibit very sharp twists and
turns (e.g., those with covers), we can formulate the in
tegral curve estimation problem as a variational principle.
Thus, given the eigen vector field v = (vl, v2, v3)T corre
sponding to the dominant eigen values, our formulation of
the variational principle involves minimizing the following
functional:
minE(p) = rmin cic'/) + '(p) v(c(p)) 2dp
(3)
where, c(p) = (x(p), ,/(/'i z(p))T, is the integral curve we
want to estimate and p E [0, 1] is the parameterization of
the curve, Q is the domain over which the curves are to be
determined and v(c(p)) is the vector field v restricted to the
curve c(p). The first term in this functional E(p) is seeking
to minimize the Li norm of the first derivative of the curve
i.e., seeking smooth curves and the second term requires
that the tangent to the smooth curve that we seek be close
to the the given dominant eigen vectors in an L2 sense. The
gradient descent i.e., the EulerLagrange expressed as an
initial boundary value problem, of the variational principle
in equation 3 is given by
ct = c'(p) kn + 3[c"(p) V(c(p))c'(/') +
3V (c(p))(c'(p) v(c(p))
where k is the curvature of the space curve, 3 is a regular
ization parameter and
(Vix Vly Vlz
V = Dv = v2x V 2z
03x V 3z/
VT = the transpose of V. The curvature k is given by
/ c'(p) N
ip \c'Qu)J,,
The above initial boundary value problem can be solved
numerically using a variety of methods. We propose to use
the CrankNicholson implicit method which is a very stable
scheme (see [26]). Each iteration of this numerical iterative
scheme requires the solution a sparse banded (tridiagonal
and positive definite) linear system which can be solved in
O(n) time, where n is the size of the linear system equal
to the number of discrete points on the space curve. As
an initial condition for solving this PDE, we use the results
obtained by simply integrating the given vector field numer
ically using a fourth order RungeKutta method [27]. The
computed integral curves are superimposed on the original
DTI data and visualized using a volume renderer. Results
of this volume visualization are presented in section 4.
3 Existence and Uniqueness Results
In this section, we will briefly outline the approach for es
tablishing the existence and uniqueness of a solution to the
weighted TVnorm minimization equation 1. We will actu
ally present the wellposedness result for a general weighted
TVnorm minimization and our minimization problem (1),
is included in this framework.
Consider the problem
k
min E(u) = min / l{g Vui +
uEBV(n,R"n)nL2 uEBV(n,Rk)nL2 i =1
IUi I,2}
2
(5)
To study the wellposedness of this problem, it is necessary
to introduce the concept of weighted TV norms for func
tions of bounded variation. Recalling definition of bounded
variation (BV) spaces [15, 16].
Definition 1: Let Q C R" be an open set and let f =
(fi,... ,fk) E L(Q,Rk). Define
with ao < g < ai for x E ,, where ao and aI are positive
constants. Define the weighted total variation norm of f(x)
with the weight function g, fo g VfI by
SglV I = glVfi
k
=: sup fi(x)div((x))dx ,
where
) =:{> = 3D( < ,... ,,) e l (,nR")
l0(x))l < g, for all x E a}.
(10)
Theorem 1 Suppose that g satisfies all the assumptions in
Definition 4 and I E L2(, Rk). Then, the minimization
problem (5) has a unique solution u E BV(Q, R") n L2.
Proof Outline 1 The uniqueness of the minimum in (5)fol
lows by the strict convexity of the functional
, +(1/:'l 1 I}.
(11)
jVf
k
Ej^h
k
: Esup
@i= e' l.In
fi(x)div(O(x
where
S=: 0{ = (1, ... d) E Co (, Rn)
1(x) < 1,
for all x E t
Definition 2: A function f E L' (, Rk) is said
have bounded variation in Q, if JI VfI < oo. We def
BV( Rk) as the space of all functions in L1 (Q, Rk) w
bounded variation.
Iff E BV(Q,Rk), Vfi (i = 1,... ,k) then Vfi is
R" valued Radon vector measure. This means
Sfi(x)div(J(x))dx
)) dx (The second term is strictly convex). The existence is proved
) as follows: Let ua, be a minimizing sequence. Then ua,
(6) is bounded in L2 because of the second term in (5). Us
ing the convexity of the functional, we may assume that un
converges in L2 to say u. Then, we use the lower semi
continuity of the norm is BV spaces to conclude that u is a
minimizer
(7) In order to use the gradient descent method for solving
(5), we consider the following evolution problem:
to
n OAu, = div(g(Vu/jV,, ')) (ui Ii), x E Q, t > 0
vith
an
for all E Co(, R").
Definition 3: fm E BV(Q, R) converges weakly in
BV(Q, Rk) to a function f E BV(Q, Rk), if for each i =
1,... ,k
Ui(x,0) = Ii(x), x E O~n
On
(12)
0, x E OQ, t > 0,
(13)
where i = 1,... k, and n is the outward unit norm to 0.
Definition 4: A function u = (U1,... ,U ) E
L2(0, T; BV(Q, Rk)nL2) is called a weak solutionof (12)
(13), if for each i = 1,... ,k Otui E L2(0,T; L2()),
ui = Ii, a.e. in Q, and ui satisfies
lim J Vf mi
) 01Vfi,
./n
for all 0 E C o(Q, R"), where f,m is the ith element of
fm.
Definition 4: Let Q C R" be an open set and f =
(fi,... ,fk) E L k(Q,Rk) and let g be continuous in Q,
\! Ju. rui)+ g vi1
 (ui Ii)(vi ui) f v.,
Q 0 Q
(14)
for a.e. s E [0,T], all v = (v,... ,Vk) E
L2(0, T; BV(Q, Rk) L2). We can now prove the follow
ing theorem.
 f VA(x) ().
Qa
Theorem 2 The problem (12)(13) has an unique solution
u E L2(0, T; BV(Q, Rk)) nL in the sense of(14). More
over, as t + oo u(., t) converges weakly in BV(Q, Rk)nL2
to a function u,, which solves (5).
Proof Outline 2
Otui = divg( Vui)(pv (u j,)
The existence here is well known. We show that the solu
tion Up is indeed in L". We obtain uniform estimates of the
norms to show that the limit as p + 1 of the asymptotic lim
its up + oo exists and is indeed a minimizer i.e., a solution
in BV of(5).
4 Experimental Results
In this section, we present two sets of experiments on, ap
plication of the proposed selective smoothing to the raw im
age data yielding smoothed tensor fields, and the computed
fiber tract maps from these smoothed data for the case of a
normal and an injured rat spinal cord respectively.
In both the experiments, we first smooth the seven 3D di
rectional images using the novel selective smoothing tech
nique outlined in section 2. Following this, the diffusion
tensor is estimated from the smoothed data using a stan
dard least squares technique. The fractional anisotropy, the
color trace of the diffusion tensor (sum of the diagonal terms
of the diffusion tensor in an RGB color space), the domi
nant eigen value as well as the color map of the direction
cosines of the eigen vector corresponding to the dominant
eigen value are computed. The latter color map depicts the
standard axis (X, Y, Z) toward which direction of diffusion
in the data is dominant and the dominant eigen value cor
responds to the magnitude of this dominant diffusion di
rection. Red corresponding to the X component, Green for
the Y component and Blue for the Z component. This color
coding will indicate the standard direction (X, Y or Z) along
which the dominant eigen vector has the strongest compo
nent. Images obtained as a result of these computations
from raw data, smoothed data using a competing smoothing
method outlined in Parker et.al., [24] and smoothed data us
ing our proposed method are depicted for the two data sets.
In addition, estimated fiber tracts from the smoothed data
using the proposed fiber tract/integral curves computation
scheme are also depicted for the data sets. The results of
smoothing for the two examples are shown in Figures 1 and
3 which are organized as follows: first row contains im
ages computed from raw (noisy) data, second row contains
images computed using methods in [24] and third row con
tains computed images using the proposed image smooth
ing technique. Note the superior performance of the pro
posed smoothing scheme in comparison to the method in
Parker et .al., [24].
Figure 2 depicts the computed 3D fiber tracts for a nor
mal rat spine and is organized as follows: On the top left,
computed fiber tracts are shown in green. These fiber tracts
are supposed to run along the length of the spinal cord in the
white matter which is exactly what our computations reveal.
The top right image depicts a volume rendered DTMR data
set with the fiber tracts overlayed in green. The bottom row
shows two different cut away views of the overlayed fiber
tracts on the volume rendered DTMR scan of this normal
rat spinal cord.
Figure 3, depicts the results obtained from an injured rat
spine by application of the proposed smoothing in compar
ison to results obtained from raw data and the competing
method in [24]. Fiber tract mapping is also shown for this
injured spine in figure 4. These results are organized as in
the figures 1 and 2. As evident, the injury has caused a large
cavity down the length of the spine and there are no fibers in
this region. Also evident is the fact that our data smoothing
results are superior to smoothing performed by the Perona
and Malik scheme [25] which was used in [24]. Also, the
visual quality of the fiber tracts is satisfactory.
In the above presented results, what is to be noted is that
we have demonstrated a proof of concept for the proposed
data smoothing and fiber tract mapping algorithms in the
case of the normal and injured rat spinal cords respectively.
The quality of results obtained is reasonably satisfactory for
visual inspection purposes but quantitative validation needs
to be performed and will be the focus of our future efforts.
5 Conclusions
In this paper, we presented a new weighted TVnorm min
imization formulation for smoothing vectorvalued data
specifically tuned to computation of smooth diffusion tensor
MR images. Existence and uniqueness of a solution for the
weighted TVnorm minimization is outlined. The smoothed
vector valued data was then used to compute a diffusion
tensor image using standard least squares technique. Fiber
tracts in 3D were computed as the integral curves of the
dominant eigen vector field obtained from the diffusion ten
sor image. The integral curve computation was formulated
in a variational framework as well and solved using efficient
numerical schemes. Finally, results of fiber tract mapping of
a normal and an injured rat spinal cord were depicted using
standard volume rendering techniques. The computed fiber
tracts are quite accurate when inspected visually. However,
quantitative validation of the computed fiber tracts is essen
tial and will be the focus of our future efforts.
Figure 1: Normal Cord, left to right: FA, Color Trace, di
rection cosines of the dominant eigen vector, and the dom
inant eigen value. First row: results computed from raw
data, Row2: results computed using PeronaMalik diffu
sion, Row3: results from the proposed smoothing.
Figure 2: Normal rat spinal cord: computed fibers over
layed on volume rendered DTMR data and its cut away
views.
Figure 3: Injured Cord, left to right: FA, Color Trace, di
rection cosines of the dominant eigen vector and the dom
inant eigen value. First row: results computed from raw
data, Row2: results computed using PeronaMalik Diffu
sion, Row3: results from the proposed smoothing.
(a) (b)
(c) (d)
Figure 4: Injured rat spinal cord: computed fibers overlayed
on volume rendered DTMR data and its cut away views.
References
[1] L.Alvarez, P. L. Lions, and J. M. Morel, "Image selec
tive smoothing and edge detection by nonlinear dif
fusion. ii," SIAMJ. 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. Magn. Reson.
B 110, 209219 (1996)
[3] P. Blomgren and T F. Chan,"Color TV: Total Varia
tion Methods for Restration of VectorValued Images,"
IEEE Transaction on Image Processing, Vol. 7, no. 3,
pp. 304309, March, 1998.
[4] V Caselles, J. M. Morel, G. Sapiro and A. Tannen
baum,IEEE Trans. on Image Processing, special issue
on PDEs and geometrydriven diffusion in image pro
cessing and analysis, Vol 7, No. 3, 1998.
[5] F.Catte, PL. Lions, J.M. Morel, and T.Coll, "Image
selective smoothing and edge detection by nonlinear
diffusion," SIAM Journal of Numerical Analysis, vol.
29, pp. 182193, 1992.
[6] T. F. Chan, G. Golub, and P. Mulet, "A nonlinear
primaldual method for TVbased image restoration,"
in Proc. 12th Int. Conf Analysis and Optimization of
Systems: Images, Wavelets, and PDE's, Paris, France,
June 2628, 1996, M. Berger et al., Eds., no. 219, pp.
241252.
[7] T. Chan and P. Mulet, "On the Convergence of the
Lagged Diffusivity Fixed Point Method in Total Vari
ation Image Restoration," September 1997, CAM TR
9746.
[8] P.Charbonnier, L.BlancFeraud, G.Aubert, and
M.Barlaud, T\\o deterministic halfquadratic regu
larization algorithms for computed imaging,," in in
Proc. of the IEEE Intl. Conf on Image Processing
(ICIP), 1994, vol. 2, pp. 168172, IEEE Computer
Society Press.
[9] 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,
2000, pp. 131149.
[10] A. Chorin, Computational Fluid Mechanics, Selected
papers, Academic Press, 1989.
[11] T. E. Conturo, R. C. McKinstry, E. Akbudak, and B.
H. Robinson "Encoding of Anisotropic Diffusion with
Tetrahedral Gradients: A General Mathematical Dif
fusion Formalism and Experimental Results" Magn.
Reson. Med. 35, 399412 (1996)
[12] T. E. Conturo, N. F. Lori, T. S. Cull, E. Akbudak, A.
Z. Snyder, J. S. Shimony, R. C. McKinstry, H. Burton
and M. E. Raichle "Tracking neuronal fiber pathways
in the living human brain" Proc. Natl. Acad. Sci. USA
96, 1042210427 (1999)
[13] F. Crick and E. Jones "Backwardness of human neu
roanatomy" Nature 361, 109110 (1993)
[14] P Douek, R. Turner, J. Pekar, N. Patronas and D.
LeBihan "MR color mapping of myelin fiber orienta
tion" J. Comput. Assist. Tomogr. 15, 923929 (1991)
[15] L.C.Evans and R.Gariepy, "Measure theory and fine
properties of functions", Studies in Advanced Mathe
matics, 1992.
[16] E. Giusti, "Minimal surfaces and functions of bounded
variation", Birkhauser, Basel, 1985.
[17] D. K. Jones, A. Simmons, S. C. R. Williams andM. A.
Horsfield, "Noninvasive assessment of axonal fiber
connectivity in the human brain via diffusion tensor
MRI" Magn. Reson. Med., 42, 3741 (1999).
[18] R.Kimmel, N.Sochen, and R.Malladi, "Images as em
bedding maps and minimal surfaces:movies, color and
volumetric medical images,," in Proc. of the IEEE
Conf on Computer Vision and Pattern F.,.... ,,,. ,,i
June 1997, pp. 350355.
[19] N. Makris, A. J. Worth, A. G. Sorensen, G. M. Pa
padimitriou, 0. 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 "Mor
phometry of in vivo human white matter associa
tion pathways with diffusionweighted magnetic res
onance imaging," Ann. Neurol., 42, 951962 (1999).
[20] R.Malladi and J.A. Sethian, "A unified approach to
noise removal, image enhancement and shape recov
ery," IEEE Trans. on Image Processing, vol. 5, no. 11,
pp. 15541568, 1996.
[21] S. Mori, B. J. Crain, V P. Chacko and P C. M. van Zijl
"Threedimensional tracking of axonal projections in
the brain by magnetic resonance imaging" Ann. Neu
rol., 45, 265269 (1999)
[22] M.Nitzberg and T.Shiota, "Nonlinear image filtering
with edge and corer enhancement,," IEEE Transac
tions on Pattern Analysis and Machine Intelligence,
vol. 14, no. 8, pp. 826832, 1992.
[23] P.J. Olver, G.Sapiro, and A.Tannenbaum, "Invari
ant geometric evolutions of surfaces and volumetric
smoothing," SIAM J. Appl. Math., vol. 57, pp. 176
194, 1997.
[24] G.J. M. Parker, J. A. Schnabel, M. R. Symms, D. J.
Werring and G. J. Baker, "Nonlinaer smoothing for
reduction of systematic and random errors in diffu
sion tensor imaging," Magn. Reson. Imag., 11, 702
710, 2000.
[25] P.Perona and J.Malik, "Scalespace and edge detection
using anisotropic diffusion," IEEE Trans. PAMI, vol.
12, no. 7, pp. 629639, 1990.
[26] L. Lapidus and G. F. Pinder, Numerical solution of
partial ditrrrential equations in science and engineer
ing, John Wiley and Sons, 1982.
[27] W.H.Press, B.P.Flannery, S.A.Teukolsky and
W.T.Vetterling, [1992], Numerical Recipes in C:
The Art of Scientific C. ii,,imw' Cambridge Univer
sity Press, Cambridge, England, second edition.
[28] L. I. Rudin, S. Osher, and E. Fatemi, "Nonlinear varia
tion based noise removal algorithms", Physica D, vol.
60, pp. 259268, 1992.
[29] G.Sapiro, "Vectorvalued active contours," in IEEE
Proc. of the CVPR. 1996, pp. 680685, IEEE Com
puter Soc. Press.
[30] G.Sapiro and D.L. Ringach, "Anisotropic diffusion of
multivalued images with applications to color filter
ing," IEEE Trans. on Image Processing, vol. 5, pp.
15821586, 1996.
[31] J.Shah, "A common framework for curve evolution,
segmentation and anisotropic diffusion," in IEEE
Conf on Computer Vision and Pattern F.,...,,,m,. a,
1996.
[32] D. M. Strong and T. F. Chan, "Relation of regulariza
tion parameter and scale in total variation based de
noising," in Proc. IEEE Workshop on Mathematical
Methods in Biomedical Image Analysis, 1996
[33] J.Weickert, "A review of nonlinear diffusion filter
ing,," in Scalespace theory in computer vision,, (Eds.)
B. ter Haar Romney, L.Florack, J. Koenderink, and M.
Viergever, Eds. 1997, vol. 1252, of Lecture Notes in
Computer Science,, pp. 328, SpringerVerlag.
[34] R. Whitaker and G. Gerig, "Vectorvalued diffusions,"
in Geometrydriven Dirtu.\i'on. in Computer Vision, B.
Romney et.al., (Eds.), Kluwer, 1994.
[35] J.Weickert, "A review of nonlinear diffusion filter
ing,," in Scalespace theory in computer vision,, (Eds.)
B. ter Haar Romney, L.Florack, J. Koenderink, and M.
Viergever, Eds. 1997, vol. 1252, of Lecture Notes in
Computer Science,, pp. 328, SpringerVerlag.
[36] A. Yezzi, Jr.,"Modified Curvature Motion for Image
Smoothing and Enhancement," IEEE Transaction on
Image Processing, Vol. 7, no. 3, pp. 345352, March,
1998.
[37] W. Young "Secondary injury mechanisms in acute
spinal cord injury" J. Emerg. Med. 13, 1322 (1993)
