Point-in-Polygon (Ray Casting)

How ray casting works

You have a polygon — a closed loop of vertices — and a query point. Is the point inside the loop or outside it? For a convex blob the answer looks obvious to a human eye, but a computer sees only a list of coordinates, and a polygon can be a star, a crescent, a comb with forty teeth, or a country outline with islands and lakes. You need a rule that works on all of them.

The trick is the Jordan curve theorem, stated computationally. Pick the query point and fire a ray from it in any fixed direction — say straight to the right, toward +∞. Count how many polygon edges that ray crosses. If the count is odd, the point is inside. If it's even, it's outside. That's the whole algorithm.

The intuition is a walk. Imagine starting at infinity — unambiguously outside — and sliding back along the ray toward your point. Every edge you pass takes you from outside to inside or back again. After an odd number of these flips you've ended on the opposite parity from where you started, so you're inside. This is the even-odd rule, also called the crossing number rule.

The implementation is one loop over the polygon's edges. For each edge — the segment from vertex i to vertex i+1, wrapping the last back to the first — you ask two questions. First: does this edge straddle the ray's horizontal line, i.e. is the query point's y between the edge's two endpoint y values? Second, if it does: is the edge's intersection with that horizontal line to the right of the query point? If both are true, the ray crosses this edge; flip a boolean. After the loop, the boolean's final state is your answer.

When to use ray casting

• One-off or low-volume queries against an arbitrary polygon. Zero preprocessing means the first query is as cheap as the millionth; you pay only the O(n) edge walk.
• Polygons that change every frame — a deforming game collision hull, a hand-drawn lasso selection — where any precomputed acceleration structure would be invalidated immediately.
• Concave, star-shaped, and multi-ring polygons. Ray casting handles non-convex shapes and holes with no extra code, which simpler tests (like "is the point on the correct side of every edge") cannot.
• GIS and hit-testing — "which county contains this GPS coordinate?", "did the click land in this SVG region?" The even-odd rule is exactly what SVG's fill-rule: evenodd implements.

Reach for something else when the polygon is fixed and you'll query it millions of times: then a O(n log n) trapezoidal decomposition or a slab structure amortizes the build cost into O(log n) queries. If the polygon is convex, binary-searching its vertices gives O(log n) per query directly. And if you only need a fast reject, a bounding-box check before ray casting culls most outside points in four comparisons.

Ray casting vs other inside tests

	Ray casting (even-odd)	Winding number	Sign-of-cross-product (convex only)	Trapezoidal map	Bounding box
Query time	O(n)	O(n)	O(n) or O(log n)	O(log n)	O(1)
Preprocessing	None	None	None (sort for log n)	O(n log n)	O(n) once
Extra space	O(1)	O(1)	O(1)	O(n)	O(1)
Works on concave	Yes	Yes	No	Yes	n/a (approximate)
Handles holes	Yes (extra rings)	Yes	No	Yes	No
Self-overlap behavior	Overlap = outside	Overlap = inside	n/a	Depends on fill rule	n/a
Per-edge cost	1-2 compares + 1 divide	1 cross product + atan/turn	1 cross product	tree descent	4 compares total
Typical use	SVG evenodd, GIS, hit-test	CAD, robust fill of overlapping paths	Convex collision hulls	Repeated queries, fixed map	Broad-phase reject

The headline contrast is ray casting vs winding number. Both are O(n) and give identical answers on a simple polygon (no self-intersections). They diverge only when the boundary crosses itself: the even-odd rule cancels overlapping regions to "outside," while winding number keeps accumulating turns and reports overlaps as "inside." Most code uses ray casting because it's cheaper per edge — comparisons beat trigonometry — and because the even-odd rule is what designers expect from SVG and PostScript fills.

What the numbers actually say

• Per query: roughly n iterations, each ~2 comparisons and ≤1 division. A 1,000-vertex county outline costs about 1,000 cheap iterations — on the order of a few microseconds in compiled code. Run that for 10 million GPS points and you're at tens of seconds of pure ray casting, which is why a bounding-box pre-reject matters: it culls the ~90% of points that fall outside the box in 4 comparisons each, often a 5–10× speedup overall.
• Bounding box rejects in 4 comparisons. If px < minX || px > maxX || py < minY || py > maxY, return false immediately and skip the entire edge loop. For scattered query points this short-circuit dominates the runtime profile.
• Preprocessing pays off past the break-even point. A trapezoidal map costs O(n log n) to build but answers in O(log n). For a 1,000-vertex polygon, log₂(1000) ≈ 10 edge-equivalent operations per query versus 1,000 for ray casting — a ~100× per-query win. The build amortizes fast: ray casting costs Q·n for Q queries while preprocessing costs n·log n + Q·log n, so the break-even is Q ≈ n·log n / (n − log n) ≈ log n queries — on the order of ~10 queries on this polygon, preprocessing wins.
• Convex shortcut: O(log n) with no build. For a convex polygon stored in angular order, binary search finds the wedge containing the point's direction and one cross product confirms it — about 10 operations for n = 1,000, no preprocessing structure to maintain.

JavaScript implementation

The canonical implementation is a few lines. poly is an array of [x, y] vertices; edges wrap with the j = i trailing-index trick so the last vertex connects back to the first.

// Crossing-number ray cast with a horizontal ray to +∞.
// poly: [[x0,y0], [x1,y1], ...]; point: [px, py]
function pointInPolygon([px, py], poly) {
  let inside = false;
  for (let i = 0, j = poly.length - 1; i < poly.length; j = i++) {
    const [xi, yi] = poly[i];
    const [xj, yj] = poly[j];
    // Does edge (j → i) straddle the horizontal ray at py?
    // Half-open rule: (yi > py) !== (yj > py) counts each vertex once.
    const straddles = (yi > py) !== (yj > py);
    if (straddles) {
      // x-coordinate where edge crosses the line y = py
      const xCross = xi + (py - yi) * (xj - xi) / (yj - yi);
      if (px < xCross) inside = !inside;
    }
  }
  return inside;
}

