# How pybind11-rdp Achieves 8000x Speedup Over the Python rdp Library

> Discover how pybind11-rdp achieves an 8000x speedup over the Python rdp library. Learn about C++ implementation Eigen vectorization and pybind11 bindings for unparalleled performance.

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

---

**pybind11-rdp delivers approximately 8,000 times faster performance than the pure-Python rdp library by implementing the Ramer-Douglas-Peucker algorithm in compiled C++ with Eigen vectorization, eliminating Python interpreter overhead through pybind11 bindings, and using zero-copy NumPy memory access.**

The `cubao/pybind11-rdp` repository reimplements the classic polyline simplification algorithm using modern C++17, leveraging the Eigen library for SIMD-optimized linear algebra and pybind11 for seamless Python interoperability. This architectural shift moves the computational bottleneck from interpreted Python loops to native machine code, transforming a O(N log N) algorithm that struggles with thousand-point datasets into one that processes millions of points in milliseconds.

## Architectural Deep-Dive: From Python Loops to Compiled C++

The performance gap stems from fundamental architectural differences between the reference Python implementation and the pybind11-rdp C++ core. While the pure-Python library iterates through points with interpreter overhead for every distance calculation, pybind11-rdp executes tight loops in compiled code with vectorized SIMD instructions.

### Zero-Copy NumPy Bridge with Eigen

