Point-in-Polygon Testing

The ray casting algorithm and spatial membership queries

Published

February 27, 2026

Before You Start

You should know
What coordinates are, and that polygons represent bounded areas made from connected edges.

You will learn
How point-in-polygon testing works, why ray casting solves the problem, and which edge cases make seemingly simple spatial queries harder than they first appear.

Why this matters
This operation sits underneath delivery zones, electoral districts, habitat boundaries, parcel assignment, and many spatial joins. It is one of the core membership tests in vector GIS.

If this gets hard, focus on…
The parity idea: each time a ray crosses a boundary, you switch between outside and inside.

Every time you ask a mapping application whether a location is within a delivery zone, a ward boundary, or a protected area, it executes a point-in-polygon test. Every time a GPS track is intersected with a land-use map to compute time spent in each zone, point-in-polygon tests run for every logged coordinate. The operation is so fundamental to GIS that it is easy to take for granted — but the algorithm that makes it work is not obvious, and getting the edge cases right (what happens when the point is exactly on the boundary?) requires careful thinking.

The ray casting algorithm provides the answer. Cast an infinite ray in any direction from the test point. Count how many times the ray crosses the polygon boundary. If the count is odd, the point is inside; if even, it is outside. The intuition: starting inside the polygon, each boundary crossing takes you outside and then back in, so an odd number of crossings means you are still inside. This topological argument — parity of crossing number — works for any polygon, however convoluted, as long as it does not self-intersect. It is one of the most elegant results in computational geometry, and implementing it robustly requires confronting exactly the kinds of degenerate cases — the ray passing through a vertex, the point on an edge — that separate working code from correct code.

1. The Question

Is this GPS coordinate inside the national park boundary?

This fundamental spatial query appears everywhere: - Wildlife ecology: Is this animal location inside its territory? - Demographics: Which census tract contains this address? - Politics: Which congressional district is this voter in? - Emergency services: Which fire station’s coverage area is closest? - Retail: Is this customer inside our delivery zone?

The mathematical question: Given a point (x, y) and a polygon defined by vertices, is the point inside the polygon?

Complications:

  • Polygons can be concave (not just rectangles)
  • Polygons can have holes
  • Edge cases: What if the point is exactly on the boundary?

2. The Conceptual Model

Intuitive Approach

Method: Cast a ray from the point to infinity. Count how many polygon edges the ray crosses.

Rule:

  • Odd number of crossings → Point is inside
  • Even number of crossings → Point is outside

Why it works:

Imagine walking from far outside the polygon toward the point: - Each time you cross a polygon edge, you alternate between inside and outside - If you cross an odd number of edges, you end up inside - If you cross an even number, you end up outside

Example:

Point A: Ray crosses 1 edge → inside
Point B: Ray crosses 2 edges → outside
Point C: Ray crosses 3 edges → inside

    ┌──────────┐
    │    A●    │ ←── 1 crossing
────┼──────────┼──→
  B●│          │     2 crossings
    │        C●│ ←── 3 crossings
    └──────────┘

Ray Direction Choice

Simplification: Cast ray horizontally to the right (positive x-direction).

Why horizontal? - Reduces 2D problem to 1D (only need to check y-coordinates) - Easier to implement - Standard convention

The test becomes: Does the horizontal ray from (x, y) cross each edge?

Ray Casting

Count Boundary Crossings, Then Check Parity

A point is inside when the ray crosses the boundary an odd number of times. Each crossing flips you from outside to inside, or back again.

test point
crossing 1
crossing 2
Cast a ray to the right from the test point.
Every boundary crossing flips your inside/outside state.
Parity rule

Odd = Inside

1 crossing, 3 crossings, 5 crossings: inside. 0, 2, 4: outside.

The ray-casting test turns a geometric membership question into a simple parity count.

3. Building the Mathematical Model

Edge Crossing Test

For each polygon edge from (x_1, y_1) to (x_2, y_2), determine if the ray from point (x_p, y_p) crosses it.

Conditions for crossing:

  1. Y-coordinate check: Point’s y must be between edge’s y-values

\min(y_1, y_2) \leq y_p < \max(y_1, y_2)

(Note: Use < for one endpoint to avoid double-counting shared vertices)

  1. X-coordinate check: Intersection point must be to the right of the point

First, find where the edge crosses the horizontal line y = y_p.

Line equation for edge:

x = x_1 + \frac{(y_p - y_1)(x_2 - x_1)}{y_2 - y_1}

This is the x-coordinate where the edge intersects the ray’s horizontal line.

Crossing occurs if:

x_{\text{intersection}} > x_p

