From: Tom Duff <td@pixar.com>
Newsgroups: comp.graphics.algorithms,sci.math
Subject: Re: Algorithm for convex hull in higher dimensions
Date: Thu, 02 Oct 1997 14:14:55 -0700

Oxford Materials wrote:
> 
> I have not heard of any algorithm for finding the convex hull
> of a set of points which works in higher dimensions, so I have
> created one which is based on the following 2 properties:
> [ad hoc, non-algorithm deleted]

Try quickhull -- it's easy, fast, accurate, versatile and already
written:  http://www.geom.umn.edu/software/qhull/
-- 
Tom Duff, KF6LWB

==============================================================================

From: Jeff Erickson <jeffe@cs.duke.edu>
Newsgroups: comp.graphics.algorithms,sci.math
Subject: Re: Algorithm for convex hull in higher dimensions
Date: Thu, 02 Oct 1997 19:14:24 -0400

Oxford Materials wrote:

> I have not heard of any algorithm for finding the convex hull
> of a set of points which works in higher dimensions

I have. :-)  There are several, of varying complexities and running
times.  For a nice survey, see the chapter "Convex Hull Computations" by
Raimund Seidel in the CRC Handbook of Discrete and Computational
Geometry.  If you want code, see any of the following web pages.  The
first two are specific programs; the last two are lists pointing to
several more programs.

	http://www.geom.umn.edu/software/qhull/
	http://netlib.bell-labs.com/netlib/voronoi/hull.html

	http://www.geom.umn.edu/software/cglist/ch.html
	http://www.cs.duke.edu/~jeffe/compgeom/code.html#topes

You've rediscovered (a variant of) the "quickhull" algorithm -- starting
with some "base" polytope, insert the point farthest from some facet
until all the points are inside.  Except possibly for the initial step,
this is the algorithm implemented in "qhull" (a the first URL above).

> The problem is that the polytope may not always remain convex
> (that is, adjacent faces may be angled towards each other) during
> this algorithm.  Thus the "furthest in the normal direction" point
> may not be associated with the current face.

Also, even if the furthest point from a facet is associated with that
facet, the point might be able to "see" other facets.  Just joining the
facet with the new point won't be enough.

One way to fix the algorithm is to fill in the concavities as soon as
they appear.  After you connect a point to a facet, creating 3 new
triangular facets, check the edges of the facet you threw away.  If an
edge is concave, then you need to "flip" the pair of triangles that
contains it:
	 ---             ---
	|\  |           |  /|
	| \ |   ====>   | / |
	|  \|           |/  |
	 ---             ---
(If you prefer to think of the convex hull as a solid object instead of
a surface, think of gluing a tetrahedron into the groove between the two
facets.)  This deletes two facets and creates two new ones.  Any points
that were associated with the two old facets either become interior
points or are now associated with one of the new facets.

Every time you create a new facet, recursively check any of its edges
that do not touch the point you're trying to add.  After each flip,
there will be two new edges to check, which may lead to more flips.  The
edges you check will expand out from the new point in a wave. 
Eventually this process will delete precisely the old facets that the
new point could see.

If you're familiar with the incremental flipping algorithm for Delaunay
trianglations, bells should be ringing in your head right about now. 
It's not just a similar algorithm; it's the SAME algorithm.

All this stuff works in higher dimensions, too, except that flips are a
little more complicated.

If you add the points in random order, instead of always adding the
"furthest" point, the running time of this algorithm (in 3d) is
O(n log n), which is the best you can do in the worst case.  However,
adding the "furthest" points first is probably going to be even faster
in practice, since you'll quickly throw away most of the internal
points.

All this assumes (implicitly) that the facets are going to be triangles
at every stage of construction.  If your point sets contains four or
more points on the same plane,you could have non-triangular facets,
which could cause you problems.  Even worse, if some four points are
almost on the same plane, your algorithm may think they're coplanar due
to lack of numerical precision, and this almost certianly WILL cause
problems.  The programs I pointed to earlier, qhull and hull, devote a
LOT of attention to this aspect of the computation.

--
Jeff Erickson                         Center for Geometric Computing
jeffe@cs.duke.edu                     Department of Computer Science
http://www.cs.duke.edu/~jeffe                        Duke University