c 2006 Gilbert Strang

5.7 Level Sets and the Fast Marching Method

The level sets of f(x, y) are the sets on which the function is constant. For example
f(x, y) = x
2
+ y
2
is constant on circles around the origin. Geometrically, a level plane
z = constant will cut through the surface z = f(x, y) on a level set. One attractive
feature of working with level sets is that their topology can change (pieces of the level
set can separate or come together) just by changing the constant.
Starting from one level set, the signed distance function d(x, y) is especially
important. It gives the distance to the level set, and also the sign: typically d > 0
outside and d < 0 inside. For the unit circle, d = r 1 = x
2
+ y
2
1 will be the
signed distance function. In the mesh generation algorithm of Section 2. , it was
convenient to describe the region by its distance function d(x, y).
A fundamental fact of calculus: The gradient of f(x, y) is perpendicular to its
level sets. Reason: In the tangent direction t to the level set, f(x, y) is not changing
2
and (grad f) t is zero. So grad f is in the normal direction. For the function x
2
+ y ,
the gradient (2x, 2y) points outward from the circular level sets. The gradient of
d(x, y) = x
2
+ y
2
1 points the same way, and it has a special property: The
gradient of a distance function is a unit vector. It is the unit normal n(x, y)
to the level sets. For the circles,
2 2
grad( x
2
+ y
2
1) = (
x
,
y
) and | grad |
2
= +
y
= 1 . (1)

