The ray casting algorithm and spatial membership queries
2026-02-27
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.
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?
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
└──────────┘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?
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:
\min(y_1, y_2) \leq y_p < \max(y_1, y_2)
(Note: Use < for one endpoint to avoid double-counting shared vertices)
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)
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 = outsidePoint 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
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
Edge 1: A(1,1) \to B(5,2)
Y-check: Is $1 < 2$? No (3 is not less than 2)
Skip this edge.
Edge 2: B(5,2) \to C(4,5)
Y-check: Is $2 < 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 < 5$? No (3 is not ≥ 4)
Skip this edge.
Edge 4: D(2,4) \to A(1,1)
Y-check: Is $1 < 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
Below is an interactive point-in-polygon tester.
<label>
Test shape:
<select id="shape-select">
<option value="square">Square</option>
<option value="triangle">Triangle</option>
<option value="concave" selected>Concave Polygon</option>
<option value="star">Star (Concave)</option>
<option value="complex">Complex Shape</option>
</select>
</label>
<label>
Show ray:
<input type="checkbox" id="show-ray" checked>
</label>
<label>
Show intersections:
<input type="checkbox" id="show-intersections" checked>
</label>
<div class="pip-info">
<p><strong>Click on the canvas to test points</strong></p>
<p>Point: <span id="point-coords">(--, --)</span></p>
<p>Ray crossings: <span id="crossing-count">--</span></p>
<p>Result: <span id="pip-result">--</span></p>
</div><canvas id="pip-canvas" width="600" height="500" style="border: 1px solid #ddd; cursor: crosshair;"></canvas>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.
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:
This is Jordan Curve Theorem (simple closed curve divides plane into two regions).
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
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
breakOptimization: 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 / countEcological 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 * 100Floating-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" edgeSolution: Use tolerance for boundary tests:
if abs(distance_to_edge) < epsilon:
treat as "on boundary" (typically: inside)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)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).
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).
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).
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.
Edge from (0, 0) to (0, 10):
Two edges meet at (0, 10).
If we use \leq for both endpoints: - First edge: $0 y_p $ → counts crossing at y_p = 10 - Second edge: $10 y_p $ → 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.