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

> Explore the douglas_simplify algorithm internals in C++ using pybind11-rdp. Understand its recursive and iterative approaches for efficient polyline simplification with Eigen.

- Repository: [cubao/pybind11-rdp](https://github.com/cubao/pybind11-rdp)
- Tags: internals
- Published: 2026-02-28

---

**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`](https://github.com/cubao/pybind11-rdp/blob/main/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`](https://github.com/cubao/pybind11-rdp/blob/main/src/main.cpp)) pre-computes vector mathematics to minimize per-point calculation overhead during the distance-scanning phase:

```cpp
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>>`:

```cpp
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`](https://github.com/cubao/pybind11-rdp/blob/main/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

```python
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

```python
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

```cpp
#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`](https://github.com/cubao/pybind11-rdp/blob/main/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`](https://github.com/cubao/pybind11-rdp/blob/main/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.