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).

How Point-in-Polygon Tests Actually Fail at Vertices

Point-in-polygon (PIP) determines whether a given point lies inside, outside, or exactly on the boundary of a polygon. The core mechanic is the ray casting algorithm: cast a ray from the point in any direction and count intersections with polygon edges. An odd count means inside; even means outside. This works for simple polygons in O(n) time per point, where n is the number of edges.

In practice, the algorithm's correctness hinges on how you handle edge cases — especially when the ray passes directly through a vertex. Standard implementations treat a ray-vertex intersection as a single edge crossing, but this miscounts when the vertex is shared by two edges on opposite sides of the ray. The result: points near vertices can be misclassified as inside or outside, producing silent logic errors in spatial queries.

Use PIP when you need to test containment in geographic boundaries, game hitboxes, or CAD regions. It matters because production systems — like geofencing for delivery zones or collision detection in physics engines — rely on exact boundary semantics. A vertex bug can cause a delivery truck to be incorrectly charged or a player to clip through a wall.

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.

Why You Shouldn't Roll Your Own — The Hidden Cost of Point-in-Polygon

Every junior I've onboarded wants to write a point-in-polygon check from scratch. It's a rite of passage. But in production, that code rarely stays isolated. It gets called inside a hot loop, exposed via an API, or used to validate geometry for a rendering pipeline.

Here's the kicker: most hand-rolled implementations fail silently on boundary cases. The ray-casting algorithm works beautifully for simple convex polygons. Throw in a concave polygon with a hole, and suddenly your 'inside' results flip at random. I've seen this bring down a GPS geofencing service in production — trucks reported as 'outside' when they were clearly inside the depot.

The real problem isn't the algorithm. It's the assumption that polygons are well-behaved. Real-world polygons have self-intersections, degenerate edges, and floating-point drift. That's why libraries like JTS (Java Topology Suite) or GEOS exist. They've spent years handling edge cases you haven't thought of yet.

Before you import a library, know this: the best point-in-polygon code is the code you don't write. If you must write it, spend equal time testing on malformed polygons — your future self (and your on-call rotation) will thank you.

// io.thecodeforge — dsa tutorial
// Scenario: A real-world GPS geofence check that 
// fails on concave polygons with self-intersections

import java.awt.geom.Path2D;

public class GeoFenceBypass {
    // This polygon has a self-intersecting edge
    static Path2D.Double buildBadFence() {
        Path2D.Double fence = new Path2D.Double();
        fence.moveTo(0, 0);
        fence.lineTo(10, 0);
        fence.lineTo(5, 10);  // self-intersects later
        fence.lineTo(10, 10); 
        fence.lineTo(0, 10);
        fence.closePath();
        return fence;
    }

    public static void main(String[] args) {
        Path2D.Double fence = buildBadFence();
        double pointX = 5, pointY = 5;
        // JTS would handle this; naive ray casting doesn't
        System.out.println("Inside? " + fence.contains(pointX, pointY));
    }
}

Graham Scan: The Sorting Sauce Behind Convex Hulls

You're debugging a geofence that uses a convex polygon as its boundary. The polygon was generated by a clustering algorithm, but vertices came out in random order. That's when you need a convex hull — the minimal convex polygon that encloses all points.

Graham Scan is the workhorse here. It's O(n log n) — dominated by sorting — and it's deceptively simple. You pick the lowest-y point (pivot), sort all other points by polar angle relative to that pivot, then sweep through, discarding any point that creates a clockwise turn (i.e., a right turn). The stack of surviving points is your convex hull.

The pitfall? The sorting step. If two points have the same angle relative to the pivot, you must keep the farthest one and discard the nearer one. If you don't, you'll get colinear vertices on the hull that break downstream algorithms (e.g., edge-crossing checks). I've seen this cause a 30-minute outage in a logistics system when a driver was incorrectly routed because the hull contained a degenerate 'spike'.

