Point in Polygon: Crossing Number and Winding Number

A polygon is a closed chain of edges, and the question of whether a point lies inside it has no answer until “inside” is defined. For a simple polygon, one that does not cross itself, the intuitive definition works and any correct test agrees with any other. For a polygon whose edges intersect, two reasonable definitions exist and they disagree on entire regions. Both are cheap to compute, both run in time linear in the number of edges, and both need the same care about rays that pass exactly through vertices.

The two definitions are the crossing number and the winding number. The crossing number counts how many times a ray from the query point crosses the boundary, and calls the point inside when that count is odd. The winding number accumulates a signed count of how many times the boundary loops around the point, and calls it inside when the total is nonzero. The first is a parity test and the second is an orientation-sensitive test, and the difference is exactly the difference between the evenodd and nonzero fill rules that 2D graphics specifications define.

Counting crossings along a horizontal ray

The classic implementation shoots a ray in the +x direction and counts edges it crosses. An edge from (x1, y1) to (x2, y2) crosses the ray when the query’s y lies in the edge’s y-range and the intersection x is greater than the query x. The range is half-open — the lower endpoint is included, the upper one is excluded — and that is what makes the test work: an edge cannot be counted twice when the ray passes through a vertex shared by two edges.

def inside_crossing(poly, p):
    px, py = p
    inside = False
    n = len(poly)
    for i in range(n):
        x1, y1 = poly[i]
        x2, y2 = poly[(i + 1) % n]
        if (y1 > py) != (y2 > py):              # edge straddles the ray line
            t = (py - y1) / (y2 - y1)
            if px < x1 + t * (x2 - x1):         # crossing lies to the right
                inside = not inside
    return inside

The half-open rule (y1 > py) != (y2 > py) is the whole trick. An edge is counted if and only if the ray’s y is in [min(y1,y2), max(y1,y2)). A vertex shared by two edges is then counted for exactly one of them, so a ray grazing a vertex contributes one crossing rather than zero or two. The same reasoning applies to a horizontal edge lying exactly on the ray: both endpoints have the same y, the comparison yields equal values, and the edge is never counted.

The same walk, with signs

The winding number replaces the toggle with an accumulator. Each upward crossing adds one, each downward crossing subtracts one, and the result is the number of times the boundary wraps around the point counterclockwise. For a simple polygon the number is 0 outside, +1 inside if the vertices are listed counterclockwise, and -1 if they are listed clockwise. Taking the absolute value, or testing for nonzero, makes the result independent of orientation.

For a convex polygon, the winding number can be computed without any ray at all, by summing the signed angles subtended by each edge, but that needs atan2 per edge and produces a floating-point total that must be compared against a tolerance. The crossing formulation is exact integer work when the inputs are integers, which is why it is preferred in practice even when the winding interpretation is what the caller wants.

Where the two definitions part company

The two rules agree on every simple polygon and disagree on some self-intersecting ones. A pentagram drawn as a single closed stroke is the standard example. The central pentagon is crossed by the boundary twice, so the crossing count there is even and the even-odd rule leaves it unfilled, while the winding number around it is two and the nonzero rule fills it. The same disagreement appears in a polygon that traces the same square twice in the same direction: inside the square the winding number is two and the crossing count is two, so the nonzero rule fills it and the even-odd rule does not. A polygon that traces the square twice in opposite directions, once in each orientation, has winding number zero inside, and the crossing count remains two: the nonzero rule leaves it unfilled because the windings cancel, and the even-odd rule leaves it unfilled because the count is even, so the two agree here.

This matters because the choice is user-visible in every vector graphics system that exposes a fill rule. The SVG specification names both nonzero and evenodd as values of the fill-rule property and defines them exactly this way, so a renderer must implement both. A geometry library that exposes only one test is making a policy decision on the caller’s behalf, and the pentagram is the shape that reveals which decision was made.

Boundary points and the tolerance question

Neither test as written has a defined answer for a point exactly on an edge. The crossing loop compares px < x_intersect and will classify an on-edge point arbitrarily depending on which side rounding lands. Applications differ on what they want: point-in-polygon for a spatial index usually wants boundary points counted as inside, while a rasteriser wants a consistent rule that assigns each boundary pixel to exactly one of two adjacent polygons, and “inside” for both would double-fill.

The usual resolution is an explicit boundary check before the main loop: for each edge, test whether the point is collinear with it and within its bounding box, using the orientation predicate from convex hulls by monotone chain. That costs one extra pass, or can be fused into the crossing loop by accumulating the cross product alongside the crossing count. With floating-point coordinates the check needs a tolerance, and the tolerance has to be chosen relative to the coordinate magnitude rather than as an absolute epsilon, because a distance of 10^-9 means different things at coordinates near 1 and coordinates near 10^9.

Cost, and what preprocessing buys

Both tests are O(n) in the number of edges per query, with no allocation beyond a few scalars. For a single query against a polygon of a few dozen vertices, that is the right implementation and anything more elaborate is wasted work. For many queries against one polygon, the linear cost becomes the bottleneck, and the fix is to index the edges by their bounding boxes so that a query only tests the edges whose boxes could contain the crossing point. A grid or a box hierarchy is the general form; a k-d tree over the edge endpoints is one concrete choice, and a uniform grid is often the faster one for polygons whose edges are evenly distributed. When the polygon is convex, containment reduces to checking that the point lies to the same side of every edge, a linear scan with no ray and no parity, and that scan can be accelerated by binary search over the edge directions to O(log n). Concave polygons do not admit that shortcut, which is one reason general polygon containment remains a linear-time operation in most libraries.

The arithmetic decides when to bother. A polygon with 10^5 edges queried 10^6 times implies 10^11 edge tests if each query is a linear scan, each a handful of comparisons, and building an index is the only way to reach that scale. Building an index costs something proportional to the edge count, so it is amortised only when the query count is large relative to the polygon size. For a coastline polygon of a few thousand vertices queried a few times per frame, the linear scan is a few thousand comparisons per query, and the setup cost of an index is never recovered. The division in the crossing test is also avoidable: comparing px against the intersection can be rewritten as a cross-product comparison against the edge’s endpoints, which keeps the arithmetic in integers and removes the only non-exact operation in the loop.

A related trap is coordinate systems. A point-in-polygon test performed on latitude and longitude treats degrees as if they were planar coordinates, and the result is wrong near the poles and across the antimeridian, where longitude wraps from 180 to -180 and a polygon spanning the date line becomes a degenerate strip. The planar test is correct only after a projection that preserves the relevant topology, and for small regions an equirectangular approximation scaled by the cosine of the latitude is usually adequate, while global queries need a spherical or ellipsoidal formulation.