Implementing polygon containment checks in C++

This page solves one concrete problem: deciding, deterministically and allocation-free, whether a streaming GPS fix falls inside a predefined geofence polygon on a constrained edge gateway — a header-only C++ predicate that runs on everything from an ARM Cortex-M7 microcontroller to an x86_64 industrial gateway. Within the Local Spatial Processing Patterns practice, and specifically as a concrete build of on-device geometry filtering, the containment check is the precise second stage that runs only on points a cheap bounding-box reject has already let through. Pushing this decision to the network edge eliminates cloud round-trip latency, spares metered backhaul, and keeps the geofence verdict available even when the uplink has been gone for hours.

The deployment context is unforgiving: field nodes cannot afford dynamic allocation, garbage-collection pauses, or unpredictable cache behavior while a modem and a sensor bus contend for the same core. The implementation below is therefore built to be allocation-free, branch-predictable, and tolerant of raw GNSS drift — the three properties that separate a containment check that survives a fanless gateway in the field from one that only ever ran on a workstation.

Algorithm selection rationale

The ray-casting algorithm (even-odd rule) remains the right choice for on-device polygon containment because its cost profile maps cleanly onto the constraint envelope: O(n) in the vertex count, with minimal state and a strictly sequential memory access pattern. Unlike winding-number approaches — which track directional crossings and carry more topological bookkeeping — ray-casting only counts edge intersections along a single horizontal projection. On embedded Linux or a bare-metal RTOS that translates directly into deterministic cache behavior and zero heap fragmentation, which is exactly what a per-packet hot path needs.

Three deployment constraints are baked into the design rather than bolted on afterward:

  1. Zero dynamic allocation. All vertex data is passed via raw pointers or fixed-size arrays. No STL containers, no heap allocation, no exceptions — the entire working set is the polygon, the cached bounding box, and a handful of stack locals.
  2. Early-exit bounding box. A cheap AABB (axis-aligned bounding box) check rejects the overwhelming majority of out-of-bounds points — typically 90%+ during real telemetry — before the intersection loop ever runs, preserving CPU cycles for high-frequency ingestion.
  3. Boundary tolerance. GNSS jitter rarely lands exactly on a vertex or edge. A configurable epsilon prevents false negatives caused by coordinate quantization and atmospheric multipath, the same drift-tolerance idea applied elsewhere to movement and threshold gating for sensor events.
Two-stage geofence test: AABB early-exit then ray-casting A concave geofence polygon shaped like a U sits inside its axis-aligned bounding box, drawn as a dashed rectangle. Point A lies to the left of the bounding box and is rejected in constant time before any loop runs. Point B lies inside the bounding box but in the notch of the U, so a horizontal ray cast to the right crosses two polygon edges; an even crossing count means the point is outside and it is rejected by ray-casting. Point C lies inside both the bounding box and the solid part of the polygon and is accepted. Two-stage reject: AABB early-exit, then ray-casting AABB (cached bounding box) geofence polygon 2 crossings (even) — outside A B C A Outside AABB rejected in O(1) — four float compares, no loop B Inside AABB, outside polygon rejected by ray-casting — even crossing count C Inside both accepted — odd crossing count, within epsilon Horizontal ray counts edges crossed to the right of the point

Ray-casting containment with an AABB early-exit and vertex tolerance.

Control flow of the point-in-polygon predicate The predicate takes a point, polygon and bounding box. It first tests whether the point is inside the AABB plus epsilon; if not, it returns false immediately. If it is inside, the loop walks each edge pair j and i. For each edge it checks whether the point sits on a vertex within epsilon, returning true if so. Otherwise it runs the crossing test, toggling the inside flag when the test is true, then advances to the next edge. When edges remain the loop repeats, and when they are exhausted the function returns the accumulated inside flag. point_in_polygon — per-call control flow no yes yes no yes no no yes Point + polygon + AABB Inside AABB + ε? Return false Loop edges j…i On a vertex? Return true Crossing test true? Toggle inside flag Next edge More edges? Return inside flag

