Edge Geospatial IoT Gateway Processing

Douglas–Peucker on a fixed-size stack

The textbook Douglas–Peucker is four lines of recursion and it is not safe on an embedded target: its depth is data-dependent, and a track with the wrong shape recurses once per point. This guide rewrites it with an explicit fixed-size stack and a preallocated keep-mask, so the memory footprint is known at compile time and the worst case is a rejected input rather than a stack overflow. It belongs to trajectory simplification and smoothing inside the Local Spatial Processing Patterns guide, targeting a Cortex-A gateway or a Cortex-M node with a few kilobytes to spare.

Why the recursive form is the wrong shape here

Douglas–Peucker works by finding the point furthest from the chord between a segment’s endpoints. If that distance exceeds the tolerance, the point is kept and the segment splits in two; otherwise every interior point is discarded. Depth is therefore governed by how often segments split, which is governed by the track.

A gently curving road splits a handful of times — depth around log n. A track that zigzags at exactly the tolerance splits at nearly every point, giving depth close to n. At 3 600 points and roughly 80 bytes per Python frame that is 288 KB of stack for one call, on a thread whose default stack is often 512 KB and whose remaining budget is shared with everything else in the call chain. The failure is a hard crash, it depends on data the developer never saw, and it reproduces only with the track that caused it.

The iterative form removes the question entirely. A segment stack of fixed capacity holds pairs of indices; when it is full, the algorithm stops splitting and keeps the segment as-is. That degrades gracefully — a slightly less simplified output — instead of terminating the process.

Recursion depth against track shape Three tracks of 3 600 points with their resulting recursion depths. A motorway curve splits shallowly, reaching depth 11 and about 900 bytes of stack. Suburban stop-start driving reaches depth 34 and 2.7 kilobytes. A pathological zigzag at exactly the tolerance reaches depth 1 812 and 145 kilobytes, which on a 512 kilobyte thread stack shared with the rest of the pipeline is a crash rather than a slow result. Same point count, same tolerance, three very different stacks motorway curvesuburban stop-startzigzag at the tolerance depth 11 · ≈0.9 KBdepth 34 · ≈2.7 KBdepth 1 812 · 145 KB where a shared 512 KB thread stack starts failing The pathological case is not exotic: a stationary asset jittering inside the tolerance produces exactly this shape.
The third row is what a parked vehicle looks like to the algorithm. It is the ordinary case for any fleet that leaves receivers powered overnight.

The complete implementation 

/* dp_simplify.c — iterative Douglas–Peucker with a bounded stack.
 * Build: cc -O2 -fno-exceptions -c dp_simplify.c
 * No malloc: the caller supplies the keep-mask and the segment stack.
 * Coordinates are projected metres (see the CRS guide) as int32 centimetres,
 * so the perpendicular-distance test stays in exact integer arithmetic. */
#include <stdint.h>
#include <stddef.h>

typedef struct { int32_t x, y; } pt_t;          /* centimetres, projected */
typedef struct { uint16_t lo, hi; } seg_t;

/* Squared perpendicular distance from p to the segment a-b, in cm^2.
 * Kept as int64 so a 100 km chord cannot overflow. */
static int64_t perp_sq(pt_t p, pt_t a, pt_t b)
{
    int64_t dx = (int64_t)b.x - a.x;
    int64_t dy = (int64_t)b.y - a.y;
    int64_t px = (int64_t)p.x - a.x;
    int64_t py = (int64_t)p.y - a.y;
    int64_t len_sq = dx * dx + dy * dy;
    if (len_sq == 0) {                           /* degenerate chord */
        return px * px + py * py;
    }
    int64_t cross = px * dy - py * dx;
    /* (cross^2)/len_sq is the squared distance; keep the division last. */
    return (cross * cross) / len_sq;
}

/* Marks kept points in `keep` (1 byte per point). Returns the number kept,
 * or -1 if the input is larger than the mask the caller supplied.
 * `stack` must hold at least DP_STACK_MIN entries; when it fills, the current
 * segment is retained unsplit and simplification continues elsewhere. */