x
r r r
2
r
2
You could think of the level set d(x, y) = 0 as a wall of re. This refront will move
normal to itself. If it has constant velocity 1 then at time T the re will reach all
points on the level set d(x, y) = T .
That wall of re example brings out an important point when the zero level set
has a corner (it might be shaped like a V). The points at distance d outside that set
(the refront at time d) will lie on lines parallel to the sides of the V, and also on a
circular arc of radius d around the corner. For d < 0 the V moves inward. It remains
a V (with no smoothing of the corner).
The central problem of the level set method is to propagate a curve like the
refront. A velocity eld v = (v
1
, v
2
) gives the direction and speed of each point for
the movement. At time t = 0, the curve is the level set where d(x, y) = 0. At later
times the curve is the zero level set of a function (x, y, t). The fundamental level
set equation in its rst form is
d
+ v grad = 0, with = d(x, y) at t = 0 . (2)
dt
In our wall of re example, v would be the unit vector in the normal direction to
the refront: v = n = grad /| grad |. In all cases it is only the normal component
F = v n that moves the curve! Tangential movement (like rotating a circle around
its center) gives no change in the curve as a whole. By rewriting v grad , the level
c 5.7. LEVEL SETS AND THE FAST MARCHING METHOD 2006 Gilbert Strang
set equation takes a second form that is more useful in computation:
grad d
v grad = v | grad | = F | grad | leads to + F | grad | = 0 . (3)
| grad | dt
We only need to know the velocity eld v (and only its normal component F )
near the current location of the level curvenot everywhere else. We are propagating
a curve. The velocity eld may be xed (easiest case) or it may depend on the local
shape of the curve (nonlinear case). An important example is motion by mean
curvature: F = . The neat property | grad | = 1 of distance functions simplies
the formulas for the normal n and curvature :
When is a
n =
grad
becomes n = grad
distance
| grad |
(4)
function
= div n becomes = div(grad ) = Laplacian of
But here is an unfortunate point for t > 0. Constant speed (F = 1) in the normal
direction does maintain the property | grad | = 1 of a distance function. Motion
by mean curvature, and other motions, will destroy this property. To recover the
simple formulas (4) for distance functions, the level set method often reinitializes
the problemrestarting from the current time t
0
and computing the distance function
d(x, y) to the current level set (x, y, t
0
) = 0. This reinitialization was the Fast
Marching Method, which nds distances from nearby meshpoints to the current
level set.
We describe this quick method to compute distances to meshpoints, and then
discuss the numerical solution of the level set equation (3) on the mesh.
Fast Marching Method
The problem is to march outward, computing distances from meshpoints to the in-
terface (the current level set where = 0). Imagine that we know these distances for
the grid points adjacent to the interface. (We describe fast marching but not the full
algorithm of reinitialization.) The key step is to compute the distance to the next
nearest meshpoint. Then the front moves further outward with velocity F = 1. When
the front crosses a new meshpoint, it will become the next nearest and its distance
will be settled next.
So we accept one meshpoint at a time. Distances to further meshpoints are tenta-
tive (not accepted). They have to be recomputed using the newly accepted meshpoint
and its distance. The Fast Marching Method must quickly take these steps recur-
sively:
1. Find the tentative meshpoint p with smallest distance (to be accepted).
2. Update the tentative distances to all meshpoints adjacent to p.
c 2006 Gilbert Strang
To speed up step 1, we maintain a binary tree of unaccepted meshpoints and their
tentative distances. The smallest distance is at the top of the tree, which identies
p. When that value is removed from the tree, others move up to form the new tree.
Recursively, each vacancy is lled by the smaller of the two distance values below
it. Then step 2 updates those values at points adjacent to p. These updated values
may have to move (a little) up or down to reset the tree. In general, the updated
values should be smaller (they mostly move up, since they have the latest meshpoint
p as a new candidate in nding the shortest route to the original interface).
The Fast Marching Method nds distances to N meshpoints in time O(N log N ).
The method applies when the front moves in one direction only. The underlying
equation is F |T | = 1 (Eikonal equation with F > 0). The front never crosses a
point twice (and the crossing time is T ). If the front is allowed to move in both
directions, and F can change sign, we need the initial value formulation (3).
Lagrangian versus Eulerian
A fundamental choice in analyzing and computing uid ow is between Lagrange
and Euler. For the minimizing function in optimization, they arrived at the same
Euler-Lagrange equation. In studying uids, they chose very dierent approaches:
Lagrange follows the path of each particle of uid. He moves.
Euler sees which particles pass through each point. He sits.
Lagrange is more direct. He tracks the front. At time zero, points on the front have
positions x(0). They move according to vector dierential equations dx/dt = V (x). If
we mark and follow a nite set of points, equally spaced at the start, serious diculties
can appear. Their spacing can get very tight or very wide (forcing us to remove or add
marker points). The initial curve can split apart or cross itself (changes of topology).
The level set method escapes from these diculties by going Eulerian.
For Euler, the x-y coordinate system is xed. He captures the front implicitly,
as a level set of (x, y, t). When the computational grid is also xed, we are con-
stantly interpolating to locate level sets and compute distance functions. Squeezing
or stretching or tangling of the front appear as changes in , not as disasters for the
mesh.
The velocity v on the interface determines its movement. When the level set
method needs v at a meshpoint o the interface, a good candidate is the value of v
at the nearest point on the interface.
Upwind Dierencing
The level set nite dierence method is properly developed in the books by its origina-
tors: Sethian [ ] and Osher and Fedkiw [ ]. Here we concentrate on an essential point:
c 5.7. LEVEL SETS AND THE FAST MARCHING METHOD 2006 Gilbert Strang
upwind dierencing. Recall from Section 1. the three simplest approximations
F, B, Cto the rst derivative d/dxForward, Backward, Centered :
(x + h) (x) (x) (x h) (x + h) (x h)
F B C
h h 2h
Which do we use in the simple convection equation d/dt + a d/dx = 0? Its true
solution is (x at, 0). The choice of nite dierences depends on the sign of a.
The ow moves left to right for a < 0. Then the backward dierence is naturalthe
upwind value (xh, t) on the left should contribute to (x, t+t). The downwind
value (x + h, t) on the right moves further downwind during the time step, and has
no inuence at x.
When the movement of the solution (and the wind) is right to left (with a > 0),
then the forward dierence will use the appropriate upwind value (x + h, t) along
with (x, t), in computing the new (x, t + t).
Notice the time-step limitation |a|t h. In time t, the wind will bring
the true value of from x + at to the point x. If a > 0 and nite dierences reach
upwind to x + h, that must be far enough to include information at x + at. So the
Courant-Friedrichs-Lewy condition is at h. The numerical waves must propagate
at least as fast as the physical waves (and in the right direction!). Downwind dier-
encing is looking for that information on the wrong side of the point x, and is doomed
to failure. Centered dierencing in space is unstable for ordinary forward Euler.
By careful choice of the right nite dierences, Osher has constructed higher-order
essentially non-oscillatory (ENO) schemes. A central idea in nonlinear problems,
where the dierential equation has multiple solutions (see Section ), is to choose
the viscosity solution. This physically correct solution appears in the limit as an
extra u
xx
diusion term goes to zero. With good dierencing the viscosity solution
is the one that appears as x 0.
At this point, the level set method does not appear in large production codes.
In research papers it has successfully solved a great variety of dicult nonlinear
problems.