In [`src/main.cpp`](https://github.com/cubao/pybind11-rdp/blob/main/src/main.cpp), the binding layer uses `Eigen::Ref<const RowVectors>` to accept NumPy arrays directly:

```cpp
using RowVectors = Eigen::Matrix<double, Eigen::Dynamic, 3, Eigen::RowMajor>;
void rdp(Eigen::Ref<const RowVectors> coords, double epsilon, bool recursive);

```

This `Eigen::Ref` adapter creates a zero-copy view into the NumPy memory buffer. The C++ algorithm reads coordinates directly from Python-managed memory without allocation overhead or data marshalling costs, a critical optimization when processing large point clouds.

### Pre-Computed Geometry with LineSegment

The `LineSegment` struct in [`src/main.cpp`](https://github.com/cubao/pybind11-rdp/blob/main/src/main.cpp) (lines 18-44) eliminates redundant arithmetic by pre-computing segment invariants:

```cpp
struct LineSegment {
    Eigen::Vector3d A, B, AB;
    double len2, inv_len2;
    
    LineSegment(const Eigen::Vector3d& a, const Eigen::Vector3d& b) 
        : A(a), B(b), AB(b - a) {
        len2 = AB.squaredNorm();
        inv_len2 = len2 == 0 ? 0. : 1. / len2;
    }
    
    double distance2(const Eigen::Vector3d& P) const {
        // Optimized point-to-segment distance using pre-computed values
        double t = (P - A).dot(AB) * inv_len2;
        t = std::max(0., std::min(1., t));
        return (A + t * AB - P).squaredNorm();
    }
};

```

By storing `AB` (the vector from A to B), `len2` (squared length), and `inv_len2` (inverse squared length), the distance calculation reduces to a few floating-point operations. This avoids the repeated subtraction and norm calculations that would occur in a naive implementation.

### Eliminating Interpreter Overhead

The pure-Python rdp library invokes Python-level distance functions within nested loops. In `pybind11-rdp`, the entire simplification algorithm executes as a single C++ function call:

```cpp
std::vector<bool> douglas_simplify_iter(
    Eigen::Ref<const RowVectors> coords,
    double epsilon_sq
) {
    std::vector<bool> keep(coords.rows(), false);
    std::queue<std::pair<int, int>> queue;
    queue.emplace(0, coords.rows() - 1);
    keep[0] = keep[coords.rows() - 1] = true;
    
    while (!queue.empty()) {
        auto [i, j] = queue.front();
        queue.pop();
        
        // Tight C++ loop - no Python interpreter involvement
        double max_dist = 0;
        int max_idx = i;
        LineSegment seg(coords.row(i), coords.row(j));
        
        for (int k = i + 1; k < j; ++k) {
            double d = seg.distance2(coords.row(k));
            if (d > max_dist) {
                max_dist = d;
                max_idx = k;
            }
        }
        
        if (max_dist > epsilon_sq) {
            keep[max_idx] = true;
            queue.emplace(i, max_idx);
            queue.emplace(max_idx, j);
        }
    }
    return keep;
}

```

This implementation uses `std::queue` for the iterative variant, avoiding recursion depth limits while maintaining O(N log N) complexity. The inner loop performs distance calculations using pre-computed segment data with no Python API calls.

## Core Implementation Details in src/main.cpp

The repository's performance characteristics are defined in [`src/main.cpp`](https://github.com/cubao/pybind11-rdp/blob/main/src/main.cpp), which contains both the recursive and iterative algorithm variants.

### Recursive vs. Iterative Algorithms

The module exposes two algorithmic approaches:

- **`douglas_simplify`** (lines 50-84): A recursive depth-first implementation that processes segments recursively. This variant may be faster for small inputs due to lower overhead but risks stack overflow with deep recursion on large datasets.
- **`douglas_simplify_iter`** (lines 86-126): The queue-based iterative implementation using `std::queue<std::pair<int, int>>`. This version provides stable memory usage and avoids recursion limits, making it suitable for production use with large point clouds.

The Python wrapper in [`src/pybind11_rdp/__init__.py`](https://github.com/cubao/pybind11-rdp/blob/main/src/pybind11_rdp/__init__.py) (lines 66-95) selects between these implementations based on the `algo` parameter, defaulting to the iterative method for safety.

### The Distance Calculation Bottleneck

The original pure-Python rdp library computes point-to-line distance incorrectly (using point-to-line rather than point-to-segment) and performs this calculation in Python loops. In `pybind11-rdp`, the `LineSegment::distance2` method implements the mathematically correct point-to-segment distance using vector projection:

```cpp
double distance2(const Eigen::Vector3d& P) const {
    double t = (P - A).dot(AB) * inv_len2;
    t = std::max(0., std::min(1., t));  // Clamp to segment
    return (A + t * AB - P).squaredNorm();
}

```

This implementation uses Eigen's vectorized dot products and avoids branching in the hot loop, allowing the compiler to generate SIMD instructions.

## Python API Layer in src/pybind11_rdp/__init__.py

The high-level Python interface provides a clean API while delegating all computation to the C++ backend. The `rdp` function in [`src/pybind11_rdp/__init__.py`](https://github.com/cubao/pybind11-rdp/blob/main/src/pybind11_rdp/__init__.py) (lines 66-95) handles input validation, coordinate conversion, and result formatting:

```python
def rdp(coords, epsilon=0.0, algo="iter", return_mask=False):
    """
    Simplify a line using the Ramer-Douglas-Peucker algorithm.
    
    Parameters
    ----------
    coords : array_like
        Input coordinates (N, 2) or (N, 3)
    epsilon : float
        Simplification threshold
    algo : str
        'rec' for recursive, 'iter' for iterative
    return_mask : bool
        If True, return boolean mask instead of simplified coords
    """
    coords = np.asarray(coords, dtype=np.float64)
    if coords.ndim != 2 or coords.shape[1] not in (2, 3):
        raise ValueError("coords must be of shape (N, 2) or (N, 3)")
    
    # Delegate to C++ implementation

    if algo == "rec":
        mask = _core.rdp_mask_recursive(coords, epsilon)
    else:
        mask = _core.rdp_mask_iterative(coords, epsilon)
    
    if return_mask:
        return mask
    return coords[mask.astype(bool)]

```

This wrapper ensures that users interact with a standard NumPy-based API while the heavy lifting occurs in compiled code.

## Benchmarking the 8000x Speedup

The repository includes a benchmark script in [`test.py`](https://github.com/cubao/pybind11-rdp/blob/main/test.py) (lines 37-84) that quantifies the performance difference between the pure-Python `rdp` library and the C++ implementation. For a dataset of 2,000 3-D points, the results typically show:

```python
import time
import numpy as np
from pybind11_rdp import rdp as cpp_rdp
from rdp import rdp as py_rdp

# Generate test data: 2000 3D points

np.random.seed(42)
coords = np.random.rand(2000, 3)

# Benchmark pure-Python implementation

start = time.time()
py_result = py_rdp(coords, epsilon=0.1)
py_time = time.time() - start

# Benchmark C++ implementation

start = time.time()
cpp_result = cpp_rdp(coords, epsilon=0.1)
cpp_time = time.time() - start

print(f"Python rdp:  {py_time:.4f}s")
print(f"C++ rdp:     {cpp_time:.6f}s")
print(f"Speedup:     {py_time/cpp_time:.0f}x")

```

Typical output demonstrates a **speedup of approximately 8,000×**, with the Python implementation taking several seconds while the C++ version completes in microseconds. This dramatic difference stems from the elimination of Python interpreter overhead in the tight inner loops of the distance calculations.

## Summary

- **Compiled C++ Core**: The algorithm resides in [`src/main.cpp`](https://github.com/cubao/pybind11-rdp/blob/main/src/main.cpp) as native machine code, bypassing Python's interpreter loop entirely during the heavy computational phases.

- **Eigen Vectorization**: Linear algebra operations use the Eigen library with SIMD optimizations, delivering cache-friendly, vectorized distance calculations that Python loops cannot match.

- **Zero-Copy Memory Bridge**: `Eigen::Ref<const RowVectors>` adapters in the pybind11 bindings allow direct memory access to NumPy arrays without marshalling overhead.

- **Algorithmic Optimizations**: Pre-computed segment data in the `LineSegment` struct and the availability of both recursive and iterative (`std::queue`-based) variants ensure optimal performance across dataset sizes.

- **Correctness & Speed**: The implementation fixes the point-to-line distance bug present in the pure-Python version while achieving ~8,000× acceleration on 2,000-point datasets.

## Frequently Asked Questions

### Why is pybind11-rdp faster than the pure-Python rdp library?

The pure-Python rdp library executes the Ramer-Douglas-Peucker algorithm within Python's interpreter, incurring significant overhead for every distance calculation and loop iteration. In contrast, pybind11-rdp implements the entire algorithm in compiled C++ within [`src/main.cpp`](https://github.com/cubao/pybind11-rdp/blob/main/src/main.cpp), using Eigen for vectorized linear algebra and eliminating interpreter dispatch entirely during the hot path. This compiled approach, combined with zero-copy NumPy access via `Eigen::Ref`, explains the approximately 8,000× speedup observed in benchmarks.

### Does pybind11-rdp support 3D coordinates?

Yes, the library natively supports both 2D and 3D coordinate arrays. The C++ implementation in [`src/main.cpp`](https://github.com/cubao/pybind11-rdp/blob/main/src/main.cpp) uses `Eigen::Vector3d` for all point operations, allowing the `LineSegment` distance calculations to work uniformly across dimensions. When passing NumPy arrays to the `rdp()` or `rdp_mask()` functions, you can provide shapes `(N, 2)` or `(N, 3)`, and the wrapper in [`src/pybind11_rdp/__init__.py`](https://github.com/cubao/pybind11-rdp/blob/main/src/pybind11_rdp/__init__.py) automatically handles the coordinate dimensions without data copying.

### What is the difference between the recursive and iterative algorithms?

The repository provides two algorithmic variants in [`src/main.cpp`](https://github.com/cubao/pybind11-rdp/blob/main/src/main.cpp): `douglas_simplify` (recursive) and `douglas_simplify_iter` (iterative). The recursive version uses depth-first recursion to process line segments, which may offer slightly lower overhead for small datasets but risks stack overflow with deep recursion on large inputs. The iterative version uses `std::queue` to manage segments explicitly, providing stable memory usage and avoiding recursion limits while maintaining O(N log N) complexity. The Python API exposes both via the `algo="rec"` or `algo="iter"` parameter, defaulting to iterative for safety.

### How does the zero-copy NumPy integration work?

The zero-copy bridge relies on pybind11's `Eigen::Ref` type mapper, which allows C++ functions to accept NumPy arrays as `Eigen::Ref<const RowVectors>` without copying data. When you call `rdp(coords, epsilon=0.1)` from Python, the NumPy array's underlying memory buffer is passed directly to the C++ implementation in [`src/main.cpp`](https://github.com/cubao/pybind11-rdp/blob/main/src/main.cpp). The Eigen reference provides a matrix-like interface over the existing memory, allowing the `LineSegment` calculations and simplification loops to read coordinates at machine speed without Python API calls or memory allocation overhead during the algorithm execution.