Understanding douglas_simplify Algorithm Internals in C++: A Deep Dive into pybind11-rdp

The douglas_simplify family implements the Ramer-Douglas-Peucker line simplification algorithm through both recursive divide-and-conquer and iterative queue-based approaches, using analytical squared-distance calculations and Eigen matrix operations to eliminate redundant polyline points while preserving shape fidelity.

The cubao/pybind11-rdp repository exposes high-performance RDP functionality to Python via pybind11, with the computational core residing in src/main.cpp. Understanding the douglas_simplify algorithm internals reveals how the library achieves sub-millisecond performance on large coordinate arrays through careful geometric optimizations and memory-efficient mask operations.

Recursive vs. Iterative Architecture

The C++ core provides two algorithmic variants controlled by a boolean flag:

• douglas_simplify (recursive): Pure divide-and-conquer recursion that processes sub-segments via function call stack
• douglas_simplify_iter (iterative): Breadth-first queue-based traversal using std::queue<std::pair<int,int>> to avoid deep call stacks

Both variants return identical mathematical results but differ in memory characteristics. The recursive approach has $O(\log n)$ stack depth for balanced inputs, while the iterative version caps memory usage at $O(n)$ heap-allocated queue entries, making it safer for pathological polylines with millions of points.

Core Geometry: The LineSegment Class

The geometric primitive LineSegment (lines 18-44 in src/main.cpp) pre-computes vector mathematics to minimize per-point calculation overhead during the distance-scanning phase:

struct LineSegment {
    Eigen::Vector3d A, B, AB;
    double len2, inv_len2;
    
    double distance2(const Eigen::Vector3d &P) const {
        double dot = (P - A).dot(AB);
        if (dot <= 0)         return (P - A).squaredNorm();   // before A
        else if (dot >= len2) return (P - B).squaredNorm();   // after B
        return (A + (dot * inv_len2 * AB) - P).squaredNorm(); // projection
    }
};

Key optimization: The distance2 method avoids expensive square-root operations by comparing squared distances directly against $\epsilon^2$. It handles three geometric cases: points beyond endpoint A, beyond endpoint B, or projecting onto the segment interior. The pre-computed inv_len2 (inverse squared length) enables constant-time projection calculations without division operations in the hot loop.

Recursive Divide-and-Conquer Implementation

The douglas_simplify function (lines 50-84) implements the classic RDP algorithm through recursive bisection:

1. Mark endpoints i and j as kept (to_keep[i] = to_keep[j] = 1)
2. Scan intermediate points to find the maximum perpendicular distance to the segment
3. If max_dist2 <= epsilon*epsilon, collapse the entire segment
4. Otherwise, recurse on sub-segments [i, max_index] and [max_index, j]

Degenerate-case safeguard: When multiple points share the identical maximal distance (dist2 == max_dist2), the algorithm selects the pivot closest to the middle index. This deterministic tie-breaking prevents bias toward segment ends and ensures stable simplification results regardless of input ordering.

Iterative Queue-Based Processing

douglas_simplify_iter (lines 86-124) mirrors the recursive logic but replaces the call stack with an explicit std::queue<std::pair<int,int>>:

std::queue<std::pair<int,int>> q;
q.emplace(0, N-1);
while (!q.empty()) {
    auto [i, j] = q.front(); q.pop();
    // ... find max_index ...
    if (max_dist2 > eps2) {
        q.emplace(i, max_index);
        q.emplace(max_index, j);
    }
}

This breadth-first approach prevents stack overflow on adversarial inputs (e.g., sawtooth patterns where every point is a local maximum) while maintaining identical output semantics to the recursive version.

Mask Generation and Coordinate Selection

The douglas_simplify_mask function (lines 28-40) orchestrates the simplification pipeline:

1. Allocates a zero-filled Eigen::VectorXi mask vector sized to the input point count
2. Invokes either recursive or iterative helper to mark kept indices
3. Returns the boolean mask for downstream processing

Two utility functions complete the data transformation chain:

• mask2indexes (lines 42-51): Compacts the sparse boolean mask into a dense Eigen::VectorXi of kept indices
• select_by_mask (lines 60-70): Extracts rows from the original coordinate matrix using the mask, returning the simplified Eigen::Matrix

The public inline wrapper douglas_simplify (lines 73-78) simply chains select_by_mask onto the mask produced by douglas_simplify_mask, providing a zero-copy path from raw coordinates to simplified output.

Python Binding Architecture

The Python façade in src/pybind11_rdp/__init__.py (lines 25-95) dispatches to the C++ core through three public functions:

