Robust Geometric Primitives

  

 

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:   

• Ignore it - Make as an operating assumption that your program will work correctly only if no three points are collinear, no three lines meet at a point, no intersections happen at the endpoints of line segments, etc. This is probably the most common approach, and what I would recommend for short-term projects if you can live with frequent crashes. The drawback is that interesting data often comes from points sampled on a grid and is inherently very degenerate.
• Fake it - Randomly or symbolically perturb your data so that it seems nondegenerate.   By moving each of your points a small amount in a random direction, you can break many of the existing degeneracies in the data, hopefully without creating too many new problems. This probably should be the first thing to try as soon as you decide that your program is crashing too often. One problem with random perturbations is that they do change the shape of your data in subtle ways, which may be intolerable for your application. There also exist techniques to ``symbolically'' perturb your data to remove degeneracies in a consistent manner, but these require serious study to apply correctly.
• Deal with it - Geometric applications can be made more robust by writing special code to handle each of the special cases that arise. This can work well if done with care at the beginning, but not so well if kludges are added whenever the system crashes. Expect to expend a lot of effort if you are determined to do it right.

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:

• Integer arithmetic - By forcing all points of interest to lie on a fixed-size integer grid, you can perform exact comparisons to test whether any two points are equal or two line segments intersect. The cost is that the intersection point of two lines may not be exactly representable. This is likely to be the simplest and best method, if you can get away with it.  
• Double precision reals - By using double-precision floating point numbers, you may get lucky and avoid numerical errors. Your best bet might be to keep all the data as single-precision reals, and use double-precision for intermediate computations.  
• Arbitrary precision arithmetic - This is certain to be correct, but also to be slow. This approach seems to be gaining favor in the research community with the observation that careful analysis can minimize the need for high-precision arithmetic, and thus the performance penalty. Still, you should expect high-precision arithmetic to be several orders of magnitude slower than standard floating-point arithmetic.  

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:

• Area of a triangle - Although it is well-known that the area A(t) of a triangle t=(a,b,c) is half the base times the height, computing the length of the base and altitude is messy work with trigonometric functions. It is better to use the determinant formula for twice the area:    

  

  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

  

  Note that these are signed volumes and can be negative, so take the absolute value first. Section explains how to compute determinants.

  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].

• Above-below-on test - Does a given point c lie above, below, or on a given line l?   A clean way to deal with this is to represent l as a directed line that passes through point a before point b, and ask whether c lies to the left or right of the directed line l. It is up to you to decide whether left means above or below.

  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).

• Line segment intersection - The above-below primitive can be used to test whether a line intersects a line segment. It does iff one endpoint of the segment is to the left of the line and the other is to the right. Segment intersection is similar but more complicated, and we refer you to implementations described below. The decision whether two segments intersect if they share an endpoint depends upon your application and is representative of the problems of degeneracy.  
• In-circle test - Does the point d lie inside or outside the circle defined by points a, b, and c in the plane? This primitive occurs in all Delaunay triangulation algorithms and can be used as a robust way to do distance comparisons. Assuming that a, b, c are labeled in counterclockwise order around the circle, compute the determinant:  

  

  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 ).