Accurate Computation of the Medial Axis of a
Polyhedron
Tim Culver John Keyser Dinesh Mano cha
DepartmentofComputer Science
University of North Carolina
Chap el Hill, NC 27599-3175
f
culver,keyser,dm
g
@cs.unc.edu
Abstract:
We present an accurate and ecient algorithm to compute the internal
Voronoi region and medial axis of a 3-D p olyhedron. It uses exact arithmetic and
representations for accurate computation of the medial axis. The sheets, seams, and
junctions of the medial axis are represented as trimmed quadric surfaces, algebraic
space curves, and algebraic numbers, respectively. The algorithm works by recursively
nding neighboring junctions along the axis. It utilizes discretization of space and
linear programming to sp eed up the search step. We also present a new algorithm for
analysis of the top ology of an algebraic plane curve, which is the core of our medial axis
algorithm. To speed up the computation, wehave designed specialized algorithms for
fast computation on implicit geometric structures. These include lazy evaluation based
on multivariate Sturm sequences, fast resultant computation, curve top ology analysis,
and oating-p oint lters. The algorithm has b een implemented and we highlight its
performance on a number of examples.
1 Intro duction
The medial axis of a solid ob ject provides useful shape information in terms of geo-
metric proximity of its boundary elements. It is a skeletal representation that can be
formulated as the locus of the center of a maximal sphere as it rolls around the ob ject
interior.
The medial axis of an ob ject was rst prop osed by Blum [Blu67] for biological shap e
measurement. It has also been used for a numb er of other applications, including path
planning, nite element mesh generation, automated injection molding simulation and
feature recognition. Many researchers have proposed the medial axis transform as
an alternate representation to B-rep and CSG in a design and interrogation system
[Wol92, Hof94 , SNTM92].
The medial axis of a solid is the closure of the locus of all p oints within the
solid that havetwo or more closest p oints on the b oundary of the solid. When the
solid is a p olyhedron, the medial axis is comp osed of bisectors of b oundary features
(vertices, edges, and faces). These bisectors are planes and quadric surfaces. The
1
medial axis consists of
sheets
,
seams
, and
junctions
. A sheet is the bisector of two
boundary elements, and may be represented as a trimmed quadric surface. A seam is
an algebraic space curve dened by the intersection of two or more sheets. A junction
point is dened by the intersection of three or more sheets.
The complexity of the medial axis is not fully understo od. The total number of
sheets, seams, and junctions can b e
O
(
n
2
), where
n
is the total number of faces, edges,
and vertices in the p olyhedron; it mayeven b e
O
(
n
3
). Many polyhedral medial axis
algorithms use either a spatial subdivision or some form of discrete representation of
the surface. Sheehy et al. [SAR95] claim that b ecause of the inherent complexity
of a continuous approach, one has to use some form of discretization. However, it is
hard to accurately and eciently compute the medial axis based on such a discrete
representation. In the last few years, a numb er of authors have proposed computing
the medial axis based on a tracing approach [Mil93, Chi92, SPB95 ]. Starting with a
known junction (for instance, a vertex of the polyhedron), they trace along its inci-
dent seams and discover the adjacent junction p oints. In particular, Sherbrookeet
al. [SPB95] have presented the rst implemented algorithm for computing the medial
axis based on a continuous representation of the boundary. However, they compute
numerical solutions of p olynomial equations using oating point arithmetic, and com-
pute piecewise linear approximation of the seam curves. The accuracy varies with the
step size and the tolerances used for computing the junction points. Since the medial
axis of a polyhedron can be very sensitive to perturbations, it can b e rather non-trivial
to design an accurate algorithm using nite precision arithmetic. For example, gure
13 shows the medial axis of a polyhedron which has four junction points very close to
each other. In practice, it is dicult to accurately compute the medial axis of sucha
polyhedron using nite precision arithmetic.
The problem of accurate and robust computation using nite precision arithmetic
is a ma jor op en problem in solid mo deling [Hof89]. One solution to the precision
problem is to use
exact arithmetic
. For example, many algorithms based on exact
arithmetic have been prop osed for reliable computation for boundary evaluation of
bo olean combinations of polyhedral mo dels [SI89, For95 , BMP94, Hof89]. However,
medial axis computation involves low-degree nonlinear algebraic curves and surfaces.
For such problems, it is dicult to compute tight bounds on the numerical error
generated due to oating-point arithmetic. As a result, it is hard to design reliable
algorithms using tolerance-based approaches. At the same time, progress is slowinthe
application of exact arithmetic to nonlinear problems, mainly because exact arithmetic
is perceived to b e impractically slow.
In this paper, we present an accurate and ecient algorithm to compute the in-
ternal Voronoi region and medial axis of a 3-D p olyhedron. Our overall approachis
similar to the 3-D tracing algorithm prop osed by [Mil93 , Chi92 , SPB95]. However, we
use exact arithmetic and representation for accurate computation of the medial axis.
The algorithm recursively computes the vertices and edges of the medial axis based on
seam-tracing. We also present a new curve top ology evaluation algorithm for reliable
tracing of the seams of the medial axis. Our contributions include:
Search algorithm:
An accurate and reliable algorithm for evaluating all the
seams and computing the \rst" junction point lying on the seam.
Representation:
Ecient and exact representations of junction p oints, seams
2
and sheets of the medial axis.
Eciency:
We use techniques based on space discretization and linear pro-
gramming to speed-up the search step and avoid enumerating all possible com-
binations of junction points. To speed up arithmetic operations on algebraic
numbers, wehave designed sp ecialized algorithms for fast evaluation and com-
putation of the junction p oints, seams and sheets. These include lazy evalua-
tion based on multivariate Sturm sequences, fast resultant computations, curve
topology evaluation, and oating-point lters.
Handling Degeneracies:
We identify the cases where degeneracies can aect
our algorithm, and propose ways to identify and resolve them.
The resulting algorithm has been implemented and we highlight its performance on
some example polyhedra. Its overall complexity is output sensitive and depends on
the combinatorial complexity of the medial axis.
Organization:
The rest of the paper is organized in the following manner. We
survey the previous work in medial axis computation and exact arithmetic in section 2.
Section 3 gives an overview of our algorithm and describ es the exact representations
of the elements of the medial axis. In section 4, we describe our new search and
curve top ology analysis algorithms. We presentanumber of techniques to improve
the eciency of our approach in section 5. Section 6 discusses degeneracies in the
medial axis computation and howwe handle them. We describ e our implementation
in Section 7 and highlight its performance on a few p olyhedra.
2 Previous work
2.1 Medial axis computation
There is considerable work on medial axis computation in solid modeling and compu-
tational geometry. The various approaches can b e classied into two main categories:
discrete metho ds based on sampling, and \continuous" or direct methods that do not
rely on the suciency of some sampling.
Sampling p oints on the surface:
Many authors have computed the medial axis
using the Voronoi diagram of a set of p oints lo cated on the polyhedron's surface. Usu-
ally the result is an approximation to the medial axis, though Sheehy et. al. [SAR95]
adaptively rene the sampling until all adjacency relationships are revealed. From
this, the exact medial axis may be constructed. However, the overall complexityof
this approach is not understo od.
Spatial Sub division:
Another popular approach is to imp ose a regular grid on
space, and compute an approximation to the medial axis consisting of a set of cells. The
approximation can be considered an \image" of the medial axis, and does not contain
a complete description of the adjacency relationships. Among implementations of this
approach, that of Vleugels and Overmars [VO95] is notable in that they guarantee
correct top ology of the approximate axis. Their top ology-resolution routine subdivides
adaptively. Thus, like the boundary-discretization method, the complexity is dicult
to formulate in terms of the input size.
2-D Continuous Algorithms:
The classic divide-and-conquer algorithm for the
two-dimensional problem has been presented by Lee [Lee82]. Avariant of that algo-
3
rithm is given in [Hel97]. Fortune has presentedasweepline algorithm [For87 ]. No
algorithms are known in the literature that extend these approaches to 3-D. Imai has
presented an incremental algorithm for 2-D polygons [Ima96].
Incremental 3-D Algorithm:
Milenkovic has presented an an incremental al-
gorithm for computing the voronoi region of a p olyhedron. It inserts the b oundary
elements one at a time. At step
k
, the entire generalized Voronoi diagram of the rst
k
sites has been constructed. However, we are not aware of any implementation of
this algorithm.
3-D Tracing Algorithms:
All the practical algorithms for computing the 3-
D medial axis are based on the
tracing
approach. Starting from a junction point,
a seam emanating from the junction is followed. The seam terminates at another
junction. Once this point is found, the algorithm recursively forks and follows all seams
emanating from that junction. The key step is the search for the seam terminator. In
the two-dimensional medial axis problem, the search is limited to smaller and smaller
subchains of the p olygon. No such property holds for the three-dimensional problem;
each terminating junction is found at the exp ense of a search of (in the worst case) the
entire polyhedron. Thus the complexity is jointly proportional to the input size
n
of
the p olyhedron and the output size
m
of the medial axis. The algorithm's worst-case
behavior is necessarily at least a factor of
n
larger than optimal
O
(
m
) time. A tight
upper bound on
m
, the size of the medial axis, is currently not known. It is known
that
m
=
O
(
n
3+
) for any
>
0 [AAS97], but we know of no example where
m
is
more than
O
(
n
2
).
Milenkovic [Mil93 ] proposed the 3-D tracing algorithm, and discusses its complex-
ity. Chiang [Chi92 ] presented an algorithm for computing the medial axis of a planar
region b ounded by piecewise
C
2
curves. It involves tracing branches using sets of
polynomial equations. Sherbrooke, Patrikalakis and Brisson [SPB95] have presented a
variation on the algorithm. They explicitly trace along the seam, creating a piecewise-
linear approximation to the curve. They applied their algorithm to generate the medial
axis of a number of p olyhedra. Reddy and Turkiyyah [RT95] also presentaversion of
the tracing algorithm, and compare it to an algorithm based on the Voronoi diagram
of a set of p oints on the b oundary.
Medial Axis of CSG Ob jects:
Dutta and Homann [DH90, DH93] and Ho-
mann [Hof94 ] have presented algorithms for computing the medial axis of CSG ob jects.
Their algorithms compute the points of closest approachbetween pairs of boundary
elements. They have also studied exact representation of bisectors arising in medial
axis computation of CSG ob jects bounded by planes, quadrics and torii.
Medial axis to B-rep Conversion:
Anumber of algorithms have been presented
in the literature to convert solids from a representation based on the medial axis to a
boundary representation [Wol92 , Chi92, Bra92, Ver94].
2.2 Exact arithmetic in geometric computation
Exact arithmetic has proven useful in the linear domain. Exact arithmetic means
that each number is determined and stored to whatever precision is necessary. By
allowing arbitrary bit-lengths anyinteger or rational numb er can be computed and
stored exactly. For linear ob jects, rational numb ers are usually all that is required to
perform a geometric operation. In the solid mo deling community, exact arithmetic has
been used successfully for applications involving linear ob jects (e.g. [For95]). Some of
4
the recently develop ed geometric libraries like LEDA and CGAL also provide supp ort
for exact arithmetic.
In the nonlinear domain, however, exact arithmetic becomes more dicult to apply
eectively. This is because computations on non-linear ob jects often require algebraic
numbers, which cannot be explicitly represented bya nite number of bits. Tech-
niques using bit-length estimates may, in the worst case, require bit-lengths which
are exponential with resp ect to the degree of the algebraic functions [Can88, Yu92].
One approach to this problem is to use eld extensions to allow representations of
each new number as it is computed. Many computer algebra systems and symbolic
libraries provide supp ort for eld extensions. However, computing with suchnumbers
can be quite exp ensive. Another approach is to represent an algebraic number as
the unique root of a polynomial (with rational co ecients) within a rational interval.
This is similar to the idea of interval arithmetic. As long as the interval containing
the root can be tightened or cut on demand, many useful queries can be p erformed
exactly using this representation. In one dimension, this can be accomplished through
the use of St urm sequences. In higher dimensions,
multivariate Sturm sequences
are
used. Milne [Mil92 ] and Pedersen [Ped91] have proposed metho ds for such computa-
tions. Keyser et. al. apply these metho ds to the problem of CSG b oundary evaluation
[KKM97]. The authors demonstrate that the geometric op erations for this problem
can b e computed in rational surfaces using only 1-D and 2-D St urm sequences. The
medial axis problem requires the use of either 3-D St urm sequences or eld extensions.
3 Algorithm Overview
Our approach builds on the tracing algorithm prop osed by [Mil93 ], [SPB95 ], and
[RT95]. However, instead of computing a piecewise linear approximation, we exactly
compute all the comp onents of the medial axis. Essentially, the algorithm constructs
the adjacency graph of the medial axis while traversing it. Our terminology is similar
to that of [SPB95].
The input p olyhedron consists of
faces
,
edges
, and
vertices
. The coordinates of the
vertices are represented in terms of rational numbers. The faces, edges, and vertices
are collectively called
boundary elements
. The medial axis consists of
sheets
,
seams
,
and
junctions
. The sheets corresp ond to trimmed quadric surfaces. The seams are
algebraic curves with rational coecients and the junctions are p oints whose co ordi-
nates corresp ond to algebraic numbers (of degree at most 8). A sheet is said to b e
governed
bytwo boundary elements, a seam by three (or more), and a junction by four
(or more). We denote the
squared distance function
associated with a point, line, or
plane by
d
2
e
(
x; y ; z
), where
e
stands for the p oint, line, or plane. The squared distance
function is a quadratic formula in
x; y ;
and
z
.
The algorithm starts by nding a single seam. For example, if there is a vertex
of the p olyhedron with three incident edges, all of which are convex, then its three
incident faces generate a seam. If no vertex of the p olyhedron is trivalent, a more
sophisticated technique is used to get the algorithm o the ground with an initial
seam. This is discussed in section 6.1.
Each remaining b oundary element
e
in the polyhedron is then considered as a
potential quadruple (
e
1
;e
2
;e
3
;e
) generating a junction point on the seam. One seeks
the element
e
4
whose candidate junction
p
1
comes rst along the seam, measured by
5