A note on coordinate systems before any code runs: ray-casting assumes a planar Cartesian space. For geofences under roughly 5 km radius, treating WGS84 lat/lon as planar introduces negligible error (under 0.1 m at mid-latitudes). Larger boundaries must be pre-projected to a local UTM zone or EPSG:3857 at ingestion — the projection trade-offs for low-memory nodes are covered under coordinate reference systems at the edge, and the Cortex-M specifics under CRS transformations on ARM Cortex-M devices.

The complete, self-contained implementation

The following header-only implementation is designed for direct inclusion in gateway firmware. Compile with -O2 -fno-exceptions -fno-rtti for deterministic latency on constrained targets. Avoid -ffast-math: it can reorder floating-point operations in ways that break the strict inequality guard that keeps the division safe.

// geofence_containment.hpp
#pragma once
#include <cmath>
#include <cstddef>

struct Point2D {
    double lat;
    double lon;
};

struct Polygon {
    const Point2D* vertices;
    std::size_t count;
};

struct AABB {
    double min_lat, max_lat;
    double min_lon, max_lon;
};

// Precompute AABB for early-exit filtering.
// Call once during geofence provisioning, not per-telemetry packet.
inline AABB compute_aabb(const Polygon& poly) {
    AABB box{poly.vertices[0].lat, poly.vertices[0].lat,
             poly.vertices[0].lon, poly.vertices[0].lon};
    for (std::size_t i = 1; i < poly.count; ++i) {
        if (poly.vertices[i].lat < box.min_lat) box.min_lat = poly.vertices[i].lat;
        if (poly.vertices[i].lat > box.max_lat) box.max_lat = poly.vertices[i].lat;
        if (poly.vertices[i].lon < box.min_lon) box.min_lon = poly.vertices[i].lon;
        if (poly.vertices[i].lon > box.max_lon) box.max_lon = poly.vertices[i].lon;
    }
    return box;
}

// Deterministic point-in-polygon check with epsilon tolerance and AABB early-exit.
// Division by zero is mathematically impossible due to the strict inequality guard.
[[nodiscard]] inline bool point_in_polygon(const Point2D& pt,
                                           const Polygon& poly,
                                           const AABB& box,
                                           double epsilon = 1e-9) {
    // Early-exit AABB check (expanded by epsilon for boundary tolerance)
    if (pt.lat < box.min_lat - epsilon || pt.lat > box.max_lat + epsilon ||
        pt.lon < box.min_lon - epsilon || pt.lon > box.max_lon + epsilon) {
        return false;
    }

    bool inside = false;
    for (std::size_t i = 0, j = poly.count - 1; i < poly.count; j = i++) {
        double yi = poly.vertices[i].lat, xi = poly.vertices[i].lon;
        double yj = poly.vertices[j].lat, xj = poly.vertices[j].lon;

        // Vertex proximity check (handles GPS quantization on corners)
        if (std::abs(pt.lat - yi) < epsilon && std::abs(pt.lon - xi) < epsilon) {
            return true;
        }

        // Ray-casting intersection logic
        // (yi > pt.lat) != (yj > pt.lat) ensures yi != yj, preventing division by zero
        if (((yi > pt.lat) != (yj > pt.lat)) &&
            (pt.lon < (xj - xi) * (pt.lat - yi) / (yj - yi) + xi)) {
            inside = !inside;
        }
    }
    return inside;
}

The compute_aabb call belongs in the provisioning path, not the hot path: derive the box once when a geofence is loaded and cache it in read-only memory. Every subsequent point_in_polygon call then pays only the four float comparisons of the reject, and the loop runs solely for the survivors.

Constraint validation table

Each design choice in the code maps to a specific hardware limit and the mitigation it carries. Size the pipeline against the real budget of the target — the broader picture is in device constraints and resource limits; the table below is the subset that bears directly on the containment predicate.