Two details carry all the robustness. The (yi > py) !== (yj > py) test is a strict-greater-than on one end and an implicit ≤ on the other — a half-open interval that includes the lower endpoint and excludes the upper. That single asymmetry is what stops a ray grazing a shared vertex from being counted twice. And the division / (yj - yi) is only reached when the edge straddles, so yj - yi is guaranteed non-zero — horizontal edges (where it would be zero) never pass the straddle test.

A common famous-problem extension is point on a polygon with holes (a GeoJSON polygon: one outer ring plus inner hole rings). The fix is to run the same crossing count across every ring and let parity do the work:

// rings[0] is the outer boundary; rings[1..] are holes.
function pointInPolygonWithHoles(point, rings) {
  let crossings = 0;
  for (const ring of rings) {
    for (let i = 0, j = ring.length - 1; i < ring.length; j = i++) {
      const [xi, yi] = ring[i], [xj, yj] = ring[j];
      if ((yi > point[1]) !== (yj > point[1])) {
        const xCross = xi + (point[1] - yi) * (xj - xi) / (yj - yi);
        if (point[0] < xCross) crossings++;
      }
    }
  }
  return (crossings & 1) === 1;  // odd total crossings ⇒ inside
}

A point inside a hole crosses the outer ring once (odd) and the hole ring once more (now even), so it correctly reads as outside — no special hole handling needed.

Python implementation

def point_in_polygon(point, poly):
    """Crossing-number ray cast. poly: list of (x, y); point: (px, py)."""
    px, py = point
    inside = False
    n = len(poly)
    j = n - 1
    for i in range(n):
        xi, yi = poly[i]
        xj, yj = poly[j]
        if (yi > py) != (yj > py):
            x_cross = xi + (py - yi) * (xj - xi) / (yj - yi)
            if px < x_cross:
                inside = not inside
        j = i
    return inside


def on_segment(p, a, b, eps=1e-9):
    """True if p lies on segment a–b (handles the 'on the boundary' case)."""
    (px, py), (ax, ay), (bx, by) = p, a, b
    cross = (bx - ax) * (py - ay) - (by - ay) * (px - ax)
    if abs(cross) > eps:
        return False                      # not collinear
    return (min(ax, bx) - eps <= px <= max(ax, bx) + eps and
            min(ay, by) - eps <= py <= max(ay, by) + eps)

The second function answers the question the bare ray cast quietly punts on: a point lying exactly on an edge. Ray casting's verdict for a boundary point is implementation-defined — it can come back inside or outside depending on which edge the ray happens to graze. If your application needs a definite "on the boundary" answer (snapping a cursor to a wall, validating a clicked vertex), test segment membership explicitly before or after the parity test.

Variants worth knowing

Winding number. Instead of parity, sum the signed angle the polygon boundary sweeps around the point (or, more cheaply, accumulate a signed crossing count using the edge's orientation). A nonzero winding number means inside. On simple polygons it agrees with ray casting; on self-overlapping ones it reports overlap regions as inside. This is SVG's fill-rule: nonzero — the default — versus ray casting's evenodd.

Sum-of-angles winding. The textbook winding number adds up atan2 angles to each vertex. It's numerically robust and conceptually clean but slow — a transcendental call per edge. The signed-crossing variant gets the same answer with only comparisons and cross products, and is what production code uses.

Convex binary search. If the polygon is convex, store vertices in CCW order, anchor at vertex 0, and binary-search for the triangular wedge (fan slice) that contains the point's direction. One cross product confirms the point is inside that triangle. O(log n) per query, zero preprocessing beyond having the vertices sorted — which they already are for a convex polygon.

Trapezoidal decomposition / slab method. Slice the plane into vertical slabs at every vertex x-coordinate; within a slab the polygon is a stack of trapezoids you can binary-search by y. Build is O(n log n), query O(log n). This is the structure of choice when you'll hammer one fixed polygon with millions of queries.

Grid / spatial hashing. Rasterize the polygon's bounding box into a coarse grid, marking each cell inside / outside / boundary. Inside and outside cells answer in O(1); only boundary cells fall back to a ray cast. A pragmatic speedup for dense, repeated queries when an exact tree is overkill.

Common bugs and edge cases

• Double-counting vertices. If a naive test uses yi >= py && yj <= py on both ends, a ray passing through a vertex counts both adjacent edges and flips parity twice (or zero times) — wrong answer. The half-open (yi > py) !== (yj > py) rule fixes this and is the single most important line to get right.
• Horizontal edges. An edge with yi == yj has no straddle and is correctly ignored — but only if your test is the strict half-open form. A <=/>= formulation divides by zero on horizontal edges.
• Point exactly on the boundary. Ray casting gives an arbitrary answer for boundary points. If correctness on the edge matters, test segment membership separately — don't trust the parity result there.
• Forgetting to close the ring. The j = i++ wrap connects the last vertex to the first. Omit it and you leave the final edge untested, silently misclassifying points near the closing edge.
• Holes in the wrong winding for winding-number code. Ray casting (even-odd) doesn't care about ring orientation, but if you switch to the winding-number rule, hole rings must run opposite to the outer ring or holes won't subtract.
• Floating-point grazing. A query point whose y is within machine epsilon of a vertex's y can flip the straddle test unpredictably. For robustness-critical code (CAD, GIS), use integer or exact arithmetic, or snap coordinates to a fixed grid before testing.