Douglas-Peucker 曲线简化算法详解

博客长期更新，本文最近一次更新时间为：2026年8月22日。

目录

• 一、算法概述
  • 1、应用背景
  • 2、算法简介
• 二、算法原理
  • 1、核心思想
  • 2、点到直线的距离
  • 3、递归分解过程
  • 4、时间复杂度分析
• 三、代码实现
  • 1、距离计算函数
  • 2、递归实现版本
  • 3、OpenCV 集成版
• 四、参数调优
  • 1、epsilon 的选择
  • 2、不同 epsilon 效果对比
• 五、工程应用
  • 1、焊缝轨迹简化
  • 2、机器人路径规划
  • 3、图像轮廓压缩
• 六、扩展与优化
  • 1、迭代版本（避免栈溢出）
  • 2、批量处理向量集合
• 七、参考文献

---

一、算法概述

1、应用背景

在计算机视觉、机器人路径规划、GIS 地理信息系统等领域，我们经常需要处理由大量离散点组成的曲线或折线。这些点往往存在大量冗余——相邻点之间的变化微乎其微，却占用了大量存储和计算资源。

典型场景：

• 视觉检测：边缘检测提取的轮廓可能有成百上千个点
• 机器人轨迹：激光扫描得到的焊缝轨迹点密度极高
• 地图数据：道路、海岸线的矢量数据量巨大

曲线简化（Curve Simplification）的目标是：在保持曲线整体形状的前提下，尽可能减少点的数量。

2、算法简介

Ramer–Douglas–Peucker 算法（简称 RDP 算法或 Douglas–Peucker 算法），由 Urs Ramer（1972）和 David Douglas & Thomas Peucker（1973）分别独立提出，是最经典的折线简化算法之一。

算法特点：

特性	说明
输入	有序点序列、距离阈值 ε（epsilon）
输出	简化后的点序列（原始点的子集）
核心思想	递归地用直线段近似曲线，移除偏差小于 ε 的中间点
时间复杂度	平均 O(n log n)，最坏 O(n²)
保真性	确保所有点到简化曲线的距离不超过 ε

---

二、算法原理

1、核心思想

Douglas-Peucker 算法采用分治策略，核心步骤如下：

1. 连接曲线的起点和终点，形成一条直线段
2. 计算曲线上所有中间点到该直线段的垂直距离
3. 找到距离最大的那个点，记其距离为 d_max
4. 如果 d_max < ε：所有中间点都可以舍弃，直接用首尾两点近似
5. 如果 d_max ≥ ε：以该点为界，将曲线分成两段，分别递归执行上述过程
6. 合并两段的结果（注意去重分界点）

原始曲线（8个点）：

  P2        P4
   ●        ●
    \     /  \
     ●   ●    ● P5
     P1  P3   |
              ● P6
               \
                ● P7

步骤1：连接 P1-P7，找最远点 P4
       d(P4, line_P1P7) > ε → 拆分

步骤2：左段 P1-P4，找最远点 P2
       d(P2, line_P1P4) < ε → 左段简化为 P1, P4

步骤3：右段 P4-P7，找最远点 P5
       d(P5, line_P4P7) > ε → 继续拆分...

最终结果：P1, P4, P5, P7（4个点，保留形状特征）

2、点到直线的距离

算法中最基础的计算是点到直线的垂直距离。

设直线由两点 A(x₁, y₁) 和 B(x₂, y₂) 确定，点 P(x₀, y₀) 到直线 AB 的距离公式：

        | (y₂ - y₁)·x₀ - (x₂ - x₁)·y₀ + x₂·y₁ - y₂·x₁ |
d = ──────────────────────────────────────────────────────
                 √((y₂ - y₁)² + (x₂ - x₁)²)

分子是向量 AB 和向量 AP 的叉积的绝对值，等于平行四边形的面积。
分母是 AB 的长度。
面积 / 底边长 = 高，即点到直线的垂直距离。

特殊情况：当 A 和 B 重合时（分母为 0），退化为点到点的欧氏距离。

3、递归分解过程

递归过程的伪代码：

