A Fast Marching Method of Computing Closest Points
Sean Mauch and David Breen
California Institute of Technology
Distance, Closest Points and the Eikonal Equation.
Let u(x,y,z) denote the signed distance from the closed surface
S. u satisfies the Eikonal equation,
The characteristics of this equation are straight lines that are normal
to S and point in the direction of increasing distance. For
each point in space, there is a line segment from the surface to the point
that is a characteristic of the entropy-satisfying solution of theEikonal
equation. A point (x,y,z) and the closest point on the surface S
are the endpoints of this line segment.
Previous Work
The Fast Marching Level Set Method.
Sethian [1], [2] has developed a fast marching level set method to solve
the Eikonal equation,
in the case that f is either always positive or negative.The
method uses an upwind, viscosity solution, finite difference scheme to
numerically solve this equation. For f(x,y,z)=1 and g(x,y,z)
= 0, the solution gives the signed distance from the surface S.
The initial condition u=0 on S is specified by giving
the valueof u on a narrow band of points around the surface S.
The distance is computed by pushing this narrow band outward. In
one step of the method:
-
The point with the smallest distance is removed from the narrow band and
it's value is frozen.
-
Points are added to the narrow band to maintain unit thickness.
-
The distance of the neighbors of the removed point are recomputed using
the finite difference scheme.
Complexity and Accuracy.
Recomputing the distance for a point is accomplished by iteratively solvinga
multivariable, quadratic equation that implements the finite difference
scheme. The narrow band is maintained as a priority queue, which allows
one to efficiently add points and find the point in the narrow band with
the smallest distance.If there are T elements in the priority queue,
then the operations of removing the element with the smallest distance
and inserting a new element takes O(log T) operations. Suppose
the distance is computedon an N*N*N grid. At any given point, the
narrow band will contain O(N2) points. Thus the complexity
of Sethian's algorithm is N3 O(log N2)=O(N3 log N). (There is
a factor of N3 because each point in the grid is removed from the
narrow band once.)
The finite difference schemes presented in [1] and [2] are first order.
Thus the error in the computed distance will be of the order of the spacing
on the grid. To increase the accuracy, one must refine the grid.
Motivation
To calculate the closest points to a surface on a regular grid, we utilize
Sethian's fast marching method, but instead of using a finite difference
scheme to compute distance, we use a heuristic algorithm to propogate closest
points information. Instead of specifying the distance for the points
in the narrow band as an initial condition, we specify the closest points
to the surface. In one step of the closest points method:
-
The point g with the smallest distance is removed from the narrow
band and it's value is frozen.
-
Points are added to the narrow band to maintain unit thickness.
-
The closest points of the neighbors with larger distances than g
are recomputed using the closest point information from g.
The closest points method is based on the following idea. The closest
point on the surface to a point in the grid is usually close to one of
the closest points of its neighbors in the grid. Thus if one knows
the closest points of the neighbors of a grid point g, one can compute
an approximate closest point for g by assuming that it is
near one of the closest points of its neighbors. This is only a heuristic,
in the figures below we see cases in two dimensions for which the heuristic
succeeds and fails.In the cases where the heuristic fails to determine
the correct closest point, it still gives a reasonable approximation of
the distance. The heuristic may fail if the characteristics from
several different portions of the surface S intersect near g.
Fortunately, if the heuristic fails at a point, this mistake is usually
not propogated outward to increasing distances. This is because ``information''
in the eikonal equation and the closest points method is propogated along
characteristics. Where characteristics collide, information goes
into the shock and is lost.
Terminology
Let the distance grid be the N*N*N grid that spans the distance
volume. We will refer to points in the distance grid with (i,j,k)
coordinates. Let the zero set grid be an M*M*M uniform
grid that spans the same Cartesian domain. We refer to the ratio
M/N as the super-sampling factor of the zero set grid.
In most cases the zero set grid is finer than the distance grid. We will
refer to points in the zero set grid with (I,J,K) coordinates.
For any grid point, the closest point is defined as the Cartesian coordinates
of the point on the CSG model surface that is closest to that grid point.
Initial Data
The fast marching algorithm takes as input: a set of points in the distance
grid that forms a narrow band around the CSG model surface and a
representation of the surface. The narrow band contains all the points
in the distance grid having the property that a neighbor of the point has
opposite inside/outside status. (In three dimensions neighbor means
one of the 26 locations surrounding each grid point. For the points
in the narrow band we must supply the (i,j,k) coordinates of the
points and their inside/outside status. The narrow band is used as a starting
point for propogating the closest point information outward and inward
to the rest of the points in the distance grid. Note that specifying
the inside/outside status of the points in the narrow band determines the
inside/outside status of the other points in the grid.
We represent the CSG model surface with a set of points that lie on
the surface. This set of points will be called the zero set,
as they are points lying on the isosurface of zero distance. The
zero set is made by first constructing a thin band of points in the zero
set grid that surrounds the CSG model surface. This set of grid points
will be called the zero band. The zero set is the set of closest
points to the grid points in the zero band. Given a point p
in the zero set that is closest to the grid point (i,j,k) in the
zero band, one can determine all the points in the zero set that lie in
a neighborhood of p by determining all the points in the zero band
in a neighborhood of (i,j,k). As input to the algorithm, we
must supply the (I,J,K) coordinates of the grid points in the zero
band and the corresponding (x,y,z) coordinates of the points in
the zero set as input to the algorithm. In the figure below the initial
data is shown pictorially in 2 dimensions.
The Algorithm
Initially, we have the closest points data in the zero band that surrounds
the surface. We use the closest points data in the zero band to determine
the closest points in the narrow band on the distance grid and then march
the narrow band outward and inward to calculate the closest points in the
rest of the distance grid. Consider a point g that neighbors
the band and whose closest point is unknown. The closest point of g
is probably close to one of the closest points of its neighbors in
the band. Thus for each neighbor of g in the band, we recompute
the distance of g by considering points that are near the closest
points of that neighbor. First we will present the marching algorithm
that moves the band outward and then inward. Next we will show the algorithm
for recomputing the distance at a point g, given the closest point
of one of its neighbors.
Let in_out(i,j,k) denote the inside/outside status for a point
in the distance grid; +1 for outside, -1 for inside. Let grid(i,j,k)
denote the computed distance on the distance grid. A value of infinity
indicates that the distance has not yet been computed. Let source(i,j,k)
denote the point in the zero set Z from which this distance was
computed.
Initially: The closest point to each (I,J,K) in the zero band
is known. For each (i,j,k) in the narrow band grid(i,j,k)
= in_out(i,j,k). For each point not in the narrow band gird(i,j,k)
and in_out(i,j,k) are set to be undefined. The closest points
of the zero band are used to generate approximate closest points for the
narrow band. Below is the fast marching, closest points algorithm.
begin
// March forward to find positive distances.
put each point with a non-negative, finite grid(i,j,k) in the
set U;
while U is not empty
remove the grid point g with the smallest
distance from U;
for each of the 26 neighbors of g
if the source of
the neighbor is unknown
add that neighbor to U;
if the distance of
te neighbor is larger than the distance of g
recompute the neighbor's distance using g's source s;
end
Next the narrow band is marched backward to compute the closest points
with negative distance. Below is the algorithm to recompute the distance
grid(i,j,k) to the distance grid point g, using a zero set
source s. Let (I,J,K) be the coordinates in the zero
band for which s is the closest point. The user chooses the
search radius parameter R. This is the radius of a cube around
the point (I,J,K) in the zero band that defines a neighborhood on
the surface around the point s. The parameter, D = 2R+1
is the diameter of the cube. When recomputing the distance, all the
points in the zero set in a neighborhood around s are considered
as possible closest points.
begin
for each grid point (l,m,n) in a D*D*D cube surrounding
(I,J,K)
t in Z is the closest point to (l,m,n);
calculate the distance from g to t;
grid(i,j,k) = minimum of the D3 computed distances;
source(i,j,k) = the source of this minimum distance, (an element
of Z);
end
From experience we have found that for most surfaces, a search radius
R of half the super-sampling factor of the zero set grid will suffice
to achieve the correct closest points information to the set Z.
Finally, note that since the zero band is of small constant thickness,
the number of points in the zero band in the D*D*D cube is O(D2).
Computational Complexity
There are N3 points in the distance grid. Each distance grid
point is removed from the narrow band once, giving us a factor of N3.
At any point in the algorithm, there are O(N2) points in the narrow
band. The cost of adding and deleting elements from the narrow band is
proportional to the logarithm of the number of points in the narrow band.
This gives us a factor of O(log N). The computational cost
of recomputing the distance for a given grid point is proportional to the
number of zero band points in a D*D*D cube neighborhood of a point
s in the zero band. This gives us a factor of O(D2). Thus
the overall computational complexity of the algorithm is O(N3 D2 log
N).
Error Analysis
First note that we can only talk about the error in the computed distance
and not the error in the closest point. Even if the distance were
computed exactly, there would be locations where the closest point would
not be unique, (i.e. locations that are equidistant from two or more points
on the surface).
We make the assumption that the zero set grid is sufficiently fine to
capture the features of the surface. Let d be largest distance
from the surface to the nearest point in the zero set, that is
If the search radius is sufficiently large so that the the closest points
to the zero set are correctly computed, then the maximum error in the computed
distance will be less than d. We will say that the zero set
grid captures the features of the surface if d is of the order of
the zero set grid spacing. For such surfaces the error in the computed
distance is of the order of the zero set grid spacing.
For smooth surfaces, the error decreases with increasing distance from
the surface. This is because the angle between the computed closest
point and the actual closest point becomes smaller with increasing distance.
Comparison to Previous Work
Similarities
-
Both Sethian's level set method and algorithm presented here are fast marching
algorithms. That is, they start with an initial narrow band that
is near the surface of zero distance and then march the narrow band (outward/inward)
to (greater/lesser) distances.
-
The two algorithms have similar computational complexity.
Differences
-
While Sethian's level set method computes distance, the algorithm presented
here computes the closest points, (from which the distance can be calculated).
-
The error in the distance using Sethian's level set method is proportional
to the spacing of the distance grid. The error for the closest point
algorithm is proportional to the spacing of the zero set grid, which is
independent of the distance grid.
-
Unlike the well-developed theory for numerical solutions to differential
equations, it is currently difficult to theoretically analyze the behavior
of the closest points algorithm. Despite the lack of theoretical
guarantees, the closest points algorithm behaves reasonably and predictably
in practice.
Bibliography
[1] J.A. Sethian. A fast marching level set method for monotonically
advancing fronts. In Proc. nat. Acad. Sci., volume 93 of 4,
pages 1591-1595, 1996.
[2] J.A. Sethian Level Set Methods: Evolving interfaces in
geometry, fluid mechanics, computer vision, and materials science.
Cambridge University Press, 1996.
If you have comments or suggestions, email me at sean@cco.caltech.edu