Line Segment Intersection (Bentley-Ottmann)

How the sweepline finds every crossing

The naive way to find intersections among n segments is to test all n(n−1)/2 pairs. That is O(n²): for 10,000 segments you run roughly 50 million intersection tests even if not a single pair actually crosses. Jon Bentley and Thomas Ottmann's 1979 algorithm gets the cost down to O((n + k) log n), where k is the number of intersections actually reported. When the arrangement is sparse — say a road network where most segments are far apart — this is the difference between milliseconds and minutes.

The trick is a plane sweep. Imagine a vertical line sweeping from left to right across the plane. At any instant it crosses some subset of the segments, and those segments have a well-defined top-to-bottom order along the sweepline. The key insight: two segments can only intersect if, at some moment just before the crossing, they are adjacent in that top-to-bottom order. So instead of testing all pairs, you only ever test segments that become neighbours — and there are far fewer of those.

Two data structures cooperate:

• The event queue — a priority queue of x-coordinates where something interesting happens, ordered by x (ties broken by y). It answers "what happens next?"
• The status structure — a balanced binary search tree holding the segments currently crossed by the sweepline, ordered by their y-value at the current x. It answers "who is next to whom right now?"

The sweep starts with all 2n endpoints in the queue and processes events in x-order. Each event is one of three kinds.

The three event types

Left endpoint. A new segment enters the sweep. Insert it into the status structure at its correct y-position. It now has an upper neighbour and a lower neighbour. Test the new segment against each of them; if either crossing lies to the right of the sweep, schedule it as a future event.

Right endpoint. A segment leaves the sweep. Remove it from the status structure. Its former upper and lower neighbours are now adjacent to each other — test that newly formed pair, because they may cross further right.

Intersection point. Two segments cross. Report the point, then swap their order in the status structure (the one that was on top is now on the bottom). After the swap each of the two has a new outer neighbour; test those two new pairs. This is what lets the sweep keep discovering crossings as it advances.

Every event does a constant number of O(log n) operations: one queue extraction, a few tree insert/delete/neighbour lookups, and at most a couple of geometric tests. There are 2n endpoint events plus k intersection events, so the total is O((n + k) log n).

When to reach for Bentley-Ottmann

• Map overlay and GIS. Overlaying two road layers, finding where rivers cross administrative boundaries, or clipping polygons all reduce to segment intersection. This was the original motivating application.
• Boolean operations on polygons. Union, intersection, and difference of polygons start by finding all edge crossings; the Greiner-Hormann and Vatti clippers depend on it.
• VLSI and CAD design-rule checks. Detecting overlapping wires or shapes on a chip layout is segment (and rectangle) intersection at scale.
• Self-intersection tests. Deciding whether a polygon is simple (non-self-crossing) is just "are there zero intersections among my edges?", answerable in O(n log n) by stopping at the first crossing.

If you only need to know whether any two segments cross — not enumerate all of them — you can stop the sweep at the first intersection event and get a clean O(n log n) with no dependence on k. Use brute force only when n is tiny (say under 20) or when you expect k to approach n², where the log-factor overhead of the sweep loses to a tight pairwise loop.

Bentley-Ottmann vs other intersection methods

	Brute force	Bentley-Ottmann	Chazelle-Edelsbrunner	Balaban (1995)	Uniform grid / bucketing
Report all k crossings	O(n²)	O((n + k) log n)	O(n log n + k)	O(n log n + k)	O(n + k) expected
Any crossing? (detect)	O(n²)	O(n log n)	O(n log n)	O(n log n)	O(n) expected
Extra space	O(1)	O(n + k)	O(n)	O(n)	O(n + grid)
Asymptotically optimal	only when k ≈ n²	no (log on k)	yes	yes	no (input-dependent)
Implementation difficulty	trivial	moderate	hard (hereditary segments)	hard	easy but tuning-sensitive
Degeneracy handling	none needed	fiddly (ties, verticals)	fiddly	fiddly	fiddly
Worst-case robustness	excellent	good	good	good	poor (adversarial clustering)

