Input description: A point p and a line segment l, or two line segments and .
Problem description: Does p lie over, under, or on l? Does intersect ?
Discussion: Implementing basic geometric primitives is a task fraught with peril, even for such simple tasks as returning the intersection point of two lines. What should you return if the two lines are parallel, meaning they don't intersect at all? What if the lines are identical, so the intersection is not a point but the entire line? What if one of the lines is horizontal, so that in the course of solving the equations for the intersection point you are likely to divide by zero? What if the two lines are almost parallel, so that the intersection point is so far from the origin as to cause arithmetic overflows? These issues become even more complicated for intersecting line segments, since there are a bunch of other special cases that must be watched for and treated specially.
If you are new to implementing geometric algorithms, I suggest that you study O'Rourke's Computational Geometry in C [O'R94] for practical advice and complete implementations of basic geometric algorithms and data structures. You are likely to avoid many headaches by following in his footsteps.
There are two different issues at work here: geometric degeneracy and numerical stability. Degeneracy refers to annoying special cases that must be treated in substantially different ways, such as when two lines intersect in more or less than a single point. There are three primary approaches to dealing with degeneracy:
Geometric computations often involve floating-point arithmetic, which leads to problems with overflows and numerical precision. There are three basic approaches to the issue of numerical stability:
The difficulties associated with producing robust geometric software are still under attack by researchers. The best practical technique is to base your applications on a small set of geometric primitives that handle as much of the low-level geometry as possible. These primitives include:
This formula generalizes to compute d! times the volume of a simplex in d dimensions. Thus 3! = 6 times the volume of a tetrahedron t=(a,b,c,d) in three dimensions is
The conceptually simplest way to compute the area of a polygon (or polyhedron) is to triangulate it and then sum up the area of each triangle. An implementation of a slicker algorithm that avoids triangulation is discussed in [O'R94].
This primitive can be implemented using the sign of the area of a triangle as computed above. If the area of t(a,b,c) > 0, then c lies to the left of . If the area of t(a,b,c) = 0, then c lies on . Finally, if the area of t(a,b,c) < 0, then c lies to the right of . This generalizes naturally to three dimensions, where the sign of the area denotes whether d lies above or below the oriented plane (a,b,c).
Incircle will return 0 if all four points are cocircular, a positive value if d is inside the circle, and negative if d is outside.
Check out the implementations described below before you attempt to build your own.
Implementations: LEDA (see Section ) provides a very complete set of geometric primitives for planar geometry, written in C++. If you are writing a significant geometric application, you should consider basing it on LEDA. At least check them out before you try to write your own.
O'Rourke [O'R94] provides implementations in C of most of the primitives discussed in this section. See Section . These primitives were implemented primarily for exposition rather than production use, but they should be quite reliable and might be more appropriate than LEDA for small applications.
A robust implementation of the basic geometric primitives in C++ using exact arithmetic, by Jonathan Shewchuk, is available at http://www.cs.cmu.edu/ quake/robust.html. Don't expect them to be very fast.
Pascal implementations of basic geometric primitives appear in [MS91]. Sedgewick [Sed92] provides fragments of the basic primitives in C++. See Section for both of them.
An alternative C++ library of geometric algorithms and data structures (although you are almost certainly better off sticking to LEDA) is Geolab, written by Pedro J. de Rezende, Welson R. Jacometti, Cesar N. Gon, and Laerte F. Morgado, Universidade Estadual de Campinas, Brazil. Geolab requires the SUN C++ compiler, but a Sparc binary and visualization environment is included along with all source code. Geolab appears to be primarily for the brave, since its robustness is uncertain and it contains little documentation, but it does provide 40 algorithms, including such advanced topics as farthest point Voronoi diagrams, nearest neighbor search, and ray shooting.
Notes: O'Rourke [O'R94] provides an implementation-oriented introduction to computational geometry, which stresses robust geometric primitives and is recommended reading.
Shewchuk [She96] and Fortune and van Wyk [FvW93] present careful studies on the costs of using arbitrary-precision arithmetic for geometric computation. By being careful about when to use it, reasonable efficiency can be maintained while achieving complete robustness. Other approaches to achieving robustness include [DS88, Hof89, Mil89].
Related Problems: Intersection detection (see page ), maintaining arrangements (see page ).