Point in Polygon, in depth: crossing number, winding number and the boundary cases that break them

Two ways to ask the question

Shoot a ray from q in any fixed direction to infinity. Every time the ray crosses the polygon's boundary it switches between inside and outside, and far away it is certainly outside. So q is inside exactly when the ray crosses the boundary an odd number of times. That is the crossing-number (or even-odd, or ray-casting) test. It needs no assumption about convexity and runs in O(n) for n edges.

The winding number asks a different question: how many times does the boundary wind around q? Walk the polygon once and track the angle from q to the current point; the total change divided by a full turn is the winding number. For a simple polygon it is plus or minus one inside and zero outside, so both tests agree. They disagree on self-intersecting polygons, where a region can be wound around twice: the winding number says 2, which the nonzero rule treats as inside, while the crossing count is even, so even-odd treats it as outside. That is exactly the difference between SVG's fill-rule="nonzero" and fill-rule="evenodd", and Canvas 2D's ctx.fill("nonzero") and ctx.fill("evenodd").

The winding number does not need angles. Cast a horizontal ray to the right and, instead of counting crossings, add plus one for each edge that crosses the ray going upward and minus one for each going downward. Whether the crossing is to the right of q is decided with the orientation predicate, so the whole test uses only multiplication and subtraction, no division and no trigonometry.

Degenerate rays and the half-open rule

A ray can hit a vertex, run along a horizontal edge, or start on the boundary. A test that handles these by accident will be wrong on some input. The standard fix is the half-open rule: treat each edge as containing its lower endpoint and not its upper one. An edge from (xi, yi) to (xj, yj) is considered only if exactly one endpoint is strictly above the ray, written (yi > y) != (yj > y).

Walk through what that does. Horizontal edges have both endpoints on the same side, so they are never counted, which is right because a ray sliding along an edge does not cross it. When the ray passes through a vertex where the boundary continues upward on both sides, a local minimum, both incident edges are counted and the parity is unchanged; the ray grazes. At a local maximum neither edge counts. When the boundary passes through the vertex, one edge goes up and one comes down from the ray's level, and exactly one of them is counted. Every case comes out right without special-casing, as long as every edge applies the same rule.

The rule also gives a consistent answer for points on the boundary: a point on the left or bottom edge of a square counts as inside, on the right or top as outside. That makes adjacent polygons partition the plane with every point in exactly one, which is precisely what pixel coverage in a rasteriser needs. It is not what a geofence needs, where a point on the fence usually has to be reported as such. That is why the boundary policy must be a separate, explicit step, covered below.

Exact code for both tests

Here are both tests in Python, written so they work with integers exactly. The crossing test avoids division by comparing signs of a cross product instead of computing the intersection's x coordinate.

def orient(ax, ay, bx, by, px, py):
    # > 0 if p is left of a->b, < 0 if right, 0 if collinear
    return (bx - ax) * (py - ay) - (by - ay) * (px - ax)

def crossing_inside(poly, q):
    x, y = q
    inside = False
    n = len(poly)
    for i in range(n):
        xi, yi = poly[i]
        xj, yj = poly[(i + 1) % n]
        if (yi > y) != (yj > y):                 # half-open: edge straddles the ray
            # crossing is right of q iff q is on the correct side of the edge
            o = orient(xi, yi, xj, yj, x, y)
            if (o > 0) == (yj > yi):              # upward edge: q left; downward edge: q right
                inside = not inside
    return inside

def winding_number(poly, q):
    x, y = q
    w = 0
    n = len(poly)
    for i in range(n):
        xi, yi = poly[i]
        xj, yj = poly[(i + 1) % n]
        if yi <= y < yj and orient(xi, yi, xj, yj, x, y) > 0:
            w += 1                                # upward crossing, q on the left
        elif yj <= y < yi and orient(xi, yi, xj, yj, x, y) < 0:
            w -= 1                                # downward crossing, q on the right
    return w

def on_boundary(poly, q):
    x, y = q
    n = len(poly)
    for i in range(n):
        ax, ay = poly[i]
        bx, by = poly[(i + 1) % n]
        if (orient(ax, ay, bx, by, x, y) == 0
                and min(ax, bx) <= x <= max(ax, bx)
                and min(ay, by) <= y <= max(ay, by)):
            return True
    return False

def classify(poly, q):
    if on_boundary(poly, q):
        return "boundary"
    return "inside" if winding_number(poly, q) != 0 else "outside"

The sign trick in crossing_inside deserves a sentence. For an edge going upward, the crossing point lies to the right of q exactly when q is to the left of the edge, which is a positive orientation. For a downward edge the sides flip. So the test is that the orientation is positive precisely when the edge goes up. Both functions take O(n) time and O(1) memory and never divide, so with integer coordinates they are exact.

Worked example: a concave pentagon

Use the concave pentagon v0 = (0,0), v1 = (6,0), v2 = (6,6), v3 = (3,2), v4 = (0,6). It is a square with a V-shaped notch cut down from the top edge to v3. The shoelace sum is 48, so the signed area is +24 and the vertices run counter-clockwise.

q1 = (3, 4), inside the notch. At y = 4, the edge v2 to v3 crosses at x = 3 + (4 - 2) x 3/4 = 4.5 and the edge v3 to v4 at x = 1.5. The edge v1 to v2 crosses at x = 6 and v4 to v0 at x = 0. Crossings to the right of x = 3 are at 4.5 and 6: two, so outside. The winding test agrees: v1 to v2 is upward with q1 on its left (+1), v2 to v3 is downward with q1 on its right (-1), total 0.

