Point in Polygon — Vertex Intersection Bug

Point-in-polygon (PIP) is one of the most frequently implemented computational geometry algorithms. GIS systems (is this GPS coordinate in this city boundary?), game physics (is the player inside this room?), UI hit-testing (was this click inside this region?), and graphics rendering all need PIP. Two main algorithms: ray casting (even-odd rule, simple) and winding number (handles self-intersecting polygons, more robust).

Ray Casting Algorithm

The ray casting algorithm (even-odd rule) is the simplest and fastest method for simple polygons. Cast a horizontal ray from the point to +infinity in the x-direction. Count every edge that crosses this ray. If the count is odd, the point is inside; even means outside. The implementation uses a crossing test that checks if the edge's y-range spans the point's y-coordinate and if the intersection x of the edge with the ray is to the right of the point.

The crossing logic hides a subtlety: when the ray passes exactly through a vertex, the edge on one side may count the vertex as a crossing while the adjacent edge may not. The standard fix is to treat the vertex as belonging to the edge with the higher y-coordinate (strict upper vertex rule). This ensures the crossing count remains correct even for vertex-aligned rays.

def point_in_polygon_ray(point, polygon) -> bool:
    """
    Ray casting algorithm: cast ray from point to +x infinity.
    Count edge crossings — odd = inside, even = outside.
    Handles vertex intersections via strict upper vertex rule.
    """
    x, y = point
    n = len(polygon)
    inside = False
    j = n - 1
    for i in range(n):
        xi, yi = polygon[i]
        xj, yj = polygon[j]
        # Strict upper vertex: count crossing only if vertex with greater y is start
        if ((yi > y) != (yj > y)) and \
           (x < (xj - xi) * (y - yi) / (yj - yi) + xi):
            inside = not inside
        j = i
    return inside

# Square polygon
square = [(0,0),(4,0),(4,4),(0,4)]
print(point_in_polygon_ray((2,2), square))   # True (inside)
print(point_in_polygon_ray((5,2), square))   # False (outside)
print(point_in_polygon_ray((0,0), square))   # False (vertex — edge case)

Winding Number Algorithm

The winding number algorithm computes how many times the polygon winds around the point. It does this by summing signed angles or, more efficiently, by incrementing/decrementing a counter when edges cross a reference line. The result is an integer: zero means outside, non-zero means inside. The winding number is more robust than ray casting because it handles self-intersecting polygons (like star shapes) correctly.

Implementation: start with wn = 0. For each edge, check if the point's y-coordinate is between the edge's y's. If the edge crosses above the point (upward crossing), increment wn; if it crosses below (downward crossing), decrement. The orientation (left turn or right turn) determines the sign. At the end, wn = 0 means outside, |wn| > 0 means inside.

def winding_number(point, polygon) -> int:
    """Winding number: non-zero = inside."""
    x, y = point
    wn = 0
    n = len(polygon)
    for i in range(n):
        x1, y1 = polygon[i]
        x2, y2 = polygon[(i+1) % n]
        if y1 <= y:
            if y2 > y and orientation((x1,y1),(x2,y2),(x,y)) > 0:
                wn += 1
        else:
            if y2 <= y and orientation((x1,y1),(x2,y2),(x,y)) < 0:
                wn -= 1
    return wn

square = [(0,0),(4,0),(4,4),(0,4)]
print(winding_number((2,2), square))  # 1 (inside)
print(winding_number((5,2), square))  # 0 (outside)

# Self-intersecting star shape (bowtie)
bowtie = [(0,0),(2,2),(0,4),(2,2),(4,4),(2,2),(4,0),(2,2)]
print(winding_number((2,2), bowtie))  # 0? Actually the winding number handles it; inside classification is ambiguous

Edge Cases in Ray Casting — Vertex and Edge Alignments

The most common production failures in ray casting stem from numerical instability when the ray passes exactly through a vertex or lies collinear with an edge. A vertex on the ray gets counted twice if both adjacent edges are considered crossing. The standard solution: treat the vertex as belonging only to the edge with the higher y-coordinate (or lower, as long as consistent).

Similarly, when the ray overlaps an edge, the algorithm must not count that edge as a crossing. The crossing test naturally excludes horizontal edges because the condition '(yi > y) != (yj > y)' is false when yi == yj == y. But floating-point equality can fail. Use an epsilon tolerance: if abs(y - yi) < EPS and abs(y - yj) < EPS, skip the edge.

Another edge case: the point lies exactly on an edge. Both algorithms are undefined — is the point inside, outside, or on the boundary? Application-specific rules decide. In GIS, a point on a boundary is often considered inside.

Performance Optimisation for Multiple Queries

When you need to test thousands or millions of points against the same polygon, the O(n) per query becomes too slow. Preprocessing the polygon with a spatial index reduces average query complexity to O(log n) or O(1). Common approaches:

1. Grid subdivision: Divide the polygon's bounding box into a grid. For each cell, store which edges cross that cell. A point's query first checks which grid cell it falls in, then tests only edges in that cell.
2. Slab decomposition: For convex polygons, sort vertices by angle around a point (e.g., centroid) and binary search to find which wedge the point falls into, then test against that edge. O(log n) per query.
3. R-tree: Build an R-tree over polygon edges. Query the point's location — the tree returns candidate edges that intersect the point's horizontal ray. Works well for complex polygons with many edges.

For vectorised batch queries, use NumPy to compute crossings for all points simultaneously. This leverages SIMD operations and can process millions of points per second against a single polygon.

import numpy as np

def batch_point_in_polygon_ray(points: np.ndarray, polygon: np.ndarray) -> np.ndarray:
    """
    Vectorised ray casting for many points against same polygon.
    points: (N, 2) array of x,y coordinates.
    polygon: (M, 2) array of vertices (closed first-last).
    Returns boolean array of shape (N,).
    """
    n = len(polygon)
    inside = np.zeros(len(points), dtype=bool)
    j = n - 1
    for i in range(n):
        xi, yi = polygon[i]
        xj, yj = polygon[j]
        # Edge spans point y? (strict upper vertex included)
        cond1 = ((yi > points[:,1]) != (yj > points[:,1]))
        # Intersection x is to the right of point x
        cond2 = points[:,0] < (xj - xi) * (points[:,1] - yi) / (yj - yi + 1e-12) + xi
        inside ^= cond1 & cond2
        j = i
    return inside

# Example: test 1 million points against a square
square = np.array([(0,0),(4,0),(4,4),(0,4)])
points = np.random.rand(1_000_000, 2) * 5
result = batch_point_in_polygon_ray(points, square)
print(f"Inside count: {np.sum(result)} out of {len(points)}")

Choosing Between Ray Casting and Winding Number

Both algorithms run in O(n) time and O(1) memory for a single query. The choice depends on polygon properties and numerical requirements.

Ray casting pros: Simple, fast, works well for convex and simple concave polygons. Easy to vectorise. Fewer floating-point operations per edge. Ray casting cons: Fails on self-intersecting polygons. Vertex intersections require careful handling. Does not distinguish between winding number magnitude (useful for some applications).

Winding number pros: Handles any polygon type, including self-intersecting. Produces an integer that indicates how many times the polygon winds around the point (useful for regions with holes or nested boundaries). More consistent near boundaries. Winding number cons: Slightly more complex — requires orientation test. Harder to vectorise (though possible). Slightly more arithmetic per edge.

In practice, most production systems use winding number because real-world polygon data is messy. For performance-critical applications with clean convex polygons, ray casting with proper edge case handling is the better choice.