Practical advice: use the cross-product to determine turn direction. Cross-product > 0 means counter-clockwise (keep), < 0 means clockwise (discard). If cross-product is zero, points are colinear — keep the farthest from pivot.

// io.thecodeforge — dsa tutorial
// Graham Scan: Build convex hull from random points

import java.util.*;

public class ConvexHullDebugger {
    static class Pt {
        int x, y;
        Pt(int x, int y) { this.x = x; this.y = y; }
    }

    // Fix: handle colinear points by keeping farthest from pivot
    static Pt[] grahamScan(List<Pt> points) {
        // finds pivot (lowest y, then lowest x)
        Pt pivot = Collections.min(points, 
            (a, b) -> a.y != b.y ? a.y - b.y : a.x - b.x);
        // sort by polar angle — tie-break by distance
        points.sort((a, b) -> {
            int cross = cross(pivot, a, b);
            if (cross != 0) return -cross; // ccw first
            return dist(pivot, a) - dist(pivot, b);
        });
        // Build hull with stack
        Stack<Pt> hull = new Stack<>();
        for (Pt p : points) {
            while (hull.size() >= 2 && 
                   cross(hull.get(hull.size()-2), 
                         hull.peek(), p) <= 0) {
                hull.pop();
            }
            hull.push(p);
        }
        return hull.toArray(new Pt[0]);
    }
}

Time & Space: Why Your Polygon Query Is Slow

Before you write a single line of Point-in-Polygon (PiP) code, you need to know what you're paying for. Ray casting runs in O(n) per query — n being vertex count. Winding number is also O(n). No free lunch. If you're querying a million points against a 1000-vertex polygon, you're looking at a billion edge-crossing checks. That's not a bug — that's the math.

Space complexity is trivial: O(1) if you're given the polygon as an array; O(n) if you store it. The real cost is in time. For convex polygons, you can drop to O(log n) with binary search — that's the optimization worth chasing. For arbitrary polygons? You're stuck with linear.

The lesson: profile before you optimize. If your polygon has more than a few hundred vertices and you're running batch queries, you need spatial indexing (grids, R-trees) or a bounding volume hierarchy. Don't blame the algorithm for your data structure.

// io.thecodeforge — dsa tutorial

public class PointInPolygonComplexity {
    // Ray casting: O(n) time, O(1) space
    static boolean rayCast(int[][] polygon, int[] point) {
        boolean inside = false;
        int n = polygon.length;
        for (int i = 0, j = n - 1; i < n; j = i++) {
            int xi = polygon[i][0], yi = polygon[i][1];
            int xj = polygon[j][0], yj = polygon[j][1];
            if ((yi > point[1]) != (yj > point[1]) &&
                point[0] < (xj - xi) * (point[1] - yi) / (yj - yi) + xi) {
                inside = !inside;
            }
        }
        return inside;
    }
    
    public static void main(String[] args) {
        int[][] square = {{0,0}, {4,0}, {4,4}, {0,4}};
        System.out.println(rayCast(square, new int[]{2,2})); // true
        System.out.println(rayCast(square, new int[]{5,5})); // false
    }
}

Data Structures That Don't Suck For Polygon Queries

Throwing vertices into a flat array is the beginner move. In production, your polygon is rarely static — it's a stream of GPS pings, game world polygons, or GeoJSON blobs. The data structure you pick determines whether your PiP check takes 2 microseconds or 200 milliseconds.

For batch queries against the same polygon, use a spatial index. Grids are cheap: divide bounding box into cells, store which cells contain edges. For each query point, check only edges in its cell. R-trees are better for dynamic polygons — they adapt. For convex polygons, an array sorted by angle + binary search is unbeatable: O(log n) query, O(n) build.

Don't forget winding number requires storing edges in order. Ray casting doesn't care about order. If you're mixing queries with polygon editing (adding/removing vertices), use a linked list for vertices and a separate edge index. Never concatenate geometry updates with geometry queries without a dirty flag.