(Ray goes to the right, so intersection must be rightward of point)

Complete Algorithm

function pointInPolygon(point, polygon):
    count = 0
    for each edge (v1, v2) in polygon:
        if (v1.y <= point.y < v2.y) or (v2.y <= point.y < v1.y):
            # Edge straddles horizontal ray
            x_intersect = v1.x + (point.y - v1.y) * (v2.x - v1.x) / (v2.y - v1.y)
            if x_intersect > point.x:
                count = count + 1
    
    return (count % 2 == 1)  # Odd = inside, Even = outside

Edge Cases

Point exactly on edge:

  • Mathematically: On boundary, could define as inside or outside
  • Convention: Usually treated as inside for GIS applications
  • Implementation: Check distance to edge separately if exactness matters

Point on vertex:

  • Same as edge case
  • Usually treated as inside

Horizontal edges:

  • Edge with y_1 = y_2 (horizontal edge)
  • Our algorithm skips these (both inequalities fail)
  • Correct behavior: Horizontal edges don’t affect inside/outside status

Degenerate polygons:

  • Self-intersecting edges
  • Zero-area polygons
  • Require preprocessing or rejection

4. Worked Example by Hand

Problem: Is point P = (3, 3) inside the quadrilateral with vertices: - A = (1, 1) - B = (5, 2) - C = (4, 5) - D = (2, 4)

Polygon edges: A \to B \to C \to D \to A

Solution

Edge 1: A(1,1) \to B(5,2)

Y-check: Is 1 \leq 3 < 2? No (3 is not less than 2)
Skip this edge.

Edge 2: B(5,2) \to C(4,5)

Y-check: Is 2 \leq 3 < 5? Yes

X-intersection:

x = 5 + \frac{(3-2)(4-5)}{5-2} = 5 + \frac{1 \times (-1)}{3} = 5 - 0.33 = 4.67

Is 4.67 > 3? Yes → Count = 1

Edge 3: C(4,5) \to D(2,4)

Y-check: Is 4 \leq 3 < 5? No (3 is not ≥ 4)
Skip this edge.

Edge 4: D(2,4) \to A(1,1)

Y-check: Is 1 \leq 3 < 4? Yes

X-intersection:

x = 2 + \frac{(3-4)(1-2)}{1-4} = 2 + \frac{(-1)(-1)}{-3} = 2 + \frac{1}{-3} = 2 - 0.33 = 1.67

Is 1.67 > 3? No (intersection is to the left)
Don’t count.

Final count: 1 (odd) → Point is INSIDE

5. Computational Implementation

Below is an interactive point-in-polygon tester.

Click on the canvas to test points

Point: (--, --)

Ray crossings: --

Result: --

Try this:

  • Click inside the polygon: Ray crosses odd number of edges → Green point (INSIDE)
  • Click outside: Ray crosses even number (including 0) → Red point (OUTSIDE)
  • Star shape: Tests concave polygon handling
  • Complex shape: Shows algorithm works for any simple polygon
  • Count crossings: Orange X marks show where ray intersects edges
  • Notice: Algorithm works regardless of polygon complexity

Key insight: The parity (odd/even) of crossings determines inside/outside, making the test robust for any polygon shape.

6. Interpretation

Why Ray Casting Works

Topological argument:

A polygon divides the plane into two regions: inside and outside. Any path from infinity to a point must cross the boundary an odd number of times if the point is inside, even number if outside.

Mathematical proof sketch:

  1. Start far outside (e.g., x = +\infty) → definitely outside
  2. Each boundary crossing toggles inside ↔︎ outside state
  3. Count crossings along path to point
  4. Odd count → ended inside; Even count → ended outside

This is Jordan Curve Theorem (simple closed curve divides plane into two regions).

Computational Complexity

Time complexity: O(n) where n is number of polygon vertices

Why: Must check every edge once

Cannot do better: Must examine every edge (counterexample: polygon could have narrow corridor)

Optimization: For many points in same polygon, precompute bounding box: - Quick rejection: If point outside bounding box → definitely outside polygon - Reduces unnecessary edge checks

Applications

Spatial joins:

Given 10,000 points and 500 polygons, assign each point to its containing polygon:

for each point:
    for each polygon:
        if pointInPolygon(point, polygon):
            assign point to polygon
            break

Optimization: Use spatial indexing (R-tree, quadtree) to avoid testing all polygons.

Zonal statistics:

Calculate average temperature within each climate zone:

for each temperature measurement point:
    for each zone polygon:
        if pointInPolygon(point, zone):
            add temperature to zone's sum
            
