of 9
Level set modeling and segmentation of diffusion
tensor magnetic resonance imaging brain
data
Leonid Zhukov
Ken Museth
David Breen
Alan H. Barr
California Institute of Technology
Department of Computer Science
Mail Code 350-74
Pasadena, California 91125-7400
E-mail:
$
zhukov,kmu,david,barr
%
@gg.caltech.edu
Ross Whitaker
University of Utah
School of Computing
3190 MEB
Salt Lake City, Utah 84112-9205
Abstract.
Segmentation of anatomical regions of the brain is one of
the fundamental problems in medical image analysis. It is tradition-
ally solved by iso-surfacing or through the use of active contours/
deformable models on a gray-scale magnetic resonance imaging
(MRI) data. We develop a technique that uses anisotropic diffusion
properties of brain tissue available from diffusion tensor (DT)-MRI to
segment brain structures. We develop a computational pipeline
starting from raw diffusion tensor data through computation of invari-
ant anisotropy measures to construction of geometric models of the
brain structures. This provides an environment for user-controlled
3-D segmentation of DT-MRI datasets. We use a level set approach
to remove noise from the data and to produce smooth, geometric
models. We apply our technique to DT-MRI data of a human subject
and build models of the isotropic and strongly anisotropic regions of
the brain. Once geometric models have been constructed they can
be combined to study spatial relationships and quantitatively ana-
lyzed to produce the volume and surface area of the segmented
regions.
© 2003 SPIE and IS&T.
[DOI: 10.1117/1.1527628]
1 Introduction
Diffusion tensor magnetic resonance imaging
1–4
~
DT-MRI
!
is a technique used to measure the diffusion properties of
water molecules in tissues. Anisotropic diffusion can be
described by the equation
]
C
]
t
5
π
ï
~
D
π
C
!
,
~
1
!
where
C
is the concentration of water molecules, and
D
is
a diffusion coefficient, which is a symmetric second-order
tensor
D
5
S
D
xx
D
xy
D
xz
D
yx
D
yy
D
yz
D
zx
D
zy
D
zz
D
.
~
2
!
Figure 1 presents a ‘‘slice’’ of the diffusion tensor volume
data of human brain used in our study. Each subimage pre-
sents the scalar values of the associated diffusion tensor
component for one slice of the dataset.
Tissue segmentation and classification based on DT-
MRI offers several advantages over conventional MRI,
since diffusion data contains additional physical informa-
tion about the internal structure of the tissue being scanned.
However, segmentation and visualization using diffusion
data is not entirely straightforward. First, the diffusion ma-
trix itself is not invariant with respect to rotations, and the
elements that form the matrix will be different for different
orientations of the sample or field gradient and therefore
cannot themselves be used for classification purposes.
Moreover, 3-D visualization and segmentation techniques
available today are predominantly designed for scalar and
sometimes vector fields. Thus, there are two fundamental
problems in tensor imaging:
~
1
!
finding an invariant repre-
sentation of a tensor that is independent of a frame of ref-
erence and constructing a mapping from the tensor field to
a scalar or vector field and
~
2
!
visualization and classifica-
tion of tissue using the derived scalar fields.
Paper MIP-13 received May 1, 2001; revised manuscript received Oct. 1, 2001;
accepted for publication Mar. 1, 2002.
1017-9909/2003/$15.00 © 2003 SPIE and IS&T.
Journal of Electronic Imaging 12(1), 125
133 (January 2003).
Journal of Electronic Imaging / January 2003 / Vol. 12(1) / 125
Downloaded From: http://astronomicaltelescopes.spiedigitallibrary.org/ on 09/30/2016 Terms of Use: http://spiedigitallibrary.org/ss/termsofuse.aspx
The traditional approaches to diffusion tensor imaging
involve converting the tensors into an eigenvalue/
eigenvector representation, which is rotationally invariant.
Every tensor can then be interpreted as an ellipsoid with
principal axes oriented along the eigenvectors and radii
equal to the corresponding eigenvalues. This ellipsoid de-
scribes the probabilistic distribution of a water molecule
after a fixed diffusion time.
Using eigenvalues and eigenvectors one can compute
different anisotropy measures
1,5– 8
that map tensor data onto
scalars and can be used for further visualization and seg-
mentation. Although eigenvalue/vector computation of the
3
3
3 matrix is not expensive, it must be repeatedly per-
formed for every voxel in the volume. This calculation eas-
ily becomes a bottleneck for large datasets. For example,
computing eigenvalues and eigenvectors for a 512
3
volume
requires over 20 CPU min on a powerful workstation. An-
other problem associated with eigenvalue computation is
stability—a small amount of noise will change not only the
values but also the ordering of the eigenvalues.
9
Since
many anisotropy measures depend on the ordering of the
eigenvalues, the calculated direction of diffusion and clas-
sification of tissue will be significantly altered by the noise
normally found in diffusion tensor datasets. Thus it is de-
sirable to have an anisotropy measure that is rotationally
invariant, does not require eigenvalue computations, and is
stable with respect to noise. The tensor invariants with
these characteristics were first proposed by Ulug and Zijl.
10
In Sec. 2 of this paper we formulate a new anisotropy mea-
sure for tensor field based on these invariants.
Visualization and model extraction from the invariant
3-D scalar fields is the second issue addressed in this paper.
One of the popular approaches to tensor visualization rep-
resents a tensor field by drawing ellipsoids associated with
the eigenvectors/values.
11
This method was developed for
2-D slices and creates visual cluttering when used in three
dimensions. Other standard computational fluid dynamics
~
CFD
!
visualization techniques such as tensor lines do not
provide meaningful results for the MRI data due to rapidly
changing directions and magnitudes of eigenvector/values
and also amount of noise present in the data. Recently
Kindlmann and Weinstein
12
developed a volume-rendering
approach to tensor field visualization using eigenvalue-
based anisotropy measures to construct transfer function
and color maps, that highlight some brain structures and
diffusion patterns.
In our work, we perform iso-surfacing on the 3-D scalar
fields derived from our tensor invariants to visualize and
segment the data. An advantage of iso-surfacing over other
approaches is that it can provide the shape information
needed for constructing geometric models, and computing
internal volumes and external surface areas of the extracted
regions. A detailed discussion of the modeling method is
presented in Sec. 3. Section 4 presents the results of tensor-
invariant calculations and model segmentation technique
with examples from a DT-MRI scan of a human head. Sec-
tion 5 then describes the quantitative analysis of obtained
geometric models.
Finally, a number of recent publications
13,14
have been
devoted to brain fiber tracking. This is a different and more
complex task than the one addressed in this paper and re-
quires data with a much higher resolution and better SNR
than the data used in our study.
2 Tensor Invariants
Tensor invariants
~
rotational invariants
!
are combinations
of tensor elements that do not change after the rotation of
the tensor’s frame of reference, and thus do not depend on
the orientation of the patient with respect to the scanner
when performing DT imaging. The well-known invariants
are the eigenvalues of diffusion tensor
~
matrix
!
D
, which
are the roots of corresponding characteristic equation
l
3
2
C
1
ï
l
2
1
C
2
ï
l
2
C
3
5
0,
~
3
!
with coefficients
Fig. 1
Slice of a tensor volume where every ‘‘element’’ of the image
matrix corresponds to one component of the tensor
D
.
Fig. 2
Isotropic
C
1
(left) and anisotropic
C
a
(right) tensor invariants
for the tensor slice shown in Fig. 1.
Zhukov et al.
126 / Journal of Electronic Imaging / January 2003 / Vol. 12(1)
Downloaded From: http://astronomicaltelescopes.spiedigitallibrary.org/ on 09/30/2016 Terms of Use: http://spiedigitallibrary.org/ss/termsofuse.aspx
C
1
5
D
xx
1
D
yy
1
D
zz
,
C
2
5
D
xx
D
yy
2
D
xy
D
yx
1
D
xx
D
zz
2
D
xz
D
zx
1
D
yy
D
zz
2
D
yz
D
zy
,
~
4
!
C
3
5
D
xx
~
D
yy
D
zz
2
D
zy
D
yz
!
2
D
xy
~
D
yx
D
zz
2
D
zx
D
yz
!
1
D
xz
~
D
yx
D
zy
2
D
zx
D
yy
!
.
Since the roots of Eq.
~
3
!
are rotational invariants, the co-
efficients
C
1
,
C
2
, and
C
3
are also invariant. In the eigen
frame of reference they can be easily expressed through the
eigenvalues
C
1
5
l
1
1
l
2
1
l
3
,
C
2
5
l
1
l
2
1
l
1
l
3
1
l
2
l
3
,
~
5
!
C
3
5
l
1
l
2
l
3
,
and are proportional to the sum of the radii, surface area
and the volume of the ‘‘diffusion’’ ellipsoid. Then instead
of using (
l
1
,
l
2
,
l
3
) to describe the dataset, we can use
(
C
1
,
C
2
,
C
3
) . Moreover, since
C
i
are the coefficients of
characteristic equation, they are less sensitive to noise, then
roots
l
i
of the same equation.
11
Any combination of the preceding invariants is, in turn,
an invariant. We consider the following dimensionless
combination:
C
1
C
2
/
C
3
. In the eigenvector frame of refer-
ence, it becomes
C
1
C
2
C
3
5
3
1
l
2
1
l
3
l
1
1
l
1
1
l
3
l
2
1
l
1
1
l
2
l
3
,
~
6
!
and we can define a new dimensionless anisotropy measure
C
a
5
1
6
S
C
1
C
2
C
3
2
3
D
.
~
7
!
It is easy to show that for isotropic diffusion, when
l
1
5
l
2
5
l
3
, the coefficient
C
a
5
1. In the anisotropic case,
this measure is identical for both linear, directional diffu-
sion (
l
1
@
l
2
'l
3
) and planar diffusion (
l
1
'l
2
@
l
3
) and
is equal to
C
a
limit
'
1
3
S
1
1
l
1
l
3
1
l
3
l
1
D
.
~
8
!
Thus
C
a
is always
;l
max
/
l
min
and measures the mag-
nitude of the diffusion anisotropy. Note that we use eigen-
value representation here only to analyze the behavior of
the coefficient
C
a
, but we use invariants (
C
1
,
C
2
,
C
3
)to
compute it using Eqs.
~
5
!
and
~
7
!
. Isotropic
C
1
and aniso-
tropic
C
a
tensor invariants maps for the data slice from Fig.
1 is shown in Fig. 2.
3 Geometric Modeling
Two options are usually available for viewing the scalar
volume datasets, direct volume rendering
15,16
and volume
segmentation
17
combined with conventional surface render-
ing. The first option, direct volume rendering, is capable of
supplying only images of the data. While this method may
provide useful views of the data, it is well known that it is
difficult to construct the exact transfer function that high-
lights the desired structures in the volume dataset.
18
Our
approach instead focuses on extracting geometric models of
the structures embedded in the volume datasets. The ex-
tracted models can be used for interactive viewing, but the
segmentation of geometric models from the volume
datasets provides a wealth of additional benefits and possi-
bilities. The models can be used for quantitative analysis of
the segmented structures, for example, the calculation of
surface area and volume; quantities that are important when
studying how these structures change over time. The mod-
els may be used to provide the shape information necessary
for anatomical studies and computational simulation, for
example,
electroencephalogram/magnetoencephalogram
~
EEG/MEG
!
modeling within the brain.
19
Creating separate
geometric models for each structure enables the straightfor-
ward study of the relationship between the structures, even
though they come from different datasets. The models can
also be used within a surgical planning/simulation/VR
environment,
20
providing the shape information needed for
collision detection and force calculations. The geometric
models can even be used for manufacturing real physical
models of the structures.
21
It is clear that there are numer-
ous reasons to develop techniques for extracting geometric
models from diffusion tensor volume datasets.
The most widely used technique for extracting polygo-
nal models from volume datasets is the Marching Cubes
algorithm.
22
This technique creates a polygonal model that
approximates the iso-surface embedded in a scalar volume
dataset for a particular iso-value. The surface represents all
the points within the volume that have the same scalar
value. The polygonal surface is created by examining every
‘‘cube’’ of eight volume grid points and defining a set of
triangles that approximates the piece of the iso-surface
within the space bounded by the eight points. While the
Marching Cubes algorithm is easy to understand and
straightforward to implement, applying it directly to raw
volume data from scanners can produce undesirable results,
as seen in top row images in Figs. 4 and 7. The algorithm
is susceptible to noise and can produce many unwanted
triangles that mask the central structures in the data. To
alleviate this problem, we utilize a deformable model ap-
proach to smooth the data and remove the noise-related
artifacts. Many types of deformable models have been pro-
posed for extracting structures from volumes.
17,23
We uti-
lize level set models as they have been shown to be flexible
and effective for segmentation.
24 –28
Level set methods
produce active deformable surfaces that may be directed
to conform to features in a volume dataset while simulta-
neously applying a smoothing operation based on local
surface
curvature.
28
Most
importantly,
they
easily
change topology during deformation and have no fixed pa-
rameterization, enabling them to represent complex
shapes.
Level set modeling and segmentatio
n...
Journal of Electronic Imaging / January 2003 / Vol. 12(1) / 127
Downloaded From: http://astronomicaltelescopes.spiedigitallibrary.org/ on 09/30/2016 Terms of Use: http://spiedigitallibrary.org/ss/termsofuse.aspx
3.1
Level Set Models
A level set model
29,30
specifies a surface as a level set
~
iso-
surface
!
of a scalar volumetric function,
f
:
U
R
, where
U
,
R
3
is the range of the surface model. Thus, a surface
S
is
S
5
$
s
u
f
~
s
!
5
k
%
,
~
9
!
and
k
is the isovalue. In other words,
S
is the set of points
s
in
R
3
that composes the
k
’ th iso-surface of
f
. The em-
bedding
f
can be specified as a regular sampling on a rec-
tilinear grid. The surfaces may propagate with
~
time-
varying
!
curvature-dependent speeds. Level set methods
provide the mathematical and numerical mechanisms for
computing surface deformations as iso-values of
f
by solv-
ing a partial differential equation on the 3-D grid (
U
) . That
is, the level set formulation provides a set of numerical
methods that describes how to manipulate the gray-scale
values in a volume, so that the iso-surfaces of
f
move in a
prescribed manner
~
see Fig. 3
!
.
There are two different approaches to defining a deform-
able surface from a level set of a volumetric function, as
described in Eq.
~
9
!
. Either one can think of
f
~
s
!
as a static
function and change the iso-value
k
(
t
) or alternatively fix
k
and let the volumetric function dynamically change in time,
i.e.,
f
(
s
,
t
) . Following the second approach, we can math-
ematically express the dynamic model as
f
~
s
,
t
!
5
k
.
~
10
!
To transform this definition into partial differential equation
that can easily be solved by standard numerical techniques,
we differentiate both sides of Eq.
~
10
!
with respect to time
t
, and apply the chain rule:
]
f
~
s
,
t
!
]
t
1
π
f
~
s
,
t
!
ï
d
s
d
t
5
0.
~
11
!
Equation
~
11
!
is sometimes referred to as a ‘‘Hamilton–
Jacobi-type’’ equation and defines an initial value problem
for the time-dependent
f
. Let d
s
/d
t
be the movement of a
point on a surface as it deforms, such that it can be ex-
pressed in terms of the position of
s
P
U
and the geometry
of the surface at that point, which is, in turn, a differential
expression of the implicit function,
f
. This gives a partial
differential equation
~
PDE
!
on
f
:
s
[
s
(
t
):
]
f
]
t
52
π
f
ï
d
s
d
t
[
2
π
f
ï
F
~
s
,D
f
,D
2
f
,...
!
,
~
12
!
where
F
is a user-defined ‘‘speed’’ term which generally
depends on a set of order-
n
derivatives of
f
,
D
n
f
, evalu-
ated at
s
, as well as other functions of
s
. Typically
F
~
x
!
combines a data term with a smoothing term, which pre-
vents the solution from fitting too closely to noise-
corrupted data. There are a variety of surface-motion terms
that can be used in succession or simultaneously in a linear
combination to form
F
~
x
!
. For the work presented in this
paper, we combine a feature attraction term and a smooth-
ing term weighted
28
by a factor
b
,
F
5
F
attr
1
b
F
curv
.
~
13
!
The first term
F
attr
is due to the attraction to the edges in the
volume. It attracts the surface models to certain gray-scale
features in the input data. For instance, the gradient mag-
nitude indicates areas of high contrast in volumes. By fol-
lowing the gradient of such gray-scale features, surface
models are drawn to minimum or maximum values of that
feature. Typically, gray-scale features, such as the gradient
magnitude are computed with a scale operator, e.g., a
derivative-of-Gaussian kernel. If models are properly ini-
tialized, they can move according to the gradient of the
gradient magnitude and settle onto the edges of an object at
a resolution that is finer than the original volume. For this
work we used the attraction force
F
attr
5
π
u
~
π
~
G
*
I
~
x
!!
u
,
~
14
!
where the volume data
I
(
x
) is convolved with a Gaussian
kernel
G
with
s
'
0.5, such that a positive sign moves sur-
faces toward maxima and the negative sign toward minima.
There are a variety of options for the curvature smooth-
ing terms in Eq.
~
13
!
, and the question of efficient, effective
higher order smoothing terms is the subject of on-going
research.
30
For the work presented in this paper the smooth-
ing term uses the mean curvature
K
M
of the level set
S
to
form a vector in the direction of the surface normal
n
:
F
curv
5
K
M
n
5
~
π
ï
n
!
n
5
π
ï
S
π
f
u
π
f
u
D
π
f
u
π
f
u
.
~
15
!
It is weighted by a factor
b
, enabling the user to control the
amount of smoothing, and is tuned for each dataset. The
level set propagation stops when the
F
attr
and
b
F
curv
terms
cancel each other, or when the number of computational
iterations reaches a user-specified value.
Fig. 3
Level set models represent curves and surfaces implicitly
using gray-scale images. For example, an ellipse is represented as
the level set of an image (top). To change the shape of the ellipse
we modify the gray-scale values of the image by solving a PDE
(bottom).
Zhukov et al.
128 / Journal of Electronic Imaging / January 2003 / Vol. 12(1)
Downloaded From: http://astronomicaltelescopes.spiedigitallibrary.org/ on 09/30/2016 Terms of Use: http://spiedigitallibrary.org/ss/termsofuse.aspx