// io.thecodeforge — dsa tutorial

import java.util.*;

public class SpatialGridIndex {
    static class Grid {
        int cellSize;
        Map<Long, List<int[]>> cells = new HashMap<>();
        
        Grid(int[][] polygon, int cellSize) {
            this.cellSize = cellSize;
            for (int i = 0; i < polygon.length; i++) {
                int[] p = polygon[i];
                long key = (p[0] / cellSize) << 32 | (p[1] / cellSize);
                cells.computeIfAbsent(key, k -> new ArrayList<>()).add(p);
            }
        }
        
        boolean contains(int[] point) {
            long key = (point[0] / cellSize) << 32 | (point[1] / cellSize);
            List<int[]> verts = cells.get(key);
            if (verts == null) return false;
            // Use ray casting on just this cell's vertices
            // (simplified — real impl needs edge stitching)
            return verts.size() > 0;
        }
    }
    
    public static void main(String[] args) {
        int[][] poly = {{0,0}, {10,0}, {10,10}, {0,10}};
        Grid g = new Grid(poly, 5);
        System.out.println(g.contains(new int[]{3,3}));
        System.out.println(g.contains(new int[]{12,12}));
    }
}

Why Your PiP Implementation Fails in Production (And How to Fix It)

You've read the theory. You wrote ray casting. Your unit tests pass. Then a user submits a point exactly on a vertex, and your function returns false. Welcome to the real world.

Failures in production PiP aren't about the algorithm — they're about edge cases. Points on vertices, points on horizontal edges, floating-point precision at the 15th decimal. Your ray casting assumes infinite precision. The CPU gives you 53 bits. The difference between "on the edge" and "1e-16 inside" is indistinguishable to the hardware, but it's a bug to your customer.

Fix: Use epsilon comparisons. Treat any distance less than 1e-9 as "on boundary" and define behavior explicitly — include or exclude? For winding number, use integer arithmetic if coordinates are integers. For ray casting, shift the ray slightly off the horizontal (angle epsilon) to avoid parallel edge degeneracy. Log failures. Instrument every query with a check: if result flips on a microperturbation, you have a degenerate case.

Don't trust floating point. Don't trust your test data. Do trust that someone will send you a point exactly on a vertex at 3 AM.

// io.thecodeforge — dsa tutorial

public class RobustPiP {
    static final double EPS = 1e-9;
    
    static boolean isOnEdge(int[] a, int[] b, double[] p) {
        double cross = (p[0] - a[0]) * (b[1] - a[1]) - 
                       (p[1] - a[1]) * (b[0] - a[0]);
        if (Math.abs(cross) > EPS) return false;
        double dot = (p[0] - a[0]) * (b[0] - a[0]) + 
                     (p[1] - a[1]) * (b[1] - a[1]);
        if (dot < 0 || dot > (b[0]-a[0])*(b[0]-a[0]) + (b[1]-a[1])*(b[1]-a[1]))
            return false;
        return true;
    }
    
    static boolean strictInside(int[][] poly, double[] p) {
        boolean inside = false;
        for (int i = 0, j = poly.length - 1; i < poly.length; j = i++) {
            int xi = poly[i][0], yi = poly[i][1];
            int xj = poly[j][0], yj = poly[j][1];
            if (isOnEdge(poly[i], poly[j], p)) return false; // boundary = out
            if ((yi > p[1]) != (yj > p[1]) &&
                p[0] < (xj - xi) * (p[1] - yi) / (yj - yi) + xi + EPS) {
                inside = !inside;
            }
        }
        return inside;
    }
    
    public static void main(String[] args) {
        int[][] poly = {{0,0}, {5,0}, {5,5}, {0,5}};
        System.out.println(strictInside(poly, new double[]{2.5,2.5})); // true
        System.out.println(strictInside(poly, new double[]{0,0}));     // false
    }
}