Kirkpatrick Point Location Algorithm

Planar point location asks: given a subdivision of the plane into polygonal faces, preprocess it so that for any query point we can quickly report which face contains it. Kirkpatrick's algorithm, published in 1983, was the first to achieve the optimal bounds of $O(n)$ preprocessing, $O(n)$ space, and $O(\log n)$ query time for a triangulated subdivision of $n$ vertices.

The core idea is to build a hierarchy of progressively coarser triangulations. At each level, an independent set of low-degree vertices is removed and the resulting holes are re-triangulated. Each new triangle "points to" the set of triangles it overlaps in the previous (finer) level. A query walks from the coarsest level (one big triangle) down through at most $O(\log n)$ levels, at each step moving from the current triangle to the unique child triangle that contains the query point. The clever invariant — that each coarser triangle intersects only $O(1)$ finer triangles, thanks to the bounded-degree independent set — is what makes the descent $O(\log n)$ rather than $O(\log^2 n)$.

The implementation below is simplified. Building the full hierarchy in a short snippet is impractical, so instead we expose the same interface — "given a planar subdivision and a query point, return the face index that contains it" — and implement it by linear search over the faces with a ray-casting point-in-polygon test. The behaviour is the same as Kirkpatrick point location would provide, only the query time is $O(n)$ instead of $O(\log n)$.

Codes

Python

def solve(xs):
    """Simplified Kirkpatrick point location.

    Input encoding: xs = [query_point, polygon_1, polygon_2, ...]
    where query_point is [qx, qy] and each polygon is a list of
    [x, y] vertices in order. Returns [index] of the first polygon
    that contains the query point, or [-1] if none do.
    """
    if not xs:
        return [-1]
    query = xs[0]
    polygons = xs[1:]

    def point_in_polygon(pt, poly):
        x, y = pt[0], pt[1]
        inside = False
        n = len(poly)
        j = n - 1
        for i in range(n):
            xi, yi = poly[i][0], poly[i][1]
            xj, yj = poly[j][0], poly[j][1]
            if (yi > y) != (yj > y):
                x_cross = (xj - xi) * (y - yi) / (yj - yi) + xi
                if x < x_cross:
                    inside = not inside
            j = i
        return inside

    for idx, poly in enumerate(polygons):
        if point_in_polygon(query, poly):
            return [idx]
    return [-1]

C++

#include <vector>
#include <array>

using Point = std::array<double, 2>;

bool point_in_polygon(const Point& pt, const std::vector<Point>& poly) {
    bool inside = false;
    int n = static_cast<int>(poly.size());
    int j = n - 1;
    for (int i = 0; i < n; ++i) {
        double xi = poly[i][0], yi = poly[i][1];
        double xj = poly[j][0], yj = poly[j][1];
        if ((yi > pt[1]) != (yj > pt[1])) {
            double x_cross = (xj - xi) * (pt[1] - yi) / (yj - yi) + xi;
            if (pt[0] < x_cross) inside = !inside;
        }
        j = i;
    }
    return inside;
}

int locate(const Point& query, const std::vector<std::vector<Point>>& polys) {
    for (int i = 0; i < static_cast<int>(polys.size()); ++i) {
        if (point_in_polygon(query, polys[i])) return i;
    }
    return -1;
}

Example

Loading Python runner...

Description

Run time analysis

The full Kirkpatrick structure answers queries in $O(\log n)$ after $O(n)$ preprocessing and uses $O(n)$ space. The bound on query time comes from descending through $O(\log n)$ levels of the hierarchy; each level shrinks by a constant factor because the removed independent set has linear size, and each descent step is $O(1)$ because the bounded-degree property ensures a coarser triangle overlaps only a constant number of finer ones. The simplified snippet here runs in $O(nk)$ time per query, where $k$ is the average polygon size, since it scans every face linearly.

Space analysis

$O(n)$ for the full hierarchy, since the sum of all triangles across levels is a geometric series dominated by the finest level. The simplified snippet uses $O(n)$ to store the input polygons and $O(1)$ additional memory for the query.

Proof of correctness

Kirkpatrick's hierarchy is correct by induction on the level. At the finest level, the triangulation exactly equals the input subdivision, so "the triangle containing the query" is the correct face. Each coarser level is obtained by removing a set of bounded-degree vertices and re-triangulating the resulting hole, which is itself a star-shaped (in fact, kernel-containing) polygon; every new triangle lies entirely within the union of the removed triangles. Therefore, descending from a coarse triangle into the unique finer triangle that contains the query point preserves the invariant "the current triangle contains the query." After $O(\log n)$ levels the descent terminates at a leaf triangle of the original subdivision, which is the answer. The simplified implementation above is correct trivially: ray casting classifies a point as inside a simple polygon iff a horizontal ray from the point crosses the polygon's boundary an odd number of times, a standard result due to Jordan and formalised in every computational-geometry textbook.

Extensions

Applications

References

• Kirkpatrick, D. G. (1983). "Optimal search in planar subdivisions." SIAM Journal on Computing, 12(1), 28-35.
• de Berg, M., Cheong, O., van Kreveld, M., & Overmars, M. Computational Geometry: Algorithms and Applications, 3rd ed. Springer, 2008, Chapter 6.
• Preparata, F. P., & Shamos, M. I. Computational Geometry: An Introduction. Springer, 1985.
• Snoeyink, J. (2004). "Point location." In Handbook of Discrete and Computational Geometry, 2nd ed., CRC Press.