Constraint Expected impact Mitigation built into the code
RAM A heap-backed geometry library can spike tens of MB per polygon on a 128–256 MB node Raw-pointer Polygon, fixed AABB, only stack locals — zero allocation per call
CPU / thermal Per-point topology on every fix saturates a throttling core during 50–100 Hz bursts AABB early-exit discards 90%+ before the O(n) loop; loop body is branch-light
Latency Heap churn and GC-style pauses push tail latency past the polling interval Allocation-free, exception-free path keeps per-call latency sub-millisecond and bounded
Power Software-emulated double math burns cycles (and battery) on single-precision FPUs Optional float quantization; AABB reject keeps the expensive division off the common path

Gotchas and edge cases

The traps here are floating-point and coordinate-system traps, not algorithmic ones. The crossing test (yi > pt.lat) != (yj > pt.lat) is doing double duty: it both selects edges that straddle the ray and guarantees yi != yj, so the /(yj - yi) division can never see a zero denominator. Re-deriving that guard incorrectly — for example switching > to >= — reintroduces the division-by-zero on horizontal edges that the strict inequality is there to prevent.

  • GNSS drift tolerance. A receiver parked next to a boundary oscillates by a few meters on multipath alone. Pick epsilon from your precision contract, not arbitrarily: 1.5e-5 degrees is roughly 1.5 m at the equator. Quantizing to integer microdegrees first, per spatial data precision standards, turns the boundary comparison into exact integer math and sidesteps float-equality entirely.
  • Coordinate system assumptions. The planar approximation only holds for small geofences. Mixing a WGS84 polygon with already-projected fixes silently produces wrong verdicts; pin one reference system at ingestion and reject anything that does not match it.
  • Compiler flag pitfalls. -ffast-math is the classic foot-gun — it lets the optimizer assume no special values and reassociate the multiply/divide, which can move a point across the boundary by an ULP at exactly the edge case epsilon exists to cover. Keep it off and prefer explicit quantization for reproducibility.
  • RTOS quirks. On a single-precision FPU (Cortex-M4, -mfpu=fpv4-sp-d16) the double fields fall back to slow software emulation. Either target a double-precision FPU such as the Cortex-M7’s -mfpu=fpv5-d16, or quantize coordinates to float. Align vertex arrays to a cache line (alignas(64)) to avoid false sharing and pipeline stalls during the loop.
  • Concave boundaries. Highly concave geofences make inside = !inside toggle often, raising branch-misprediction rates. Watch them with hardware performance counters and, during provisioning, pre-flatten complex boundaries into convex pieces or a triangulated mesh to lower instruction-cache pressure.

Integrating with the gateway pipeline

In the telemetry loop, the predicate sits behind the parser and ahead of the alerting path. Feed raw NMEA or binary fixes directly into Point2D structs — no intermediate string parsing — and reuse the cached AABB:

// High-frequency telemetry loop integration
void process_geofence_stream(const Point2D* telemetry_buffer,
                             std::size_t buffer_len,
                             const Polygon& geofence,
                             const AABB& geofence_aabb) {
    for (std::size_t i = 0; i < buffer_len; ++i) {
        if (point_in_polygon(telemetry_buffer[i], geofence, geofence_aabb, 1.5e-5)) {
            // ~1.5e-5 epsilon ≈ 1.5 m at equator
            trigger_edge_alert(telemetry_buffer[i]);
        }
    }
}

For Python-fronted gateways, do not serialize coordinates as JSON across the boundary. Pack lat/lon pairs into fixed-width binary frames (struct.pack('<dd', lat, lon)) and hand the buffer to the compiled predicate through ctypes or a memory view, which keeps the ingress path allocation-free and the GIL out of the hot loop. The same fixed-width framing is what makes the downstream delta sync for GPS coordinate streams cheap once a fix has cleared the geofence and earned its place on the uplink.

Validate against field conditions before trusting the epsilon: feed coordinates oscillating ±2 m around polygon edges to confirm the tolerance behaves, and cross-reference verdicts against the OGC Simple Features Specification test suites when migrating to multi-polygon or hole-containing geometries.