Cairn
algorithms · computational geometry · O(log n) expected query, O(n log n) expected construction

back to Geometry

Trapezoidal Map

Both Point in Polygon and Slab Decomposition name this page in their own Complexity sections and leave it unbuilt. Slab decomposition already gets point location down to two binary searches, but it pays for that with a preprocessing structure that can be genuinely quadratic — a slab's edge list is rebuilt from every edge that spans it, and an adversarial polygon can put nearly every edge in nearly every slab (see that page's own measured 32-tooth/64-tooth blowup). A trapezoidal map decomposes the plane the same conceptual way — extend a vertical segment up and down from every polygon vertex until it hits another edge or the bounding box, carving the plane into trapezoids each bounded by at most one edge above, at most one edge below, and two vertical sides — but it never rebuilds a slab's whole edge list from scratch. It's built by inserting the polygon's edges one at a time, in random order, into a bounding box that starts as a single trapezoid. Each new edge splits only the trapezoids it actually crosses, and a small search structure — a DAG of "which side of this vertex" and "which side of this edge" comparisons — gets grafted in wherever a split happens. A point-location query just walks that DAG from the root: no binary search over slabs, no per-slab edge list, one O(log n) expected walk straight to the trapezoid containing the query point.

One deliberate simplification from the textbook algorithm, said honestly up front: the standard construction (de Berg, Cheong, van Kreveld, and Overmars, Computational Geometry: Algorithms and Applications, ch. 6) also maintains direct trapezoid-to-trapezoid neighbor pointers, so that finding which trapezoids a new edge crosses takes O(1) work per trapezoid instead of a fresh search each time. This build drops those neighbor pointers entirely and finds crossed trapezoids the same way it answers a query — by walking the search DAG from the root — trading an extra O(log n) factor per insertion step (invisible at this page's scale, a handful of edges) for having exactly one piece of logic that needs to be correct instead of two. That trade paid off directly during development: an early version with neighbor pointers had a real, silent bug (an unmerged trapezoid several insertions later would end up with the wrong extent because a neighbor pointer wasn't retargeted correctly somewhere upstream); dropping them for the DAG-only approach below removed that whole bug class by removing the code that could have it. See Pitfalls for the bug that turned up in the simpler version instead, and the numbers from stress-testing it.

Try it

A seven-vertex arrow — the same silhouette as Point in Polygon and Slab Decomposition's own arrow, with every vertex nudged a few pixels off whatever x-coordinate it used to share with another vertex. That's not cosmetic: a trapezoidal map's vertical extensions need a well-defined single edge immediately above and below every x-coordinate, and the original arrow has three genuinely vertical edges and several repeated x's, which is exactly the degenerate case this build doesn't handle (see Pitfalls). The dashed cell outlines below are the actual trapezoidal decomposition — 19 trapezoids from 7 edges, built by inserting the edges in one fixed random order so the page renders the same decomposition every time. Pick a point below and watch the search DAG walk down to the trapezoid that contains it.

trapezoid:
Pick a point above.

Why it works