q2 = (5, 4). Only the crossing at x = 6 is to the right: one crossing, inside. Winding: v1 to v2 contributes +1; v2 to v3 crosses at x = 4.5, left of q2, so q2 is on the edge's left and the downward edge does not count. Winding number +1, inside, and positive because the polygon is counter-clockwise.

q3 = (1, 2), where the ray runs exactly through v3. Apply the half-open rule edge by edge. v0 to v1 lies at y = 0, both endpoints at or below the ray: skipped. v1 to v2 goes from 0 to 6: counted, crossing at x = 6. v2 to v3 goes from 6 to 2: 6 is above and 2 is not, counted, crossing at x = 3. v3 to v4 goes from 2 to 6: counted, crossing at x = 3. v4 to v0 crosses at x = 0, left of q3. Three crossings to the right: inside, which is correct, since the notch only touches y = 2 at the single point v3. The vertex is a local minimum of the notch, so counting both of its edges is exactly the grazing behaviour the rule promises. A naive test that counted the vertex once would have said outside.

Running classify on the three queries returns outside, inside and inside, and classify(poly, (3, 2)) returns boundary because q sits on v3. Put those four cases in a unit test before touching any optimisation.

Floating point and geographic data

With floating-point coordinates, the orientation of a point very close to an edge can come out with the wrong sign, and two adjacent edges sharing a vertex can disagree about whether the ray hits it. The symptoms are points flickering in and out of a shape as you zoom, or a GPS fix on a shared border belonging to neither zone or to both. Three remedies, in order of preference:

• Snap to integers. Map coordinates to a fixed grid (for example, longitude and latitude in units of 10^-7 degrees, which fits in 32-bit integers). Watch the products: full-range longitudes in those units reach 1.8 x 10^9, a difference can reach 3.6 x 10^9, and the cross product can then overflow a signed 64-bit integer. The predicate is exact in 64 bits only while coordinate differences stay below about 2^30, for example after subtracting a local origin near the polygon; otherwise use 128-bit intermediates or arbitrary-precision integers, which Python gives you for free.
• Robust predicates. Adaptive-precision orientation tests in the style of Shewchuk's predicates evaluate in doubles and fall back to exact arithmetic only when the result is too close to zero to trust.
• Epsilon bands, carefully. Treating abs(orient) < eps as on-boundary is acceptable for UI hit-testing, where a tolerance in pixels is a feature, but eps must scale with coordinate magnitude, and it breaks the guarantee that adjacent polygons partition the plane.

Geography adds its own traps. Polygons that cross the antimeridian at 180 degrees longitude break planar tests unless you split them or shift longitudes. Edges on a map are geodesics, not straight lines in latitude and longitude, which matters for large polygons. For city-scale geofences planar math on projected coordinates is fine; for continent-scale shapes use a geodesic library rather than this code.

Holes, many polygons and many points

One point against one polygon is O(n). Real workloads are many points against many polygons, and the cost model changes. Holes first: represent a polygon with holes as an outer ring plus inner rings and run the crossing test over all rings together; even-odd parity handles holes automatically. With the winding rule, orient holes opposite to the outer ring so their contributions cancel.

The bounding-box reject is the cheapest win and should always come first: most queries in a geofencing workload are nowhere near most fences. For the many-polygon case, index the boxes in a spatial tree; the k-d tree article covers the point-index side, and the same pruning idea applies to boxes. For convex polygons the O(log n) test uses the same orientation predicate as Graham scan: find the wedge from v0 that contains q by binary search, then do one orientation test against the far edge.

Failure modes

Failure modes that show up in production code:

• Closing vertex duplicated. GeoJSON rings repeat the first vertex at the end. The modulo loop then sees a zero-length edge, harmless in the half-open crossing test but able to confuse code that divides by the edge's height. Strip it or guard for it.
• Inclusive comparisons on both endpoints. Writing yi <= y <= yj double-counts every vertex the ray passes through. This is the single most common bug.
• Computing the intersection x with division. xi + (y - yi) * (xj - xi) / (yj - yi) works in floats but adds a rounding step; with integers it truncates. The orientation form is exact.
• Wrong rule for self-intersecting input. User-drawn lassos self-intersect. Decide whether overlapping loops count as inside (nonzero) or alternate (even-odd), and match whatever renders the shape, or users will click on a filled area and miss.
• Boundary policy left implicit. If business logic cares about the fence line, call the boundary test explicitly; never rely on which way the half-open rule happens to fall.

What to do next

1. Implement orient, winding_number and on_boundary with integer coordinates and reproduce the three worked queries plus the boundary case at v3.
2. Add a randomised test: generate simple polygons and points, compare crossing parity with the winding number, and include points placed exactly on vertices and edges.
3. Decide your boundary policy (inside, outside or reported separately) and write it in the function's contract.
4. Choose the fill rule (nonzero or even-odd) to match the renderer or the data source, and test a self-intersecting shape.
5. Before scaling up, add a bounding-box reject, then a spatial index over polygon boxes if you have many polygons.
6. If coordinates are geographic, snap to an integer grid and handle the antimeridian before trusting results.