Segment Intersection, in depth: Shamos-Hoey detection, Bentley-Ottmann reporting, degeneracies and exact arithmetic

The pairwise test

Everything rests on the orientation predicate: the sign of the cross product (b - a) x (c - a) says whether c is left of, right of, or on the line through a and b. Two closed segments intersect if each one's endpoints lie on opposite sides of the other's line, or if a collinear endpoint lies within the other segment's bounding box. With integer coordinates this is exact:

def orient(a, b, c):
    v = (b[0] - a[0]) * (c[1] - a[1]) - (b[1] - a[1]) * (c[0] - a[0])
    return (v > 0) - (v < 0)

def on_box(a, b, p):
    return (min(a[0], b[0]) <= p[0] <= max(a[0], b[0])
            and min(a[1], b[1]) <= p[1] <= max(a[1], b[1]))

def intersects(s, t):
    """Closed segments: shared endpoints and collinear overlaps count."""
    (a, b), (c, d) = s, t
    o1, o2, o3, o4 = orient(a, b, c), orient(a, b, d), orient(c, d, a), orient(c, d, b)
    if o1 != o2 and o3 != o4:
        return True
    return ((o1 == 0 and on_box(a, b, c)) or (o2 == 0 and on_box(a, b, d))
            or (o3 == 0 and on_box(c, d, a)) or (o4 == 0 and on_box(c, d, b)))

With floating-point coordinates the sign can come out wrong for nearly collinear points, and a sweep built on inconsistent signs can corrupt its own data structure. Robust Orientation Predicate shows how to get exact signs from floats cheaply.

Why a sweep works

Imagine a vertical line sweeping left to right. At any position, the segments it crosses can be ordered bottom to top by their y value at the line. This ordered set is the status. It only changes at a few events: when the line reaches a left endpoint (insert), a right endpoint (delete) or an intersection point (two segments swap order).

The key lemma: if segments s and t cross at a point p, and no other intersection lies to the left of p, then just before the sweep reaches p, s and t are adjacent in the status. Nothing can sit between them without crossing one of them first. So it is enough to test pairs that become adjacent, and a pair becomes adjacent only at an event. Each event creates at most two new adjacencies, so the work is a constant number of pair tests per event instead of n tests per segment.

Shamos-Hoey: does any pair intersect?

Shamos-Hoey answers the decision question. It processes only endpoint events, 2n of them, and stops at the first intersection found. Because it stops before the sweep could pass any intersection, the status never needs to handle a swap. At an insert it tests the new segment against its neighbours; at a delete it tests the two segments that become neighbours. The version below has been tested against a brute-force oracle on 20,000 random inputs, including vertical segments, shared endpoints and collinear overlaps:

from bisect import bisect_left
from fractions import Fraction

def y_at(seg, x):
    (x1, y1), (x2, y2) = seg
    if x1 == x2:                                   # vertical: use its lower end
        return Fraction(min(y1, y2))
    return Fraction(y1) + Fraction(y2 - y1, x2 - x1) * (x - x1)

def any_intersection(segs):
    """Return one intersecting pair (i, j), or None."""
    segs = [tuple(sorted(s)) for s in segs]          # left endpoint first
    events = []
    for i, (p, q) in enumerate(segs):
        events.append((p[0], 0, i))                  # inserts sort before deletes at equal x
        events.append((q[0], 1, i))
    events.sort()
    status = []                                      # segment ids, bottom to top at sweep x
    for x, kind, i in events:
        if kind == 0:
            ys = [y_at(segs[j], x) for j in status]
            k = bisect_left(ys, y_at(segs[i], x))
            for j in status[max(k - 1, 0):k + 1]:    # neighbours below and above
                if intersects(segs[i], segs[j]):
                    return (j, i)
            status.insert(k, i)
        else:
            k = status.index(i)
            if 0 < k < len(status) - 1 and intersects(segs[status[k - 1]], segs[status[k + 1]]):
                return (status[k - 1], status[k + 1])
            status.pop(k)
    return None

Two choices make degenerate cases work. Processing inserts before deletes at the same x means segments that only touch at an endpoint are both in the status when they meet. Treating a vertical segment by its lower end, and testing the entry directly above its insertion point, catches any segment it crosses. For clarity this version rebuilds the y keys on every insert, which is O(n); a balanced search tree such as an AVL tree with a comparator that evaluates y at the current sweep x brings each event to O(log n) and the total to O(n log n).

Bentley-Ottmann: report every intersection

Bentley-Ottmann reports every intersection. It keeps the same status but adds intersection points to the event queue as they are discovered, and at each one swaps the two segments and tests their new outer neighbours:

Q = priority queue of events ordered by (x, y); add all endpoints
T = balanced BST of segments ordered by y at the sweep line
while Q not empty:
    p = Q.pop_min()
    U = segments whose left end is p;  L = segments whose right end is p
    C = segments in T that contain p in their interior
    if |U| + |L| + |C| > 1: report p with U, L, C
    remove L and C from T
    insert U and C into T, ordered just to the right of p   # reverses C
    if U and C are empty:
        test the new neighbours s_below and s_above of p; add their crossing right of p to Q
    else:
        test the lowest of U and C against its lower neighbour, and the highest against its upper
        neighbour; add any crossing right of p (or directly above p) to Q