The asymptotically optimal algorithms (Chazelle-Edelsbrunner, Balaban) remove the log n factor on the output term, but they are markedly harder to implement and rarely worth it in practice. Bentley-Ottmann is the one in CGAL, in textbooks, and in nearly every production GIS pipeline because the constant factors are small and the code is comprehensible.

What the numbers actually say

• 10,000 segments, sparse (k ≈ 5,000): brute force runs ~50,000,000 pair tests; Bentley-Ottmann does ~25,000 events (2n endpoints + k crossings) × ~14 tree steps ≈ 350,000 operations — well over a 100× reduction in work.
• The log-factor penalty is real. When k is large — k ≈ n²/4 — the sweep does O(n² log n) work, strictly worse than the brute force's O(n²). Profile your expected k before assuming the sweep wins.
• Output can dwarf input. With n segments you may have up to n(n−1)/2 intersections — for n = 1,000 that is ~500,000 points. The O(n + k) space term, not the O(n) status structure, is usually what blows your memory budget.
• Events are cheap individually. Each event touches O(log n) tree nodes — about 14 for n = 10,000 — plus one or two cross-product orientation tests, each a handful of multiplies. The whole sweep is memory-bandwidth-bound long before it is CPU-bound.

JavaScript implementation

A compact, didactic version. It uses a sorted array as the status structure (O(n) updates, fine for teaching) and a binary-heap-style event queue; swap both for a balanced BST and a proper priority queue to reach the true bound.

// Orientation via cross product: >0 left turn, <0 right, 0 collinear.
const cross = (o, a, b) =>
  (a[0] - o[0]) * (b[1] - o[1]) - (a[1] - o[1]) * (b[0] - o[0]);

// Proper segment intersection point (returns null if they don't cross).
function intersect(s1, s2) {
  const [p, p2] = s1, [q, q2] = s2;
  const d1 = cross(q, q2, p), d2 = cross(q, q2, p2);
  const d3 = cross(p, p2, q), d4 = cross(p, p2, q2);
  if (((d1 > 0) !== (d2 > 0)) && ((d3 > 0) !== (d4 > 0))) {
    const r = [p2[0] - p[0], p2[1] - p[1]];
    const s = [q2[0] - q[0], q2[1] - q[1]];
    const denom = r[0] * s[1] - r[1] * s[0];
    const t = ((q[0] - p[0]) * s[1] - (q[1] - p[1]) * s[0]) / denom;
    return [p[0] + t * r[0], p[1] + t * r[1]];
  }
  return null;            // collinear / endpoint touches omitted for brevity
}

