Line Segment Intersection, in depth: classifying every outcome, exact intersection points, polygon validation and a measured grid broad phase

Five answers, not two

Treat each segment as a closed set: its endpoints belong to it. Write the first as P(t) = p + t*d with d = q - p and t in [0, 1], and the second as R(u) = r + u*e with e = s - r. An intersection is any pair (t, u) in the unit square with P(t) = R(u). Depending on the geometry there are five answers, and each one means something different downstream.

Overlay splits segments at a crossing but not at a shared endpoint; a polygon validator accepts a touch between consecutive edges and nothing else; snapping merges an overlap into one edge. None of these can be done with a boolean.

Solving for the parameters

Take the 2D cross product a x b = a.x*b.y - a.y*b.x. Let w = r - p. Then t*d - u*e = w. Crossing both sides with e removes u and gives t = (w x e) / (d x e). Crossing with d removes t and gives u = (w x d) / (d x e). Two divisions by one denominator solve the whole problem when d x e is not zero.

Worked example: p = (0, 0), q = (4, 2), r = (1, 3), s = (3, -1). Then d = (4, 2), e = (2, -4) and w = (1, 3). The denominator d x e = 4*(-4) - 2*2 = -20. Next, w x e = 1*(-4) - 3*2 = -10, so t = 1/2, and w x d = 1*2 - 3*4 = -10, so u = 1/2. Both lie strictly inside (0, 1), so this is a proper crossing at p + d/2 = (2, 1).

When d x e = 0 the segments are parallel. If w x d is also zero, they lie on the same line. Project the second segment onto the first: t0 = (w . d) / (d . d) and t1 = t0 + (e . d) / (d . d). Clamp the interval [min, max] to [0, 1]. An empty interval means disjoint, a single value means they touch, and a longer interval is the overlap. For (0,0)-(4,0) against (2,0)-(6,0) the projection is [1/2, 3/2], which clamps to [1/2, 1], so the overlap is (2,0)-(4,0).

A tested exact classifier

from fractions import Fraction

def classify(p, q, r, s):
    """Classify closed segments pq and rs. Inputs must be exact (int or Fraction)."""
    d = (q[0] - p[0], q[1] - p[1])
    e = (s[0] - r[0], s[1] - r[1])
    if d == (0, 0) or e == (0, 0):
        raise ValueError("zero-length segment; clean the input first")
    w = (r[0] - p[0], r[1] - p[1])
    denom = d[0] * e[1] - d[1] * e[0]                 # d x e
    if denom == 0:
        if w[0] * d[1] - w[1] * d[0] != 0:            # parallel, different lines
            return ("disjoint", None)
        dd = d[0] * d[0] + d[1] * d[1]                # collinear: project rs onto pq
        t0 = Fraction(w[0] * d[0] + w[1] * d[1], dd)
        t1 = t0 + Fraction(e[0] * d[0] + e[1] * d[1], dd)
        lo, hi = max(0, min(t0, t1)), min(1, max(t0, t1))
        if lo > hi:
            return ("disjoint", None)
        a = (p[0] + lo * d[0], p[1] + lo * d[1])
        if lo == hi:
            return ("touch", a)
        return ("overlap", (a, (p[0] + hi * d[0], p[1] + hi * d[1])))
    t = Fraction(w[0] * e[1] - w[1] * e[0], denom)   # (w x e) / (d x e)
    u = Fraction(w[0] * d[1] - w[1] * d[0], denom)   # (w x d) / (d x e)
    if not (0 <= t <= 1 and 0 <= u <= 1):
        return ("disjoint", None)
    pt = (p[0] + t * d[0], p[1] + t * d[1])
    if 0 < t < 1 and 0 < u < 1:
        return ("cross", pt)
    return ("touch", pt)

Zero-length segments are rejected: they are almost always a data error, such as a duplicated vertex. Running the classifier on eight test pairs gave:

Exactness: integers, rationals and floats

With integer inputs, every quantity above is exact. The intersection point is rational, with denominator d x e. With b-bit coordinates, the denominator needs about 2b bits and the point's numerators about 3b. That is why the code keeps Fraction rather than dividing into floats. If you must store the point in the input's integer grid, rounding it moves it off both segments. That is the problem snap rounding solves, and you should use a library such as GEOS or CGAL for it rather than invent one.

Floating-point input is where most production bugs live. To measure it, I took 100,000 random pairs a, b in the unit square, put a third point c = a + t*(b - a) on the line in double precision, and compared the float orientation determinant with the exact rational one. Because c is rounded, exact arithmetic said it was strictly off the line in all 100,000 cases. The float determinant returned exactly zero in 18,629 of them, about 19%. It never returned the opposite sign in that run. In a float classifier, an endpoint misjudged as lying on the other segment's line makes a parameter exactly 0 or 1, so pairs that miss or cross by an ulp come back as a spurious touch or T-junction.

Three fixes work, in order of preference. First, quantise to an integer grid at the data's real precision (nanometres in a chip layout, 1e-7 degrees in GIS). Second, use exact or adaptive predicates; Robust Orientation Predicate shows how to get exact signs from floats cheaply. Third, use an explicit tolerance policy, and accept that epsilon tests are not transitive: a near b and b near c does not give a near c. That breaks any algorithm, including every sweep, that sorts by them.

