stolfi@jumbo.dec.com (Jorge Stolfi) (05/24/88)
*++-+++-++-++--++++-+-+++-+----+-++----+--++--+-+-+-+--+++--+++-++--+-
Mitch Norcross asks:
I'm trying to find if a point (x,y,z) lies inside or outside of
a closed polyhedron.
You can use the three-dimensional equivalent of Jordan's theorem for
plane curves. Informally, a point p is inside a closed orientable
surface Q if and only if a path from p to infinity crosses the surface
an odd number of times.
The path can be a straight ray from p to a fixed point at infinity,
say (infinity, 0, 0), except that you must worry about degenerate
cases (p on the surface, ray tangent to the surface, ray passing through
a vertex or edge, etc.). To figure out what to do in those cases,
you may think of the point p as having coordinates (x+eps, y+eps^2,
z+eps^3), where eps is an infinitesimally small quantity (positive
but smaller than any positive real number). This trick will move
p out of the surface, and will prevent the ray from passing through
edges or vertices. The coordinates of the interesection points will
be polynomials in eps instead of real numbers, but these are easy
to compute and compare (note that a > b*eps > c*eps^2 + d*eps^3, for
any positive reals a,b,c,d). (If this sounds too cryptic, I will galdly
elaborate.)
This algorithm takes time proportional to the complexity of the
polyhedron, and works for any orientable closed surface, even if it is
non-convex, with holes and handles, or with multiple separate components.
If the polyhedron is convex, the following "walking" algorithm is
usually faster. Pick a reference vertex A and a starting face F that
is not incident to A. Test whether p is inside the pyramid with vertex
A and base F; this requires testing p against the plane of F and the
planes determined by A and every two consecutive vertices of F). If
p passes all tests, it is inside the polyhedron; return "YES". If
p is on the wrong side of F's plane, it is outside the polyhedron;
return "NO". Finally, if p is on the wrong side of the plane ABC,
where BC is an edge of F, take the face G on the other side of edge
BC; if G is incident to A return "NO", else set F := G and repeat.
Again, the hard part is taking care of boundary cases; the eps-trick
above should work here, too. For suitably random points and polyhedra,
this algorithm will stop in time proportional to the square root of
the polyhedron's size.
If that is not fast enough, you may consider some algorithm based on
dividing space into cells. But at this point things start being too
much application-dependent, and I don't feel like going into details.
For a convex polyhedron it is theoretically possible to do the test in
time proportional to the logarithm of the polyhedron size: project
everything stereographicaly from a vertex and use an optimal point
location algorithm. If this doesn't mean much to you, don't worry;
those algorithms aren't practical anyway.
How about an algorithm for finding a point which lies inside
the polyhedron?
If the polyhedron is convex, take the barycenter (average) of any
four distinct vertices.
If it is not, pick a face F. Pick a point q inside F (recursively
apply the same algorithm, with one dimension less). Pick a ray from q
to infinity that is not coplanar with F and is directed inwards relative
to the orientation of F. Find all intersections between this ray
and the polyhedron's faces. Find the intersection p that is nearest
to q. Return (p+q)/2.
Hope this helps,
Jorge Stolfi
stolfi@src.dec.com, ..{..decvax,ucbvax,allegra}!decwrl!stolfi
DISCLAIMER: Digital Equipment Corporation shall not be held legally
or morally responsible if the above turns out to be a pile of garbage.
*-+-++++---+-++-+-+-++---++-++---+-++--+++-++++--++-+++++---++-+-+-+++bs@linus.UUCP (Robert D. Silverman) (05/26/88)
In article <13077@jumbo.dec.com> stolfi@src.dec.com (Jorge Stolfi) writes:
:*++-+++-++-++--++++-+-+++-+----+-++----+--++--+-+-+-+--+++--+++-++--+-
:
:Mitch Norcross asks:
:
: I'm trying to find if a point (x,y,z) lies inside or outside of
: a closed polyhedron.
:
:You can use the three-dimensional equivalent of Jordan's theorem for
:plane curves. Informally, a point p is inside a closed orientable
etc.
A FAST way to do this numerically is to treat it as a linear programming
problem and simply determine whether the point is a feasible point for
the problem by reducing the constraint matrix to row-echelon form.
the constraint matrix is construced from the planes defining the polyhedron.
Bob Silverman