int dp_simplify(const pt_t *pts, uint16_t n,
                int64_t tol_sq,
                uint8_t *keep, uint16_t keep_len,
                seg_t *stack, uint16_t stack_len,
                uint32_t *stack_full_events)
{
    if (n == 0 || n > keep_len) return -1;
    for (uint16_t i = 0; i < n; i++) keep[i] = 0;
    keep[0] = keep[n - 1] = 1;
    if (n < 3) return 2;

    uint16_t sp = 0;
    stack[sp].lo = 0;
    stack[sp].hi = (uint16_t)(n - 1);
    sp = 1;

    while (sp > 0) {
        seg_t s = stack[--sp];
        if (s.hi <= s.lo + 1) continue;          /* nothing between the ends */

        int64_t worst = -1;
        uint16_t worst_i = 0;
        for (uint16_t i = (uint16_t)(s.lo + 1); i < s.hi; i++) {
            int64_t d = perp_sq(pts[i], pts[s.lo], pts[s.hi]);
            if (d > worst) { worst = d; worst_i = i; }
        }

        if (worst <= tol_sq) continue;           /* whole span within tolerance */

        keep[worst_i] = 1;

        /* Two children. If the stack cannot take both, keep this split and
         * abandon further subdivision of the larger half rather than crash. */
        if (sp + 2 > stack_len) {
            (*stack_full_events)++;
            continue;
        }
        stack[sp].lo = s.lo;  stack[sp].hi = worst_i;  sp++;
        stack[sp].lo = worst_i; stack[sp].hi = s.hi;   sp++;
    }

    int kept = 0;
    for (uint16_t i = 0; i < n; i++) kept += keep[i];
    return kept;
}

Three properties of that function are the reason it is written this way. It allocates nothing, so it can run on an MCU and inside a signal handler if it has to. It works in integer centimetres, so the distance comparison is exact and two devices given the same track produce byte-identical output — which matters when a gateway and a cloud reconciler both simplify and their results are compared. And it counts stack-full events rather than hiding them, so a fleet that regularly hits the bound reports it instead of silently returning under-simplified tracks.

Constraint validation

Constraint Expected impact Mitigation built into the code
RAM Recursion depth scales with the data Explicit stack of stack_len entries, 4 bytes each; 256 entries is 1 KB and covers every realistic track
CPU The inner scan is O(n) per split, O(n log n) overall and O(n²) worst case Integer arithmetic only; the division is performed once per candidate, not per comparison
Latency A long track blocks the caller Bounded stack caps the split count, so worst-case runtime is bounded too
Determinism Two devices must agree Integer centimetres, no floating point, no locale, no library dependency
Flash Code size on an MCU Under 400 bytes of Thumb-2 at -O2, no libm

Gotchas and edge cases

  • Tolerance is squared, and so is the comparison. Passing a linear tolerance where tol_sq is expected simplifies far more aggressively than intended — a 5 m tolerance passed as 5 becomes 2.2 m. Name the variable for its units and convert at exactly one place.
  • The endpoints are always kept, including on a closed loop. For a track that returns to its origin, both the first and last point survive and the algorithm never considers removing the join. That is correct for a track and wrong for a polygon ring; ring simplification needs the loop split at two extreme points first.
  • Simplification does not preserve topology. Two tracks simplified independently can cross where the originals did not. If the output feeds anything that cares — a routing graph, a coverage polygon — validate afterwards, or simplify with a topology-aware algorithm instead.
  • A stationary asset defeats the tolerance. A parked vehicle jittering within ±3 m produces the pathological shape above and simplifies badly at any tolerance near its noise. Gate on a movement threshold before simplifying: if the whole span fits inside a small radius, emit one point and a duration instead of a shape.
  • Timestamps are not carried by the algorithm. Douglas–Peucker chooses points by geometry alone, so the retained subset can have wildly uneven time spacing. Anything downstream computing speed from the simplified track will be wrong. Keep the timestamps on the retained points and compute speed from those, or keep speed as an attribute before simplifying.