Using the answer: is this polygon simple?

The most common use of a classifier is checking that a polygon ring is simple. Edges i and i+1 must meet only at their shared vertex. All other pairs must be disjoint. Anything else is a validity error, and GIS libraries report the same conditions as self-intersection errors.

def simple_polygon_errors(poly):
    n = len(poly)
    edges = [(poly[i], poly[(i + 1) % n]) for i in range(n)]
    errs = []
    for i in range(n):
        for j in range(i + 1, n):
            kind, g = classify(*edges[i], *edges[j])
            adjacent = j == i + 1 or (i == 0 and j == n - 1)
            if kind == "disjoint" or (adjacent and kind == "touch"):
                continue                      # the shared vertex is allowed
            errs.append((i, j, kind, g))
    return errs

On four test rings, the square (0,0), (4,0), (4,4), (0,4) returned no errors. The bowtie (0,0), (4,4), (4,0), (0,4) returned one crossing, between edges 0 and 2 at (2, 2). A ring with a spike, (0,0), (4,0), (4,4), (4,2), (0,4), returned an overlap between consecutive edges 1 and 2 along (4,2)-(4,4), plus a touch at (4, 2). That is exactly why the adjacency exemption must allow only a touch, not an overlap. A ring that pinches through (2, 2) twice returned four touches at that vertex. Whether a pinch is valid depends on your data model, and the classifier gives you the information to decide. The double loop is O(n^2), which is fine for rings of a few hundred vertices. For larger rings, feed the edges through the broad phase below.

A broad phase without a sweep

A uniform grid is the simplest broad phase. Put each segment into every cell its bounding box covers, then test only pairs that share a cell. Two segments whose boxes overlap in several cells would be tested several times. The reference-cell rule fixes this without a hash set: test a pair only in the cell that contains the low corner of the overlap of their two boxes. That corner lies in exactly one cell, so each pair is tested once.

from collections import defaultdict

def grid_pairs(segs, cell):
    buckets = defaultdict(list)
    for i, (p, q) in enumerate(segs):
        x0, x1 = sorted((p[0], q[0])); y0, y1 = sorted((p[1], q[1]))
        for cx in range(x0 // cell, x1 // cell + 1):
            for cy in range(y0 // cell, y1 // cell + 1):
                buckets[(cx, cy)].append(i)
    hits = []
    for (cx, cy), ids in buckets.items():
        for a in range(len(ids)):
            p, q = segs[ids[a]]
            for b in range(a + 1, len(ids)):
                r, s = segs[ids[b]]
                ox = max(min(p[0], q[0]), min(r[0], s[0]))
                oy = max(min(p[1], q[1]), min(r[1], s[1]))
                if (ox // cell, oy // cell) != (cx, cy):
                    continue                  # reference-cell rule
                if classify(p, q, r, s)[0] != "disjoint":
                    hits.append((ids[a], ids[b]))
    return hits

Measured: cell size and the long-segment trap

I generated segments with integer endpoints in a square of side U = 2^20. Each segment got a random start and offsets drawn uniformly from [-U/100, U/100], so typical segments are short relative to the scene, as in road networks. Both methods call the same exact classifier in single-threaded CPython, so compare ratios, not absolute times.

Every grid run at n = 2,000 returned exactly the brute-force pair set. Three lessons follow. First, cell size has a sweet spot near the typical segment length: U/400 tested fewer pairs but was twice as slow as U/100, and larger cells drift toward brute force. Second, the right cell size moves with density: from 20,000 to 80,000 segments the candidate count grew 16 times while n grew 4 times, because the cells were holding too many segments. Third, a few long segments wreck the grid. Making 1% of the segments span the scene multiplied the work by six and the time by fifteen. Pull long segments into a separate list and test them against an R-tree of the short ones, or walk only the cells the segment actually crosses with a DDA rather than every cell of its bounding box.

Choosing a method

In practice, use a grid or an R-tree for broad phases on real data, and the sweep when you need a bound that does not depend on the data's distribution. Whichever you choose, the narrow phase is the classifier above, and its exactness decides whether the result is correct. Background on the vector primitives is in 2D Geometry Primitives.

Failure modes

• Boolean-only tests in overlay code split segments at shared endpoints and leave zero-length edges.
• Epsilon comparisons inside sorts or trees corrupt ordering silently.
• Rounding the intersection point back to the input grid without snap rounding. The rounded point can lie on neither segment and creates new crossings.
• An adjacent-edge exemption that allows overlap lets spikes pass validation.
• Fixed grid cell size on growing data. Candidates grow faster than n, as the 80,000-segment row shows.
• Integer overflow. In fixed-width languages, d x e on 32-bit coordinates needs 64 bits, and the point's numerators need more. Python hides this; Java and C++ do not.

What to do next

1. Run classify on the eight pairs above, then on your data's awkward cases: shared vertices, repeated points, collinear chains.
2. Choose and document your number model: an integer grid or exact predicates on floats.
3. Run simple_polygon_errors over a sample of your polygons, and count how many fail and in which category.
4. For bulk checks, size grid cells near the median segment length and split off long segments.
5. If you need worst-case bounds or have many long segments, read the sweep-line article next.