function douglasPeucker(points, epsilon):
    if points.size() < 2:
        return points

    start = points[0]
    end = points[-1]

    // 找最远点
    maxDist = 0
    index = 0
    for i from 1 to points.size()-2:
        dist = perpendicularDistance(points[i], start, end)
        if dist > maxDist:
            maxDist = dist
            index = i

    if maxDist > epsilon:
        // 递归拆分
        left  = douglasPeucker(points[0..index], epsillon)
        right = douglasPeucker(points[index..-1], epsillon)
        // 合并（左段简化的最后一个点 = 右段简化的第一个点，去重）
        return left[0:.-1] + right
    else:
        // 直接用首尾两点近似
        return [start, end]

4、时间复杂度分析

最好情况：当曲线近似直线时，每次递归都能快速收敛，时间复杂度为 O(n log n)。这类似于快速排序的理想情况，每次划分都能将问题规模减半。

最坏情况：当曲线为高度曲折的锯齿形，且 ε 设置得非常小时，每次递归只能排除一个点，导致递归深度达到 n，时间复杂度退化为 O(n²)。这种情况在实际应用中较少见，但理论上存在。

平均情况：对于大多数真实世界的曲线（如自然轮廓、平滑轨迹），算法表现接近 O(n log n)。这是因为曲线通常具有局部平滑性，递归划分能有效减少计算量。

空间复杂度分析

递归版本的 Douglas-Peucker 算法在调用栈上的空间消耗与递归深度成正比：

• 最坏情况空间复杂度：O(n)，当递归深度达到 n 时（对应最坏时间复杂度情况）。
• 平均情况空间复杂度：O(log n)，对应平均时间复杂度 O(n log n) 的递归深度。

对于包含大量点的曲线（例如数万甚至数百万个点），递归版本可能因栈溢出而失败。为此，工程中常采用迭代版本（使用显式栈模拟递归过程）：

// 迭代版本：使用显式栈模拟递归，避免栈溢出
std::vector<Point> douglasPeuckerIterative(const std::vector<Point>& points, double epsilon) {
    if (points.size() <= 2) {
        return points;
    }

    std::vector<bool> keep(points.size(), false);
    keep.front() = keep.back() = true;

    // 显式栈，存储待处理的子段区间 [start, end]
    std::vector<std::pair<size_t, size_t>> stack;
    stack.emplace_back(0, points.size() - 1);

    while (!stack.empty()) {
        auto [start, end] = stack.back();
        stack.pop_back();

        if (end - start <= 1) {
            continue;
        }

        // 查找距离首尾直线最远的点
        const Point& pStart = points[start];
        const Point& pEnd = points[end];
        double maxDist = 0.0;
        size_t maxIndex = start;

        for (size_t i = start + 1; i < end; ++i) {
            double dist = perpendicularDistance(points[i], pStart, pEnd);
            if (dist > maxDist) {
                maxDist = dist;
                maxIndex = i;
            }
        }

        // 若最远距离超过阈值，则拆分并继续处理子段
        if (maxDist > epsilon) {
            keep[maxIndex] = true;
            stack.emplace_back(start, maxIndex);
            stack.emplace_back(maxIndex, end);
        }
        // 否则该子段中间点全部舍弃，仅保留首尾点
    }

    // 收集保留的点
    std::vector<Point> result;
    result.reserve(std::count(keep.begin(), keep.end(), true));
    for (size_t i = 0; i < points.size(); ++i) {
        if (keep[i]) {
            result.push_back(points[i]);
        }
    }
    return result;
}

递归 vs. 迭代空间占用对比

版本	空间复杂度	适用场景	备注
递归	O(n)（最坏） O(log n)（平均）	点集规模较小（通常 < 10⁴） 代码简洁，易于理解	递归调用栈由系统管理，深度过大时可能栈溢出
迭代	O(n)（显式栈） 实际占用通常小于递归栈	大规模点集（≥ 10⁴） 需要避免栈溢出的生产环境	显式栈的内存分配更可控，可预先分配固定大小数组

建议：对于未知规模或大规模点集，优先使用迭代版本；对于小规模数据或教学演示，递归版本更直观。- 不会栈溢出，支持超大量点集

• 通常比递归略快（减少函数调用开销）
• 内存使用可控

2、批量处理向量集合