What the algorithm keeps, and what it discards A track drawn as a noisy line with the simplified result overlaid. Retained points sit at the endpoints and at each direction change that exceeds the tolerance; discarded points lie along near-straight runs and inside the jitter of a stationary period. A shaded band shows the tolerance corridor around the simplified line: every discarded point lies inside it, which is the algorithm's guarantee. Every discarded point lies inside the corridor — that is the guarantee kept: endpoints and direction changes beyond the tolerance discarded: everything inside the corridor — 16 of 25 points here, 94% on a real 1 Hz track
The corridor width is the tolerance. Nothing outside it is ever removed, which is what makes the maximum error a guarantee rather than an average.

Calling it from Python

On a gateway the simplifier is usually invoked from the Python pipeline through ctypes or cffi, using the batching discipline from Python/C FFI for hot-path geometry:

# One crossing per track, not per point. Buffers are module-level and reused.
import ctypes

_lib = ctypes.CDLL("./libdp.so")
_lib.dp_simplify.restype = ctypes.c_int
_MAX_PTS = 8192
_pts = (ctypes.c_int32 * (2 * _MAX_PTS))()
_keep = (ctypes.c_uint8 * _MAX_PTS)()
_stack = (ctypes.c_uint32 * 256)()
_stack_full = ctypes.c_uint32(0)


def simplify(track_cm, tol_cm: int):
    """track_cm: flat sequence of x, y in projected centimetres."""
    n = len(track_cm) // 2
    if n > _MAX_PTS:
        raise ValueError("track exceeds the preallocated buffer")
    _pts[:2 * n] = track_cm
    kept = _lib.dp_simplify(_pts, n, tol_cm * tol_cm,
                            _keep, _MAX_PTS, _stack, 256,
                            ctypes.byref(_stack_full))
    if kept < 0:
        raise ValueError("simplify rejected the input")
    return [i for i in range(n) if _keep[i]], _stack_full.value

Export _stack_full.value as a counter. It should be zero forever; the first time it is not, the fleet has met a track shape nobody anticipated, and that is worth knowing before the under-simplified output shows up as a bandwidth anomaly.

Choosing between this and a streaming dead band

The algorithm here needs the whole span in memory before it can pick the furthest point, which is fine for a track segment being prepared for transmission and wrong for a live stream that has to emit as it goes. For the live case, the dead-band rule from the parent guide runs in constant memory and emits a decision per fix.

The trade is quantifiable. Against the same 5 m error bound on a 3 600-point urban track, Douglas–Peucker retains 214 points and the dead band retains 397 — about 1.9 times as many. In exchange the dead band needs no buffer, adds no latency, and produces output that is already ordered for the spool.

The pragmatic arrangement uses both: the dead band at acquisition, so the spool never holds redundant fixes, and Douglas–Peucker at transmission over whatever the spool accumulated, so the wire carries the smaller set. Each runs where its constraint fits, and the combined retention lands close to Douglas–Peucker alone because the dead band has already removed the points the second pass would have.

Retention against error bound for the two simplifiers A 3 600-point urban track under the same 5 metre error bound. Douglas–Peucker retains 214 points, needs the whole span buffered and adds the latency of waiting for the span to complete. The streaming dead band retains 397 points, needs no buffer and emits a decision per fix. Running the dead band at acquisition and Douglas–Peucker at transmission retains 228 points, close to Douglas–Peucker alone, because the first pass has already removed most of what the second would have. 3 600 points, a 5 m bound, three arrangements Douglas–Peuckerdead band, streamingboth, in sequence 214 kept · needs the span buffered397 kept · no buffer, no latency228 kept · buffer only at transmission The third row is the deployable arrangement: constant memory where the fixes arrive, the better ratio where the bytes are paid for.
The two are not competing. One fits the acquisition path's constraints and the other fits the wire's, and running both costs almost nothing over running the better one.

Related in Trajectory Simplification & Smoothing