for each zone:
    average = sum / count

Ecological habitat analysis:

Which GPS-collared animal locations fall inside protected areas?

protected_locations = []
for each GPS fix:
    if pointInPolygon(fix, park_boundary):
        protected_locations.append(fix)

percent_protected = len(protected_locations) / total_fixes * 100

7. What Could Go Wrong?

Numerical Precision Issues

Floating-point arithmetic can cause problems:

# Edge from (0, 0) to (10, 10)
# Point at (5, 5.0000000001)

Is point exactly on edge? No (numerically)
Result: Might be classified as outside when visually "on" edge

Solution: Use tolerance for boundary tests:

if abs(distance_to_edge) < epsilon:
    treat as "on boundary" (typically: inside)

Polygon with Holes

Standard ray casting doesn’t handle holes.

Example: Donut-shaped polygon (outer ring + inner hole)

Solution: Track polygon as multiple rings: - Outer ring: clockwise vertices - Holes: counter-clockwise vertices

Algorithm modification:

count = 0
for each ring in polygon:
    crossings = ray_cast(point, ring)
    if ring is outer:
        count += crossings
    else:  # hole
        count -= crossings

return (count % 2 == 1)

Self-Intersecting Polygons

Example: Figure-8 or bowtie shape

Ray casting still works but definition of “inside” becomes ambiguous.

Typical behavior: Uses winding number (how many times polygon winds around point).

Our algorithm: Treats each crossing independently (may not match visual intuition for complex self-intersections).

Degenerate Cases

Zero-area polygons:

  • All vertices collinear → no interior
  • Treat all points as outside

Duplicate vertices:

  • Clean up polygon beforehand
  • Remove consecutive duplicate points

Nearly-horizontal edges:

If edge has y_1 \approx y_2:

\frac{(x_2 - x_1)}{(y_2 - y_1)} \to \text{very large or undefined}

Solution: Skip edges with |y_2 - y_1| < \varepsilon (effectively horizontal).

8. Extension: Winding Number Algorithm

Alternative approach: Count signed crossings (direction matters).

Winding number = net number of times polygon winds counter-clockwise around point.

Formula:

w = \frac{1}{2\pi} \oint_{\partial P} d\theta

Where \theta is angle from point to polygon boundary.

Discrete version:

For each edge, calculate change in angle from point to edge:

winding_number = 0
for each edge (v1, v2):
    angle1 = atan2(v1.y - point.y, v1.x - point.x)
    angle2 = atan2(v2.y - point.y, v2.x - point.x)
    winding_number += angle_difference(angle1, angle2)

return (winding_number != 0)

Advantages:

  • Handles self-intersecting polygons better
  • More geometrically meaningful

Disadvantages:

  • Slower (requires trigonometric functions)
  • More complex implementation

Ray casting is standard for simple polygons (faster, simpler).

9. Math Refresher: Line Intersection

Finding X-Coordinate of Intersection

Given edge from (x_1, y_1) to (x_2, y_2) and horizontal line y = y_p:

Line equation (parametric form):

x(t) = x_1 + t(x_2 - x_1) y(t) = y_1 + t(y_2 - y_1)

Where t \in [0, 1] (on edge).

Find t where y(t) = y_p:

y_1 + t(y_2 - y_1) = y_p

t = \frac{y_p - y_1}{y_2 - y_1}

Substitute into x equation:

x = x_1 + \frac{y_p - y_1}{y_2 - y_1}(x_2 - x_1)

This is the intersection x-coordinate used in the algorithm.

Why Use < Instead of \leq?

Edge from (0, 0) to (0, 10):

Two edges meet at (0, 10).

If we use \leq for both endpoints: - First edge: 0 \leq y_p \leq 10 → counts crossing at y_p = 10 - Second edge: 10 \leq y_p \leq 20 → counts crossing at y_p = 10

Double counting! Same vertex counted twice.

Solution: Use y_1 \leq y_p < y_2 (include lower, exclude upper).

Result: Shared vertices counted exactly once.

Summary

  • Point-in-polygon test determines if point is inside polygon boundary
  • Ray casting algorithm: Cast ray to infinity, count edge crossings
  • Odd crossings → inside; Even crossings → outside
  • Complexity: O(n) where n = number of polygon vertices
  • Implementation: Check y-range, calculate x-intersection, count if right of point
  • Works for any simple polygon (convex or concave)
  • Edge cases: Handle boundaries, holes, and degenerate polygons carefully
  • Applications: Spatial joins, zonal statistics, habitat analysis, demographic queries
  • Alternative: Winding number algorithm for self-intersecting polygons