对于多段曲线（如多条焊缝、多个轮廓），可以封装批量处理接口：

// 批量简化，返回每段的简化结果
template <typename PointT>
std::vector<std::vector<PointT>> simplifyAll(
    const std::vector<std::vector<PointT>>& polylines,
    float epsilon)
{
    std::vector<std::vector<PointT>> results;
    results.reserve(polylines.size());

    for (const auto& poly : polylines)
        results.push_back(douglasPeucker(poly, epsilon));

    return results;
}

---

三、代码实现

1、距离计算函数

首先实现点到直线距离的计算函数，这是算法的核心基础。

#include <cmath>
#include <vector>
#include <utility>

/**
 * 计算点到直线的垂直距离
 * 
 * @param point 待计算的点，格式为 std::pair<double, double> 或自定义点结构
 * @param start 直线起点
 * @param end 直线终点
 * @return 点到直线的垂直距离
 */
template<typename Point>
double perpendicularDistance(const Point& point, const Point& start, const Point& end) {
    // 如果起点和终点重合，退化为点到点的欧氏距离
    if (start == end) {
        double dx = point.first - start.first;
        double dy = point.second - start.second;
        return std::sqrt(dx * dx + dy * dy);
    }
    
    // 计算分子：向量叉积的绝对值
    double numerator = std::abs(
        (end.second - start.second) * point.first -
        (end.first - start.first) * point.second +
        end.first * start.second -
        end.second * start.first
    );
    
    // 计算分母：直线段的长度
    double denominator = std::sqrt(
        std::pow(end.second - start.second, 2) +
        std::pow(end.first - start.first, 2)
    );
    
    return numerator / denominator;
}

// 使用示例：
struct Point2D {
    double x, y;
    bool operator==(const Point2D& other) const {
        return x == other.x && y == other.y;
    }
};

// 为 Point2D 结构特化模板
double perpendicularDistance(const Point2D& point, const Point2D& start, const Point2D& end) {
    if (start == end) {
        double dx = point.x - start.x;
        double dy = point.y - start.y;
        return std::sqrt(dx * dx + dy * dy);
    }
    
    double numerator = std::abs(
        (end.y - start.y) * point.x -
        (end.x - start.x) * point.y +
        end.x * start.y -
        end.y * start.x
    );
    
    double denominator = std::sqrt(
        std::pow(end.y - start.y, 2) +
        std::pow(end.x - start.x, 2)
    );
    
    return numerator / denominator;
}

2、递归实现版本

下面是一个完整的递归实现，包含数据生成、算法调用和结果输出。

#include <iostream>
#include <vector>
#include <cmath>
#include <algorithm>
#include <random>
#include <chrono>

// 使用 std::pair 表示二维点
using Point = std::pair<double, double>;

// 距离计算函数（使用 std::pair）
double perpendicularDistance(const Point& point, const Point& start, const Point& end) {
    if (start == end) {
        double dx = point.first - start.first;
        double dy = point.second - start.second;
        return std::sqrt(dx * dx + dy * dy);
    }
    
    double numerator = std::abs(
        (end.second - start.second) * point.first -
        (end.first - start.first) * point.second +
        end.first * start.second -
        end.second * start.first
    );
    
    double denominator = std::sqrt(
        std::pow(end.second - start.second, 2) +
        std::pow(end.first - start.first, 2)
    );
    
    return numerator / denominator;
}

/**
 * Douglas-Peucker 算法的递归实现
 * 
 * @param points 有序点序列
 * @param epsilon 距离阈值，控制简化程度
 * @return 简化后的点序列（原始点的子集）
 */