• rdp(): Generic interface with algo parameter ("rec" or "iter")
• rdp_rec(): Forces recursive mode
• rdp_iter(): Forces iterative mode

The return_mask parameter toggles between returning the simplified coordinate array or the boolean mask vector. When algo="rec", the binding passes recursive=true to douglas_simplify_mask, triggering the stack-based implementation; "iter" selects the queue-based path.

Practical Code Examples

Recursive Simplification with NumPy Arrays

import numpy as np
from pybind11_rdp import rdp

# Generate noisy line data

pts = np.column_stack((np.linspace(0, 10, 1000), 
                       np.sin(np.linspace(0, 10, 1000))))

# Simplify with epsilon=0.1 using recursive algorithm

simplified = rdp(pts, epsilon=0.1, algo="rec")
print(f"Original: {pts.shape}, Simplified: {simplified.shape}")

Underlying path: rdp() → douglas_simplify_mask(recursive=True) → douglas_simplify() → select_by_mask().

Iterative Mode with Mask Return

from pybind11_rdp import rdp
import numpy as np

pts = [[0, 0], [1, 0.1], [2, -0.05], [3, 0], [4, 0.2]]
mask = rdp(pts, epsilon=0.15, algo="iter", return_mask=True)

# mask is boolean array indicating kept points

kept_coordinates = np.asarray(pts)[mask]
print(f"Kept indices: {np.where(mask)[0]}")

Underlying path: rdp() → douglas_simplify_mask(recursive=False) → douglas_simplify_iter() → mask2indexes() (if indexes requested).

Direct C++ Integration

#include <Eigen/Core>
#include "main.cpp"  // Contains douglas_simplify definitions

int main() {
    Eigen::Matrix<double, Eigen::Dynamic, 3, Eigen::RowMajor> coords(5, 3);
    coords << 0,0,0,  1,0.1,0,  2,-0.05,0,  3,0,0,  4,0.2,0;
    
    // Simplify recursively with epsilon=0.15
    auto result = douglas_simplify(coords, 0.15, true);
    // result is RowMajor matrix containing kept points only
}

Summary

• The douglas_simplify algorithm in cubao/pybind11-rdp implements Ramer-Douglas-Peucker through mathematically equivalent recursive and iterative C++ paths
• Squared-distance calculations in LineSegment::distance2 eliminate square-root operations and enable direct comparison against $\epsilon^2$
• The recursive variant (lines 50-84) uses divide-and-conquer with deterministic middle-index tie-breaking for degenerate cases
• The iterative variant (lines 86-124) uses std::queue to process sub-segments breadth-first, preventing stack overflow on large inputs
• Mask-based architecture separates the simplification logic (douglas_simplify_mask) from coordinate extraction (select_by_mask), enabling efficient boolean indexing
• Python bindings in src/pybind11_rdp/__init__.py expose these internals through the algo parameter while handling NumPy array conversions via pybind11

Frequently Asked Questions

Why does the algorithm use squared distances instead of actual distances?

The LineSegment::distance2 method computes squared Euclidean distances to avoid expensive sqrt operations in the performance-critical inner loop. Since the algorithm only compares distances against the squared epsilon threshold (epsilon*epsilon), maintaining squared units throughout preserves correctness while significantly improving throughput on large point clouds.

When should I use the iterative algorithm over the recursive one?

Use algo="iter" when processing extremely long polylines (millions of points) or data with high-frequency oscillations that could trigger deep recursion trees. The iterative queue-based approach in douglas_simplify_iter consumes heap memory proportional to queue depth rather than call stack depth, preventing segmentation faults on pathological inputs while producing identical geometric results.

How does the algorithm handle cases where multiple points have the same maximum distance?

The implementation includes a deterministic tie-breaking mechanism in the else if (dist2 == max_dist2) block: when multiple points share the identical maximal squared distance, it selects the point closest to the middle index of the current segment range. This prevents algorithmic bias toward segment endpoints and ensures stable, repeatable simplification results regardless of input point ordering.

Can I use the C++ API directly without Python?

Yes. Include src/main.cpp in your project and link against Eigen. The public interface consists of douglas_simplify(), douglas_simplify_mask(), and douglas_simplify_indexes(), all accepting Eigen::Matrix objects with RowMajor storage. The functions return Eigen matrices or vectors, making the library suitable for embedded C++ applications requiring fast polyline decimation without Python overhead.

{
  "mcpServers": {
    "instagit": {
      "command": "npx",
      "args": ["-y", "instagit@latest"]
    }
  }
}