function bentleyOttmann(segments) {
  // Normalize each segment so endpoint[0] is the left (smaller x) end.
  const segs = segments.map(([a, b]) => (a[0] < b[0] || (a[0] === b[0] && a[1] < b[1])) ? [a, b] : [b, a]);

  const events = [];      // {x, y, type, i, j}  type: 0=left 1=right 2=cross
  for (let i = 0; i < segs.length; i++) {
    events.push({ x: segs[i][0][0], y: segs[i][0][1], type: 0, i });
    events.push({ x: segs[i][1][0], y: segs[i][1][1], type: 1, i });
  }
  const cmpEvent = (a, b) => a.x - b.x || a.y - b.y || a.type - b.type;
  events.sort(cmpEvent);

  const status = [];      // segment indices, ordered top-to-bottom at sweep x
  const result = [];
  const seen = new Set(); // dedupe intersection events

  // y of segment i at the sweep line x
  const yAt = (i, x) => {
    const [[x1, y1], [x2, y2]] = segs[i];
    return x1 === x2 ? y1 : y1 + (y2 - y1) * (x - x1) / (x2 - x1);
  };

  const schedule = (i, j, x) => {
    if (i == null || j == null) return;
    const pt = intersect(segs[i], segs[j]);
    if (pt && pt[0] >= x - 1e-9) {
      const key = pt[0].toFixed(6) + ',' + pt[1].toFixed(6) + ',' + Math.min(i, j) + ',' + Math.max(i, j);
      if (!seen.has(key)) {
        seen.add(key);
        events.push({ x: pt[0], y: pt[1], type: 2, i, j });
        events.sort(cmpEvent);  // re-sort: replace with a heap insert in production
      }
    }
  };

  while (events.length) {
    const ev = events.shift();
    const x = ev.x;
    const posOf = (i) => status.indexOf(i);

    if (ev.type === 0) {                 // left endpoint
      let lo = 0;
      while (lo < status.length && yAt(status[lo], x) > yAt(ev.i, x)) lo++;
      status.splice(lo, 0, ev.i);
      schedule(status[lo - 1], ev.i, x);
      schedule(ev.i, status[lo + 1], x);
    } else if (ev.type === 1) {          // right endpoint
      const p = posOf(ev.i);
      const above = status[p - 1], below = status[p + 1];
      status.splice(p, 1);
      schedule(above, below, x);
    } else {                             // intersection
      result.push([ev.x, ev.y]);
      const pi = posOf(ev.i), pj = posOf(ev.j);
      const lo = Math.min(pi, pj), hi = Math.max(pi, pj);
      [status[lo], status[hi]] = [status[hi], status[lo]];   // swap order
      schedule(status[lo - 1], status[lo], x);
      schedule(status[hi], status[hi + 1], x);
    }
  }
  return result;
}

Two things to flag. First, the seen set deduplicates intersection events — without it the same crossing gets scheduled from both neighbour tests and reported twice. Second, the events.sort() after each insert is the part that betrays the bound: a real implementation pushes into a heap or balanced tree in O(log n), not re-sorts the whole array.

Python implementation

The same algorithm using heapq for the event queue, which is the natural way to keep the O(log n) insert. The status structure is still a list here for clarity.

import heapq

def cross(o, a, b):
    return (a[0]-o[0])*(b[1]-o[1]) - (a[1]-o[1])*(b[0]-o[0])

def intersect(s1, s2):
    (p, p2), (q, q2) = s1, s2
    d1, d2 = cross(q, q2, p), cross(q, q2, p2)
    d3, d4 = cross(p, p2, q), cross(p, p2, q2)
    if (d1 > 0) != (d2 > 0) and (d3 > 0) != (d4 > 0):
        rx, ry = p2[0]-p[0], p2[1]-p[1]
        sx, sy = q2[0]-q[0], q2[1]-q[1]
        denom = rx*sy - ry*sx
        t = ((q[0]-p[0])*sy - (q[1]-p[1])*sx) / denom
        return (p[0]+t*rx, p[1]+t*ry)
    return None

def y_at(seg, x):
    (x1, y1), (x2, y2) = seg
    if x1 == x2:
        return y1
    return y1 + (y2-y1)*(x-x1)/(x2-x1)