std::vector<Point> douglasPeuckerRecursive(const std::vector<Point>& points, double epsilon) {
    // 基本情况：点数少于3个，直接返回
    if (points.size() <= 2) {
        return points;
    }
    
    // 找到距离首尾直线最远的点
    const Point& start = points.front();
    const Point& end = points.back();
    double maxDist = 0.0;
    size_t maxIndex = 0;
    
    for (size_t i = 1; i < points.size() - 1; ++i) {
        double dist = perpendicularDistance(points[i], start, end);
        if (dist > maxDist) {
            maxDist = dist;
            maxIndex = i;
        }
    }
    
    // 如果最大距离小于阈值，直接返回首尾两点
    if (maxDist <= epsilon) {
        return {start, end};
    }
    
    // 递归处理左右两段
    std::vector<Point> leftSegment(points.begin(), points.begin() + maxIndex + 1);
    std::vector<Point> rightSegment(points.begin() + maxIndex, points.end());
    
    std::vector<Point> leftResult = douglasPeuckerRecursive(leftSegment, epsilon);
    std::vector<Point> rightResult = douglasPeuckerRecursive(rightSegment, epsilon);
    
    // 合并结果，注意去重（左段的最后一个点 = 右段的第一个点）
    std::vector<Point> result;
    result.reserve(leftResult.size() + rightResult.size() - 1);
    result.insert(result.end(), leftResult.begin(), leftResult.end() - 1);
    result.insert(result.end(), rightResult.begin(), rightResult.end());
    
    return result;
}

/**
 * 生成测试曲线：正弦波叠加噪声
 */
std::vector<Point> generateTestCurve(int numPoints = 100) {
    std::vector<Point> points;
    points.reserve(numPoints);
    
    std::random_device rd;
    std::mt19937 gen(rd());
    std::normal_distribution<> noise(0.0, 0.1);
    
    for (int i = 0; i < numPoints; ++i) {
        double x = 10.0 * i / (numPoints - 1);
        double y = std::sin(x) + noise(gen);
        points.emplace_back(x, y);
    }
    
    return points;
}

/**
 * 输出结果到控制台和文件（CSV格式）
 */
void outputResults(const std::vector<Point>& original, 
                   const std::vector<Point>& simplified, 
                   double epsilon,
                   const std::string& filename = "results.csv") {
    
    std::cout << "\n=== 算法执行结果 ===\n";
    std::cout << "原始点数: " << original.size() << "\n";
    std::cout << "简化后点数: " << simplified.size() << "\n";
    std::cout << "压缩率: " << 100.0 * (1.0 - static_cast<double>(simplified.size()) / original.size()) << "%\n";
    std::cout << "距离阈值 ε: " << epsilon << "\n";
    
    // 计算最大误差
    double maxError = 0.0;
    for (const auto& point : original) {
        double minDist = std::numeric_limits<double>::max();
        for (size_t i = 0; i < simplified.size() - 1; ++i) {
            double dist = perpendicularDistance(point, simplified[i], simplified[i + 1]);
            if (dist < minDist) {
                minDist = dist;
            }
        }
        if (minDist > maxError) {
            maxError = minDist;
        }
    }
    
    std::cout << "最大误差: " << maxError << "\n";
    std::cout << "质量验证: " << (maxError <= epsilon ? "通过" : "失败") << "\n";
    
    // 输出到CSV文件（可用于外部可视化）
    std::ofstream outFile(filename);
    if (outFile.is_open()) {
        outFile << "type,x,y\n";
        outFile << "# Original points\n";
        for (const auto& p : original) {
            outFile << "original," << p.first << "," << p.second << "\n";
        }
        outFile << "# Simplified points\n";
        for (const auto& p : simplified) {
            outFile << "simplified," << p.first << "," << p.second << "\n";
        }
        outFile.close();
        std::cout << "结果已保存到: " << filename << " (可用Excel/Python matplotlib等工具可视化)\n";
    }
}

/**
 * 主程序：完整的测试流程
 */
