Trapezoidal Maps and Randomized Point Location

Point location is the query “which face of this planar subdivision contains this point”. The subdivision is fixed, whether a map with countries, a triangulated terrain, or a set of polygon boundaries, and the queries arrive many times. A linear scan over all edges per query is correct and, for a subdivision of a few thousand edges queried millions of times, far too slow. The target is O(log n) per query with a structure of O(n) size, and there are several ways to reach it. The trapezoidal map is the one that generalises most directly to arbitrary subdivisions, including disconnected and non-convex ones, and it is the one whose construction is easiest to argue about.

The idea is to partition the plane into vertical slabs by drawing a vertical line through every vertex. Within a slab no vertex is present, so the edges crossing the slab are ordered top to bottom and never cross. Adding a vertical extension from each vertex to the edge above and below splits each slab into trapezoids: cells bounded above and below by edges and on the left and right by vertical segments. Each trapezoid belongs to exactly one face of the subdivision, so locating a point reduces to finding its trapezoid.

The structure is a DAG, not a tree

The search structure is not the map. It is a directed acyclic graph whose internal nodes are of two kinds. An x-node holds the x-coordinate of a vertex and sends the query left or right depending on which side it falls. A y-node holds an edge and sends the query above or below depending on which side of that edge’s supporting line it falls. Leaves are trapezoids.

A query walks from the root, taking one branch per node, and arrives at a leaf in time proportional to the path length. The path length is the depth of the leaf in the DAG, and the construction aims to make every root-to-leaf path short.

The two node types correspond exactly to the ways a new segment cuts the trapezoids it crosses. When a segment is inserted, it passes through a sequence of trapezoids. In a trapezoid it crosses from one vertical wall to the other without terminating, the split is a y-node into an upper and a lower piece, giving two leaves. Where an endpoint of the segment lies inside a trapezoid, the split additionally creates a piece to the left or right of that endpoint, so the leaf is replaced by a subtree with three leaves; when both endpoints lie inside the same trapezoid, the result has four leaves. Those are the only cases, and each replaces one leaf with a small subtree of two, three, or four leaves, so the total number of nodes stays linear in the number of trapezoids created.

Randomized incremental construction

Insert the segments one at a time in uniformly random order. Each insertion touches only the trapezoids the new segment crosses, which on average is a constant number. The reason is the backward analysis: a trapezoid that exists after all n segments are inserted was created by the last of the segments bounding it, and in a random order each of those segments is equally likely to be last, so the expected number of trapezoids created at step i is O(1). The total number of trapezoids ever created is O(n), and the total expected construction time is O(n log n) — the log coming from locating the endpoints and the crossed trapezoids in the current structure.

The same backward argument bounds the expected query time. A point’s search path in the final structure has length equal to the number of insertions that split the trapezoid containing it. At step i, the probability that the i-th inserted segment splits that particular trapezoid is bounded by 4/i, because a trapezoid is delimited by at most four segments and, in a random order, each of them is equally likely to be the last one inserted; the expected number of splits affecting any one point is therefore the harmonic sum up to n, which is O(log n). The bound holds for every query point, not just on average over query points, which is what makes the structure usable as a worst-case guarantee in expectation.

def locate(node, p):
    while node.kind != 'leaf':
        if node.kind == 'x':
            node = node.left if p.x < node.value else node.right
        else:
            node = node.above if side(node.edge, p) > 0 else node.below
    return node.trapezoid

The side test is the orientation predicate again, and it must be exact for the same reason as everywhere else: the leaf a query lands in determines the reported face, and a sign error near an edge reports the wrong face. The structure does not tolerate approximate comparisons, because the faces are defined by the edges and the query’s answer is a discrete label.

Why a DAG rather than a search tree

An ordinary binary search tree over the trapezoids would need the query’s position relative to a global ordering to descend, and there is no global ordering of trapezoids that a point can be compared against. The DAG encodes the history of splits instead: each node’s test is meaningful only in the region that reaches it, which is exactly the region its subtree covers. Different branches can apply different tests to the same coordinate, and that is the point.

The cost is that the DAG has more nodes than the final map has trapezoids, by a constant factor of a few. A leaf can acquire several incoming edges because inserting a segment merges trapezoids that came from different predecessors into one, and that merging is what keeps the node count linear rather than letting each path carry its own copy of the history.

Degeneracies and the bounding box

Two vertices with the same x-coordinate produce two x-nodes with the same value, and the search then depends on the tie-break. The standard fix is to sort vertices by x and break ties by y, and to treat the vertical line through a vertex as belonging to that vertex alone, so no two vertical walls coincide. Without that convention the slab decomposition has zero-width slabs, and the trapezoids in them are degenerate.

The plane is unbounded, so the construction encloses everything in a large bounding box whose edges count as subdivision edges. The box must be large enough to contain all vertices and all query points, which is a requirement the structure cannot enforce after construction; queries outside the box fall through to the outer trapezoids and report the unbounded face, which is correct only if the caller understands the box as part of the subdivision.

Vertical edges are the other classic difficulty. A vertical edge has no unique y-order against the slab’s other edges, so it cannot be treated as an ordinary segment. The usual solution is the same symbolic shear used in Bentley–Ottmann: replace each x-coordinate by x + εy for an infinitesimal ε, which makes every edge non-vertical and preserves all incidences.

Where this fits against the alternatives

Kirkpatrick’s hierarchy gives O(log n) worst-case query time with a triangulated subdivision and a hierarchy of progressively coarser triangulations, at the cost of a more elaborate construction and a structure that only applies to triangulations. The trapezoidal map gives O(log n) expected time for any subdivision and is simpler to implement, which is why it is the usual choice when the subdivision is not already triangulated.

For a triangulation, the Delaunay triangulation supports point location by walking from a known triangle toward the query, which visits O(√n) triangles in expectation for uniformly distributed points and covers only a short distance when successive queries are close together. The DAG wins when queries are adversarial or uniformly distributed over a large map; the walk wins when successive queries are correlated, which is the common case in tracking and simulation.

The k-d tree solves a superficially similar problem and is not a substitute. It indexes points, not faces, and a query for “which polygon contains this point” would need the polygon edges inserted as objects, at which point the tree answers “which edges are near this point” and the caller still has to do the containment test. The trapezoidal map answers the face question directly, in one descent, with no post-processing.

Memory is the practical constraint. The structure stores one node per split, which for n segments is a few times n nodes, each holding a coordinate, two child pointers, and a node kind. For a map with 10^6 edges that is on the order of 10^8 bytes, and it is built once and queried indefinitely, which is the right trade for any map that is read more than it is edited. Subdivisions that change frequently are better served by rebuilding the structure after batches of edits than by attempting incremental updates, since the randomized analysis assumes the insertion order and incremental updates destroy that assumption.