def bentley_ottmann(segments):
    segs = [s if s[0][0] < s[1][0] or (s[0][0]==s[1][0] and s[0][1]<s[1][1])
            else (s[1], s[0]) for s in segments]

    pq = []                       # (x, y, type, i, j)  type 0=left 1=right 2=cross
    for i, s in enumerate(segs):
        heapq.heappush(pq, (s[0][0], s[0][1], 0, i, -1))
        heapq.heappush(pq, (s[1][0], s[1][1], 1, i, -1))

    status, result, seen = [], [], set()

    def schedule(i, j, x):
        if i is None or j is None:
            return
        pt = intersect(segs[i], segs[j])
        if pt and pt[0] >= x - 1e-9:
            key = (round(pt[0], 6), round(pt[1], 6), min(i, j), max(i, j))
            if key not in seen:
                seen.add(key)
                heapq.heappush(pq, (pt[0], pt[1], 2, min(i, j), max(i, j)))

    while pq:
        x, y, typ, i, j = heapq.heappop(pq)
        if typ == 0:                      # left endpoint
            lo = 0
            while lo < len(status) and y_at(segs[status[lo]], x) > y_at(segs[i], x):
                lo += 1
            status.insert(lo, i)
            schedule(status[lo-1] if lo > 0 else None, i, x)
            schedule(i, status[lo+1] if lo+1 < len(status) else None, x)
        elif typ == 1:                    # right endpoint
            p = status.index(i)
            above = status[p-1] if p > 0 else None
            below = status[p+1] if p+1 < len(status) else None
            status.pop(p)
            schedule(above, below, x)
        else:                             # intersection
            result.append((x, y))
            pi, pj = status.index(i), status.index(j)
            lo, hi = min(pi, pj), max(pi, pj)
            status[lo], status[hi] = status[hi], status[lo]
            schedule(status[lo-1] if lo > 0 else None, status[lo], x)
            schedule(status[hi], status[hi+1] if hi+1 < len(status) else None, x)
    return result

Note the asymmetry between the three handlers. Left endpoints test two new pairs (the segment against both neighbours); right endpoints test exactly one (the pair that closes up behind the departing segment); intersections test two (the new outer neighbours after the swap). Getting these neighbour tests exactly right — and only these — is the whole correctness argument.

Variants worth knowing

Detection-only sweep. Shamos and Hoey's 1976 algorithm — the ancestor of Bentley-Ottmann — answers only "do any two segments cross?" in O(n log n) with no k term. Stop the sweep at the first intersection event. This is the right tool for simple-polygon testing.

Chazelle-Edelsbrunner (1992). The first O(n log n + k) reporting algorithm, removing the log factor from the output term. It maintains a more elaborate structure of "hereditary segments" and is the asymptotically optimal answer, but the constant factors and code complexity keep it out of most libraries.

Balaban's algorithm (1995). Also O(n log n + k), using a divide-and-conquer "staircase" decomposition instead of a single sweep. Conceptually cleaner than Chazelle-Edelsbrunner and easier to make robust, though still rarely implemented.

Red-blue (bichromatic) intersection. When segments come in two non-self-intersecting sets (red and blue) and you only want red-blue crossings, specialized algorithms run in O(n log n + k) and skip same-color tests entirely. This is exactly the map-overlay case: two clean layers crossing each other.

Randomized incremental construction. Instead of sweeping, insert segments one at a time into a trapezoidal decomposition. Expected O(n log n + k) and the basis of the robust intersection routines in CGAL's arrangement package.

Common bugs and edge cases

• Floating-point in the orientation test. The cross-product sign decides everything. Near-collinear inputs flip the sign and the sweep mis-orders the status structure, missing or duplicating crossings. Use integer coordinates or exact/adaptive predicates (Shewchuk) for anything production.
• Vertical segments. A vertical segment has a single x for both endpoints, breaking the x-ordered event sweep. Standard fixes: rotate the plane by a tiny angle, or define a consistent tie-break (process vertical segments specially at their x).
• Double-reporting intersections. A crossing is reachable from both neighbour tests on either side. Without a "seen" set keyed on the point plus the segment pair, you report it twice — or worse, swap the segments twice and corrupt the order.
• Three or more segments through one point. The clean "swap two segments" model breaks when several cross at the same point. You must gather all segments passing through the event point and reverse their whole sub-order at once.
• Scheduling intersections to the left of the sweep. Only schedule a crossing whose x is at or right of the current sweep x. A crossing already behind the line was either handled or is an artifact of an order error — re-adding it loops forever.
• Assuming the sweep always wins. When k approaches n², the O((n + k) log n) sweep does more work than an O(n²) brute force. For dense arrangements, measure before you commit.