int main() {
    std::cout << "Douglas-Peucker 曲线简化算法 (C++ 递归实现)\n";
    std::cout << "===========================================\n";
    
    // 1. 生成测试数据
    std::cout << "\n步骤1: 生成测试曲线...\n";
    auto testPoints = generateTestCurve(100);
    std::cout << "生成 " << testPoints.size() << " 个测试点\n";
    
    // 2. 设置阈值并执行算法
    double epsilon = 0.2;
    std::cout << "\n步骤2: 执行 Douglas-Peucker 算法 (ε=" << epsilon << ")...\n";
    
    auto startTime = std::chrono::high_resolution_clock::now();
    auto simplified = douglasPeuckerRecursive(testPoints, epsilon);
    auto endTime = std::chrono::high_resolution_clock::now();
    
    auto duration = std::chrono::duration_cast<std::chrono::microseconds>(endTime - startTime);
    std::cout << "算法执行时间: " << duration.count() << " μs\n";
    
    // 3. 输出结果
    std::cout << "\n步骤3: 输出结果...\n";
    outputResults(testPoints, simplified, epsilon);
    
    // 4. 显示简化后的点
    std::cout << "\n步骤4: 简化后的点序列:\n";
    std::cout << "索引\tX\t\tY\n";
    std::cout << "----\t--------\t--------\n";
    for (size_t i = 0; i < simplified.size(); ++i) {
        std::cout << i << "\t" << simplified[i].first << "\t" << simplified[i].second << "\n";
    }
    
    // 5. 性能测试：多次运行取平均
    std::cout << "\n步骤5: 性能测试 (运行100次取平均)...\n";
    const int numRuns = 100;
    long long totalTime = 0;
    
    for (int i = 0; i < numRuns; ++i) {
        auto start = std::chrono::high_resolution_clock::now();
        auto _ = douglasPeuckerRecursive(testPoints, epsilon);
        auto end = std::chrono::high_resolution_clock::now();
        totalTime += std::chrono::duration_cast<std::chrono::microseconds>(end - start).count();
    }
    
    double avgTime = static_cast<double>(totalTime) / numRuns;
    std::cout << "平均执行时间: " << avgTime << " μs\n";
    std::cout << "每秒可处理点数: " << (testPoints.size() / (avgTime / 1e6)) << " points/s\n";
    
    std::cout << "\n=== 测试完成 ===\n";
    
    std::cout << R"(
#include <fstream>
#include <iostream>
#include <sstream>
#include <vector>

// 读取 CSV 并按类型分组，使用 matplotlib 风格输出（需配合 gnuplot 或自绘）
int main() {
    std::ifstream in("results.csv");
    std::string line;
    std::vector<std::pair<double, double>> orig, simp;

    while (std::getline(in, line)) {
        if (line.empty() || line[0] == '#') continue;
        std::stringstream ss(line);
        std::string type, xs, ys;
        std::getline(ss, type, ',');
        std::getline(ss, xs, ',');
        std::getline(ss, ys, ',');
        double x = std::stod(xs), y = std::stod(ys);
        if (type == "original") orig.emplace_back(x, y);
        else if (type == "simplified") simp.emplace_back(x, y);
    }

    // 输出 gnuplot 绘图脚本
    std::cout << "plot '-' with lines title '原始曲线', '-' with points title '简化点'\\n";
    for (auto& p : orig) std::cout << p.first << " " << p.second << "\\n";
    std::cout << "e\\n";
    for (auto& p : simp) std::cout << p.first << " " << p.second << "\\n";
    std::cout << "e\\n";
    return 0;
}
    )" << std::endl;
    
    return 0;
}

4、使用建议

1. 阈值选择：ε 值越大，简化程度越高，但可能丢失细节。通常通过试验确定，或根据应用需求设定（如地图显示精度要求）。

2. 性能考虑：

   • 对于小规模数据（< 1000点），递归版本足够
   • 对于中等规模数据（1000-10000点），建议使用迭代版本
   • 对于大规模数据（> 10000点），优先使用 OpenCV 等优化库

3. 内存管理：

   • 递归版本可能栈溢出，使用迭代版本更安全
   • 对于批量处理，考虑分块处理或使用生成器

4. 精度控制：

   • 算法保证所有点到简化曲线的距离 ≤ ε
   • 可通过后处理验证简化质量

5. 实际应用：

   • 图像处理：轮廓简化、边缘检测后处理
   • GIS：地图数据压缩、路径规划
   • 机器人：轨迹平滑、运动规划

七、参考文献

[1] Ramer, U. An iterative procedure for the polygonal approximation of plane curves. Computer Graphics and Image Processing, 1(3):244-256, 1972.

[2] Douglas, D. H., Peucker, T. K. Algorithms for the reduction of the number of points required to represent a digitized line or its caricature. Cartographica: The International Journal for Geographic Information and Geovisualization, 10(2):112-122, 1973.

[3] OpenCV: approxPolyDP

[4] Wikipedia: Ramer-Douglas-Peucker algorithm