There are 2n + k events, each costing O(log n) in the queue and the tree, so the total is O((n + k) log n) time. A careful implementation keeps only intersection events between currently adjacent pairs in the queue, which bounds memory at O(n). Grouping all segments through one point into a single event, as above, is what makes many segments meeting at a point a single report instead of a cascade of pairwise swaps.

The bound is output-sensitive, which cuts both ways. When k is small it is close to sorting. When k approaches n^2, as in a dense grid of long lines, the log factor makes it slower than brute force. In that regime, just test all pairs.

Worked example

Take four integer segments: s0 from (0, 0) to (10, 2), s1 from (1, 6) to (9, 7), s2 from (2, 1) to (8, 6) and s3 from (3, 5) to (7, 2). The events, sorted by x, are inserts at x = 0, 1, 2, 3 and deletes at x = 7, 8, 9, 10.

1. x = 0: insert s0. The status is [s0].
2. x = 1: s1 is at y = 6 and s0 is at y = 0.2, so s1 goes above. Test s1 against s0: no intersection.
3. x = 2: s2 starts at y = 1, between s0 (y = 0.4) and s1 (y = 6.125). Test s2 against both: no intersection. The status is [s0, s2, s1].
4. x = 3: s3 starts at y = 5. s2 is at about 1.83 and s1 at 6.25, so s3 goes between them. Testing s3 against s2 finds an intersection, and the algorithm returns (s2, s3).

The crossing itself is at (5, 3.5): solving 1 + 5(x - 2)/6 = 5 - 3(x - 3)/4 gives x = 5. The decision is made at x = 3, when the pair first became adjacent. Bentley-Ottmann would instead queue the point (5, 3.5), continue, swap s2 and s3 there and test s3 against s0 and s2 against s1, finding nothing more.

Degenerate cases

Decide up front whether intersections at shared endpoints count. For polygon simplicity they must not count between consecutive edges, but must count between any other pair. For wiring checks, touching is usually a violation.

Numeric precision

With integer coordinates bounded by 2^b in magnitude, the orientation determinant needs about 2b + 2 bits, so 64-bit integers are exact for coordinates up to about 2^30. Bentley-Ottmann is harder: intersection points are rational, and comparing segments at an intersection's x-coordinate needs several times as many bits as the input coordinates. Use exact rationals (Python's Fraction, or a multiprecision library), filtered predicates that use floats when the sign is clearly safe, or snap rounding, which rounds intersection points to a grid while preserving topology. Libraries such as CGAL and GEOS already do this; reaching for them is often the right production choice.

Practical alternatives and trade-offs

Sweeps are not the only option, and for typical data they are not always fastest. Two simple methods report all pairs and are easy to get right:

def all_pairs_prune(segs):
    """Sort by left x; test each segment only against those still overlapping in x."""
    xlo = [min(p[0], q[0]) for p, q in segs]
    xhi = [max(p[0], q[0]) for p, q in segs]
    active, out = [], []
    for i in sorted(range(len(segs)), key=xlo.__getitem__):
        active = [j for j in active if xhi[j] >= xlo[i]]
        out += [(min(i, j), max(i, j)) for j in active if intersects(segs[i], segs[j])]
        active.append(i)
    return out

This sort-and-prune pass costs O(n log n) plus a test for each pair that overlaps in x, which is small when segments are short relative to the scene, as in road networks and meshes. A uniform grid that buckets each segment into the cells it crosses does the same in two dimensions. Both degrade on long segments; both also parallelise trivially, which the sweep does not. For rectangle-shaped queries, Segment Trees for Geometry covers the related sweep with cover counts.

Failure modes

• Comparator inconsistency. A floating-point comparator that says a < b and b < a at different moments corrupts a balanced tree silently.
• Missing the swap. In Bentley-Ottmann, forgetting to test the new outer neighbours after a swap loses intersections far to the right.
• Duplicate events. The same pair can become adjacent several times; deduplicate intersection events or you report a crossing more than once.
• Wrong question. Many uses only need yes or no; Shamos-Hoey is simpler and stops early.
• No oracle. Always test against the O(n^2) brute force on small random inputs with tiny coordinate ranges, which generate degenerate cases quickly.

What to do next

1. Decide the question: any intersection, all pairs, or all points, and whether touching counts.
2. Implement the exact pairwise test with integer or exact coordinates and a brute-force oracle.
3. If you only need a yes or no, implement Shamos-Hoey and stress-test it with small coordinate ranges.
4. For all pairs, try sort-and-prune on real data first and measure the number of x-overlapping pairs.
5. Use Bentley-Ottmann, or a library such as CGAL or GEOS, when n is large and k is small.
6. Log n, k and running time in production so you notice when the data moves into a regime the algorithm handles badly.