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 Thursday 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.
Thanks, guys.
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 singlecondclause: 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.
It’s much easier in projective geometry.
Cartesian coordinates of points to projective points:
(x,y) => (x,y,1)
The projective line ax+by+cz=0 through two projective points:
(y0z1-y1z0, z0x1-z1x0, x0y1-x1y0)
The projective point where two lines intersect:
same equation as the previous one, with the variable names changed
Projective back to Cartesian:
(x,y,z) => (x/z,y/z)
If you get a divide by zero in the last step, the lines are parallel.
I would find the intersection point of the two entire lines (not segments), then check the distance from the midpoint of each segment to the intersection point. If the distance is both cases is less than 1/2 the segment length, then the segments intersect.
(I think that would work. Didn’t test it though.)
Also, if you leave the cartesian coordinate system and use polar coordinates, you can specify a line in it’s normal form, as is done with the Hough transform:
http://planetmath.org/encyclopedia/HoughTransform.html
In normal form, a line is defined as:
x cos(theta) + y sin(theta) = rho.
That will eliminate the problem of infinite slopes.
Just a thought, since you have two line segments, can you figure out whether the endpoints of one are both in front of or behind the other? If they aren’t, then you know they don’t cross.
I guess what I’m suggesting is that maybe first you can figure out whether one line segment crosses the line made by extending the other, as an easier problem?
My other thought would be to project the one line onto the other, and see if the projected… Uh, I’m not sure where I’m going with that. Wait, if the projection doesn’t overlap (which is easy to check), then the line segments don’t cross?
I don’t really know if that gets you any closer, but those would be the first two things I would try.
I sympathize with your plight: I’m surprised though that you didn’t encounter yet another problem: precision !! After all, how do you really know that the two endpoints of a segment coincide. If they are floats, and came from some other application upstream, it’s not clear that x == y suffices as a check.
Well, projective geometry is easier than Euclidean geometry in part because it has no good way to define “between”, so you can’t really even talk about line segments. (This is similar to the reason that algebra is easier over the complex numbers than over the real numbers; in fact, maybe it’s the exact same reason.) (Also, even the projective geometry is not so easy as all that; the functions from two points to a line, and from two lines to a point, only work if the points/lines are distinct.)
I’m going to try for an algorithm without many case splits. This is inspired by your description of the algorithm in Sedgewick, but it does give you the intersection point.
I’ll use some special notation: “.” is the dot-product operator; if A is a point then Ax is the x coordinate; and if A,B,C are points then Area(ABC) is the signed parallelogram area: (Ax-Bx)*(Cy-By) - (Ay-By)*(Cx-Bx) (then |Area(ABC)| is the area of a parallelogram where 3 of the vertices are A,B,C; and Area(ABC) = -Area(CBA)).
We’re testing whether the line segments AB and CD intersect.
First, if A=B and C=D, then we actually have two points; check whether A=C. (first case split)
Now, assume without loss of generality that A!=B. (second case split)
Next, compute Area(CAB) and Area(BAD). If the product of these numbers is negative, then C and D are on the same side of AB; report no intersection, and exit. (third case split)
If Area(CAB) and Area(BAD) are both 0, then A,B,C,D are all collinear; jump to the CollinearIntersection(AB,CD) routine below. (fourth case split)
Now we know that the lines AB and CD are neither identical nor parallel, so we can find a unique intersection point E. (Note that I said “lines”, not “line segments”; we still don’t know if the line segments intersect.) This intersection point is:
E = C + (D-C)*Area(CAB)/(Area(CAB) + Area(BAD))
(Due to the above case splits, we know that we are not dividing by 0 here.)
Now, A,B,E are collinear; jump to the CollinearIntersection(AB,EE) routine immediately below.
Now we are testing CollinearIntersection(AB,EF) (where EF is either CD or EE, depending on how we got here); we know that A!=B, and that A,B,E,F are all collinear.
Define UNIT=(B-A).(B-A).
Compute dE = (E-A).(B-A)/UNIT, and dF = (F-A).(B-A)/UNIT.
If dE and dF are both less than 0, or both greater than 1, then report no intersection, and exit. (fifth case split)
Now, sort the 4 numbers 0,1,dE,dF; remove the least and the greatest number; and call the remaining middle two numbers dG and dH. (several more case splits, depending on how you count)
Compute G = A + dG*(B-A) and H = A + dH*(B-A). The intersection is the line segment GH (which may be a single point).
Finished! After 5 case splits (into 6 cases), plus a 4-element sort (which you can decide to count as several case splits, or not, as you please).
Of course, in the real world, there’s more to worry about geometric algorithms than just this. Are you computing over the rationals (in which case this algorithm should work perfectly), or over floating-point numbers? If you are computing over floating-point numbers, and you pass in coordinates which are collinear, is it acceptable to sometimes say that there is no intersection, or that there is a single-point intersection, rather than that there is a line-segment intersection? (If you want an exact answer using floating-point arithmetic, then you’ve got a lot more work ahead of you!)
For one thing: what Eppstein said.
For another. Why not go parametric? You represent your lines as a pair of functions R->R. To make it work with being segments, let them be represented as functions [0,1]->R .Thus, if a line goes from (x1,y1) to (x2,y2), then this is represented by the pair
x: t -> x1+t*(x2-x1)
y: t -> y1+t*(y2-y1)
Now, the vertical line has t -> x1, a constant function. The point has both functions constant. And finding the intersection of the lines described by (x1(t),y1(t)), (x2(t),y2(t)), respectively, is the matter of finding t1 and t2 such that x1(t1)=x2(t2) and y1(t1)=y2(t2). This yields a system of linear equations, easily solved, and with a final check for whether the solution stays within the constraints given.
The method is not hard to adapt to rays and lines - just adjust permissible solutions to be positive or all of them, respectively.
To elaborate a little on my answer, so that it applies to line segment intersection and not just line intersection (and yes, line segments are perfectly well defined in projective geometry, only there are two of them between every pair of points so you have to be clear which one you mean):
- If the calculation for the line between one of the pairs of endpoints comes up (0,0,0), the two points are equal. In this case, there is only an intersection if one of these points lies on the line through the other two points, which can be tested by taking the dot product of the point’s coordinates with the line’s.
- If the calculation for the point where the two lines meet comes up (0,0,0), the two lines are coincident. In this case you can test for an intersection by comparing the order of the points’ coordinates.
- As I said in my original comment, if you get a divide by zero converting back to Cartesian coordinates, the lines are parallel. In this case there can be no intersection.
- Finally, once you have the intersection point in Cartesian coordinates, you can test the ordering of its coordinates against the original points to determine whether it lies between them.
But, if you want to perform exact comparisons such as this (do you really have (0,0,0) or just a vector very close to it?) you probably need to be using multiprecision integers, and avoiding floating point, or your results will not be reliable. Or, to put it another way: if you think approximate calculation with floating point is good enough, why are you worrying about cases like two coincident lines that only make sense for exactly defined numbers?
(Oops, my last reply got corrupted by less-than-signs reinterpreted as HTML. Here’s a repost.)
I wrote some lecture notes that include the line intersection formula, as well as a bunch of other geometric constructors. You can download it from
http://www.cs.berkeley.edu/~jrs/meshpapers/robnotes.ps.gz
(in gzipped PostScript form; see Section 3.6). I believe this is the simplest possible formula. The notes also contain some tips on getting a numerically stable answer despite floating-point roundoff.
The notes don’t address how to tell whether the segments intersect, so I’ll address that here (and I should add it to the notes next week). If your segments are ab and cd, and they’re not parallel (see the notes for how to tell that), then they intersect if and only if
Orient2D(c,d,a) x Orient2D(c,d,b) <= 0 AND
Orient2D(a,b,c) x Orient2D(a,b,d) <= 0.
The function Orient2D() is described in the notes, and I have C code for computing it exactly from floating-point input at
http://www.cs.cmu.edu/~quake/robust.html
If the lines are identical, then you determine the segment intersection by comparing x-coordinates (y-coordinates if the lines are vertical).
You might also be interested in the final chapter of the recent O’Reilly book Beautiful Code: Leading Programmers Explain How They Think, which focuses on the related problem of determining whether three points are collinear. The description of both the problem and the development of the eventual solution is lucid, engaging, and insightful. Still, I doubt there’s much you could learn from the author.
David: Thanks for the suggestion. I’ll see if I can learn anything from that author. (I’m an autodidact, but I had a very good teacher.)