The DAG's depth on any one root-to-leaf path is exactly the number of times the trapezoid containing that path's point got split over the course of construction, plus one. Fix a query point q that isn't on any edge. Before any edges are inserted, q sits in the single bounding-box trapezoid — depth 1. Each time an inserted edge happens to cross the trapezoid q currently lives in, that trapezoid gets split and q's depth increases by one (or two, if q's trapezoid also needed a left or right remainder). The question "how many of the n edges, inserted in a uniformly random order, split q's trapezoid?" is answered by the standard backward analysis for randomized incremental algorithms: look at the final structure and remove the last-inserted edge. q's final trapezoid merges back into a larger one exactly when the removed edge was one of the (at most four) edges bounding it. Since the insertion order is uniformly random, any one of the k edges present at step k is equally likely to have been inserted last, so the probability the k-th edge caused a split at q is at most 4/k. Summed over all n edges, the expected number of splits along q's path is at most 4·(1/1 + 1/2 + ... + 1/n) = O(log n) — the harmonic series, the same sum that gives quickselect and randomized quicksort their expected O(n log n) and O(n) bounds. That bound holds for every point q, not just the ones this demo happens to query, which is what makes "expected O(log n) query" a property of the structure itself rather than a property of any particular test point.

Reference implementation

function insertSegment(dag, seg) {
  // seg.p1 is lexicographically left of seg.p2 (x, then y, to order
  // vertical-adjacent points consistently without a real shear transform)
  const crossed = findCrossedTrapezoids(dag, seg); // walks the DAG, see below
  crossed.forEach((old, idx) => {
    const first = idx === 0, last = idx === crossed.length - 1;
    const sliceLeft = first ? seg.p1 : old.leftp;
    const sliceRight = last ? seg.p2 : old.rightp;
    const top = makeTrap(sliceLeft, sliceRight, old.top, seg);
    const bottom = makeTrap(sliceLeft, sliceRight, seg, old.bottom);
    let subtree = makeYNode(seg, leaf(top), leaf(bottom));
    if (first && !pointsEqual(old.leftp, seg.p1)) {
      const leftRemainder = makeTrap(old.leftp, seg.p1, old.top, old.bottom);
      subtree = makeXNode(seg.p1, leaf(leftRemainder), subtree);
    }
    if (last && !pointsEqual(old.rightp, seg.p2)) {
      const rightRemainder = makeTrap(seg.p2, old.rightp, old.top, old.bottom);
      subtree = makeXNode(seg.p2, subtree, leaf(rightRemainder));
    }
    spliceIntoDag(old.dagNode, subtree); // in place, so nothing else needs updating
  });
}

function locate(dag, pt, forSeg) {
  let node = dag.root;
  while (node.type !== 'leaf') {
    if (node.type === 'x') {
      node = pt.x >= node.point.x ? node.right : node.left;
    } else {
      const seg = node.seg;
      const onEndpoint = pointsEqual(pt, seg.p1) || pointsEqual(pt, seg.p2);
      const isAbove = onEndpoint
        ? orient(seg.p1, seg.p2, forSeg.p2) > 0   // shared vertex: use the new
        : orient(seg.p1, seg.p2, pt) > 0;          // edge's own direction to break the tie
      node = isAbove ? node.above : node.below;
    }
  }
  return node.trap;
}

function findCrossedTrapezoids(dag, seg) {
  const path = [locate(dag, seg.p1, seg)];
  let trap = path[0];
  while (lexLess(trap.rightp, seg.p2)) {
    const x = trap.rightp.x;
    const y = yOnSegment(seg, x); // seg's OWN y here, not trap.rightp's y
    trap = locate(dag, { x, y }, seg);
    path.push(trap);
  }
  return path;
}

Pitfalls

A vertex's own x-coordinate has to be compared as a plain number, not as a (x, y) tie-broken pair — using the wrong one turns the crossed-trapezoid walk into an infinite loop. findCrossedTrapezoids above deliberately re-queries the DAG at x = trap.rightp.x — exactly the vertex that ends the current trapezoid — with a y that comes from the segment being inserted, which generally has nothing to do with that vertex's own y. An earlier version of locate's x-node comparison used the same lexicographic (x, then y) order everywhere in the codebase, on the reasoning that it was already the order used to pick each segment's left endpoint. At an x-node exactly matching trap.rightp, that sent the query left back into the trapezoid it just came from whenever the vertex's own y happened to exceed the query's y — undoing the rightward walk and looping forever, since nothing about the situation ever changes on the next iteration. Caught immediately: the very first stress-test run (2,000+ random simple polygons, cross-checked against an independent brute-force oracle — see below) crashed the harness's own runaway-loop guard on 6 of the first 7 polygons tried, before the harness's crash cap stopped it early. The fix is the one used in the reference code above: x-nodes compare pt.x >= node.point.x alone, full stop, no y tie-break. After the fix: 0 crashes and 0 mismatches across five separate random seeds (roughly 14,500 polygons, 436,000 point-location queries in total), plus a separate run up to 34-vertex polygons (2,000 trials, 60,000 queries, still 0 mismatches).

General position is assumed, the same simplifying assumption every textbook presentation of this algorithm makes: no polygon edge is perfectly vertical, and no two vertices share an x-coordinate that also happens to matter (see "Try it" above for why this demo's own arrow is nudged off the original arrow's shared x's). The lexicographic (x, then y) point order sidesteps the "which vertex is really first" question without a real shear transform, but a genuinely vertical edge has no well-defined single y per x at all, which this implementation doesn't attempt to handle. Production-quality implementations resolve this with symbolic perturbation — treat every coordinate as infinitesimally rotated so no two x's ever tie exactly — rather than skipping the case; that's real added machinery this build doesn't include.

A point-location query landing exactly on a polygon vertex or exactly on an edge has no single mathematically correct trapezoid — it's genuinely on the boundary between two or more. This build's locate resolves that case by always taking the "above" branch when there's no segment direction available to disambiguate (the exactly on a vertex preset above demonstrates it directly). That's a consistent, defensible convention — the same kind of "pick one side and document it" call Slab Decomposition's own Pitfalls makes for its slab-boundary rule — but it's still an arbitrary choice, not a derived one, and a caller that cares about exact-boundary semantics needs to know it's there.

Correctness was checked against an independently-built oracle, not just against this page's own demo polygon. The oracle brute-forces the same decomposition a completely different way — slice the plane into vertical slabs at every vertex x (the same idea Slab Decomposition uses), classify each slab's regions inside/outside by parity, then merge horizontally-adjacent slab regions that share the same bounding edges into full trapezoids — and every region's inside/outside tag was cross-checked a third way, against a plain crossing-number point-in-polygon test. All three agreeing across the runs cited above is a stronger signal than any one of them agreeing with itself.

Complexity

Query: O(log n) expected, walking the search DAG from root to leaf — see Why it works for the backward-analysis argument, and the demo's own measured numbers (19 trapezoids from 7 edges, average 3.7 steps and a maximum of 8 across 5,000 random query points, against log2 19 ≈ 4.2). Construction: O(n log n) expected time and O(n) expected space with the textbook's neighbor-pointer approach; this build's DAG-only simplification (see above) costs an extra expected O(log n) factor per insertion step, so O(n log² n) expected time overall — still a genuine improvement over Slab Decomposition's O(n²) worst case, just not quite the textbook bound. Space stays O(n) expected either way: each of the n insertions creates at most 4 new trapezoids, regardless of how those trapezoids get found.

This site's guide, Choosing a Geometry Algorithm, places this entry as the worst-case-resistant alternative to Slab Decomposition for the same repeated point-in-polygon query, trading a simpler build for a randomization-backed guarantee.