Euclid famously said, “There is no royal road to geometry.” Among the non-royal roads, the computational pathway is notably muddy, rutty and potholed. Over the weekend I needed to write a program for a simple geometric task—finding the intersection of two line segments in the plane. It’s
Wednesday now, and I finally have my program. Let me tell you some of my adventures along the way.
First steps. This is a low-key, off-the-cuff, back-of-the-envelope project. The program doesn’t have to be blazingly fast; it isn’t going to be doing real-time animation in a video game. Nor is it going to trace the trajectories of aircraft in the holding pattern at JFK; no lives will be lost if I make a goof. Still, I would like to get correct answers, and not waste cpu cycles too extravagantly.
In Euclid’s way of doing geometry—with straightedge and compass—there’s not much to be said about algorithms for finding intersections. You construct the lines and see where they cross. Most computers, however, are not adept with straightedge and compass; the problem has to be encoded somehow in a more-algebraic language. Descartes knew how to do this. The cartesian equivalent of the Euclidean procedure goes something like this: Given two segments specified by the x and y coordinates of their end points, write the equation y = mx + b for each of the infinite lines that pass through these pairs of points. Then try to solve the two simultaneous equations to find a point (x,y) that lies on both lines; this point, if it exists, must be the intersection of the lines. Now go back to check whether the intersection lies within each of the segments.
Seems straightforward—but watch out for those potholes. For starters, in the equation y = mx + b the slope m is defined as Δy/ Δx. What if the segment is vertical? Then Δx = 0, and the slope is undefined. More insidiously, what if the two end points of a segment are actually the same point, so that the segment has zero length? Again the slope is undefined; now it’s 0/0. Then there’s the matter of parallel line segments. In Euclid’s world, parallel lines just don’t intersect; that’s pretty much how “parallel” is defined. But in dealing with segments, it does makes sense to ask about the intersection of two segments lying along the same line. In this case the intersection could be empty, or it could be a single point of overlap, or it could be a segment. The method of slopes and intercepts and simultaneous equations won’t work in this situation.
Down the garden path. I like a good puzzle, and I try not to cheat, but sometimes it’s just crazy to pretend you’re the first person on Earth who ever faced a problem. So I did some scouting around on the Net, as well as in filing cabinets and on bookshelves.
If you do a Web search for “segment intersection algorithm,” you’ll find lots of references to works by eminent computer scientists: Jon Louis Bentley, Bernard Chazelle, Herbert Edelsbrunner, and others. Following those leads, you’ll find some eminently clever algorithms, the result of a long series of inventions and refinements in a sort of algorithmic arms race that has gone on for more than 20 years. Unfortunately, these powerful tools solve the wrong problem—or at least they don’t solve my problem. They deal with the task of efficiently identifying all the intersections among a large set of line segments; the subtask of finding the specific point of intersection between two given segments is a minor detail generally left as an exercise for the reader.
Furthermore, most of the papers describing these algorithms take a fairly cavalier attitude to the various singularities and degeneracies mentioned above. I quote from Chazelle and Edelsbrunner 1992 (An optimal algorithm for intersecting line segments in the plane, Journal of the Association for Computing Machinery 39:1–54):
For the ease of exposition, we shall assume that no two endpoints have the same x or y coordinates. This, in particular, applies to the two endpoints of the same segment, and thus rules out the presence of vertical or horizontal segments…. Our rationale is that the key ideas of the algorithm are best explained without having to worry about special cases every step of the way. Relaxing the assumptions is very easy (no new ideas are required) but tedious. That’s for the theory. Implementing the algorithm so that the program works in all cases, however, is a daunting task. There are also numerical problems that alone would justify writing another paper. Following a venerable tradition, however, we shall try not to worry too much about it.
Out on the Web I did find a few odd bits of code that address the specifics of the two-segments problem. One author finesses the inconvenience of vertically oriented segments by setting Δy/0 equal to 1010, with the comment, “close enough to infinity.” I wasn’t seriously tempted to follow this example. It’s not that I want a better approximation to infinity; I’d agree that 1010 is probably close enough. But this is one of those cases where infinity itself is not a very good approximation to the result of dividing by zero. Even if you believe that parallel lines meet at infinity, single isolated points don’t get there even in the limit. Also, despite venerable tradition, I do worry that a finite infinity may invite numerical trouble somewhere down the road.
The road not taken. What’s most annoying about the whole vertical-line mess is that it’s not a fundamental geometric constraint but just an artifact of how lines are represented. It’s only by convention that we measure slope with respect to the y axis; we could choose another orientation. Which suggests a way around the problem. If one of the segments turns out to be vertical, rotate the whole coordinate frame before attempting to test for intersections, and then rotate it back afterwards. The rotation is just a matrix multiplication, so it’s not a big computational burden, and you can do it only when needed. I seriously considered this strategy, and even now I wonder if it isn’t the right way to go about it. I finally decided against it because it doesn’t solve the problem of a segment that’s a single point: That kind of degenerate line has no well-defined slope in any coordinate frame.
Perhaps this fact is an argument for outlawing single-point segments altogether. I considered that too, but it seemed a bit pusillanimous, legislating a problem out of existence just because I found it inconvenient to solve. Mathematically, a one-point line segment may be pathological, but computationally it’s just an instance of aliasing—a very common occurrence, where two things turn out to be the same thing. Besides, it would be handy to have an intersection routine that can also test whether a given point is on a line.
On the road again. Lacking a brilliant Gordian-knot insight that would solve the problem with a single stroke, I had no choice but to push on into the tangled maze of if-then-else case analysis. Is either segment vertical? Is either segment a single point? Are the segments parallel (i.e., identical slopes)? If they are parallel, are they also collinear (identical y intercepts)? But be careful: If they’re parallel, they could both be vertical, with undefined y intercepts. Each of these cases could require a separate calculation of the intersection point. I’m not even sure how many cases there are—at least a dozen, but it depends on how you count.
The procedure I eventually wrote (after exploring several more blind alleys) is not quite as ugly as I feared it would be, but I still have the firm sense that there’s a better way. If not a royal road, then at least a trail that doesn’t require four-wheel drive. I was able to get the number of cases down to five, though only by exploiting a little Lisp hocus-pocus. For convenience of reference in what follows I label the segments red and blue.
- Case 1. Both slopes exist, and they are not equal. This is the general case of two lines that are not vertical and not parallel. We can solve the simultaneous equations.
- Case 2. The red slope is well-defined but the blue slope does not exist. The blue segment could be either a vertical segment or a single point. In either case, if the intersection exists, its x coordinate must be that of the blue segment.
- Case 3. The blue slope is well-defined but the red slope does not exist. Ibid, mutatis mutandis.
- Case 4. Both segments have identical slopes and also identical y intercepts; hence they are parallel and collinear. This is where the hocus-pocus comes in: The Lisp predicate
(equalp (slope red-segment) (slope blue-segment))
returns true if both slopes are numeric and are equal, and it also yields true if both slopes are
nil, the Lisp idiom for things that don’t exist. By this trickery, with a similar test on the y intercepts, I pack several cases into a single
condclause: two collinear vertical segments, two collinear non-vertical segments, and various combinations of vertical segments and single points.
- Case 5. None of the above. The only possibility remaining—unless I am mistaken—is equal slopes and different y intercepts: The lines are parallel but not collinear, and so there is no intersection.
All that analysis merely determines whether or not the lines intersect; we still have to check whether or not the intersection lies within the segments. That involves still more case analysis, since the process is a little different for parallel segments than for others. I found a way of doing it that’s not grotesquely ugly.
The road back. Belatedly, after I had a working procedure, I did some more scouting in the literature and discovered a few paths worth exploring.
Robert Sedgewick’s Algorithms textbook suggests a quite different and ingenious method of detecting intersections of segments. Suppose you walk along the red segment, and when you come to the end, you have to turn counterclockwise to reach one end point of the blue segment, and clockwise to reach the other blue end point. Now try the same experiment walking on the blue segment; if again you must turn counterclockwise in one case and clockwise in the other to see the red end points, then the two segments intersect. Conversely, if in either case both end points are reached by turning in the same direction, there can be no intersection. (To be comprehensive, we also have to consider the case where an end point is straight ahead or directly behind.) Nifty! Unfortunately, the procedure only detects the existence of an intersection; I see no easy way to make it yield the coordinates.
Joseph O’Rourke’s Computational Geometry text (I looked at the C version) also has a discussion of intersection-finding. O’Rourke suggests working with parametric equations for the two segments—defining x and y as functions of distance along the segment. This avoids the problem with undefined slopes but encounters another singularity with parallel lines. Overall, the computation does not appear to be notably simpler.
Stuck in the mud. The animus behind this entire rant is a feeling that I must be missing something obvious, that a problem like finding the intersection of two line segments shouldn’t be this hard. The difficulty I’m talking about is not computational but conceptual. There are lots of hard computational problems—graph coloring, say, or factoring integers—for which one can write a very tidy and perspicuous program. True, the program may have a running time that exceeds your lifespan, but it’s easy to describe what needs to be done. Programs for geometric problems, in contrast, seem often to be efficient but hideous, with tangled logic, an abundance of special cases, and hidden numerical perils. Why is that? Nature seems to have no trouble at all detecting intersections or collisions. If two wires cross on a circuit board, you can count on blowing a fuse no matter what the slopes of the conductors. Why can’t we compute the same result so effortlessly?
Or maybe I have indeed missed something obvious—or even something subtle. If anyone has a better intersection algorithm, please pass it on.
Update 2006-09-15. My thanks to all of the readers who so promptly weighed in with good ideas. Over the weekend I’m going to try turning a couple of those suggestions into working code, and I’ll report back.