Douglas-Peucker折线抽稀算法 Python 通过距离阈值控制和通过节点数量控制两种实现

折线抽稀

当折线点的数量过多时，往往需要对折线简化，且简化的折线能最大程度的保留原始折线的几何特征，这就是折线抽稀。主要的抽稀算法有Ramer-Douglas-Peucker、Visvalingam-Whyatt等，本文介绍Ramer-Douglas-Peucker算法。

Douglas-Peucker原理

具体原理可以参考Wiki或者其他博客，在此不过多介绍简单来说，Douglas-Peucker算法首先连接首尾两端，然后寻找折线上距离该线最远的点并计算距离，当距离大于距离阈值时，则保留改点，然后对首点和该点之间的折线部分和该点和尾点之间的折线部分重复上述过程；当距离小于距离阈值时，则用连线代替该部分折线。

Python 实现

调库实现

使用 https://rdp.readthedocs.io/en/latest/ 提供的库，代码如下

from rdp import rdp
import matplotlib.pyplot as plt
coords = [[1, 2], [2, 4], [3, 3], [4, 6],[7,9],[8,5],[9,2]]
simplify1 = rdp(coords,epsilon=1)
plt.plot([i[0] for i in coords],[i[1] for i in coords],color='r',marker='*',markersize=10)
plt.plot([i[0] for i in simplify1],[i[1] for i in simplify1],color='b',marker='*',markersize=10)
plt.show()

结果如下，红色为原始折线，蓝色为抽稀后，距离阈值为1

递归实现

该部分代码参考https://blog.csdn.net/rover2002/article/details/107026278
解释一下计算点到直线距离的部分，设目标点p(x, y)，线段（直线）两点p₁(x₁,y₁)，p₂(x₂,y₂)，则点p到p1 p2所在直线的距离为
$$d=\frac{|\overrightarrow{p_{1}p}\times \overrightarrow{p_1{p}_2}|}{|\overrightarrow{p_1{p}_2}|}$$
用坐标表示为
$$d=\frac{|(x-x_1)(y_2-y_1)-(x_2-x_1)(y-y_1)|}{\sqrt{(x_2-x_1{)}^2+(y_2-y_1{)}^2}}$$
代码：

import math

class Point(object):
    def __init__(self, id, x, y):
        self.id = id
        self.x = x
        self.y = y
    def __str__(self):
        return '{} {} {}'.format(self.id,self.x,self.y)

class DPCompress(object):
    def __init__(self, pointList, tolerance):
        self.Compressed = list()
        self.pointList = pointList
        self.tolerance = tolerance
        self.runDP(pointList, tolerance)

    def calc_height(self, point1, point2, point):
        """
        计算point到[point1, point2所在直线]的距离
        点到直线距离：
        A = point2.y - point1.y;
        B = point1.x - point2.x;
        C = point2.x * point1.y - point1.x * point2.y
        Dist = abs((A * point3.X + B * point3.Y + C) / sqrt(A * A + B * B))
        """
        tops = abs(point1.x * point.y + point2.x * point1.y + point.x * point2.y
                - point1.x * point2.y - point2.x * point.y - point.x * point1.y
                )

        bottom = math.sqrt(
            math.pow(point2.y - point1.y, 2) + math.pow(point2.x - point1.x, 2)
        )

        height = tops / bottom
        return height

    def DouglasPeucker(self, pointList, firsPoint, lastPoint, tolerance):
        """
        计算通过的内容
        DP算法
        :param pointList: 点列表
        :param firsPoint: 第一个点
        :param lastPoint: 最后一个点
        :param tolerance: 容差
        :return:
        """
        maxDistance = 0.0
        indexFarthest = firsPoint
        for i in range(firsPoint, lastPoint):
            distance = self.calc_height(pointList[firsPoint], pointList[lastPoint], pointList[i])
            if (distance > maxDistance):
                maxDistance = distance
                indexFarthest = i
    #    print('max_dis=', maxDistance)

        if maxDistance > tolerance and indexFarthest != firsPoint:
            self.Compressed.append(pointList[indexFarthest])
            self.DouglasPeucker(pointList, firsPoint, indexFarthest, tolerance)
            self.DouglasPeucker(pointList, indexFarthest, lastPoint, tolerance)

    def runDP(self, pointList, tolerance):
        """
        主要运行结果
        :param pointList: Point 列表
        :param tolerance: 值越小，压缩后剩余的越多
        :return:
        """
        if pointList == None or pointList.__len__() < 3:
            return pointList

        firspoint = 0
        lastPoint = len(pointList) - 1

        self.Compressed.append(pointList[firspoint])
        self.Compressed.append(pointList[lastPoint])

        self.DouglasPeucker(pointList, firspoint, lastPoint, tolerance)

    def getCompressed(self):
        self.Compressed.sort(key=lambda point: int(point.id))
        return self.Compressed

测试：

coords = [[1, 2], [2, 4], [3, 3], [4, 6],[7,9],[8,5],[9,2]]
pointList = []
for i,coord in enumerate(coords):
    pointList.append(Point(i,coord[0],coord[1]))
dp = DPCompress(pointList,1) # 距离阈值为1
simplify2 = []
for p in dp.getCompressed():
    simplify2.append([p.x,p.y])
plt.plot([i[0] for i in coords],[i[1] for i in coords],color='r',marker='*',markersize=10)
plt.plot([i[0] for i in simplify2],[i[1] for i in simplify2],color='b',marker='*',markersize=10)
plt.show()

输出：

控制结果点数的Douglas-Peucker算法

通常我们更需要控制抽稀的折线点数，例如从20个点抽稀到10个点，通过上述算法设置距离阈值很难满足要求，因此重新设计了可以输出指定数量结果的Douglas-Peucker算法。相当于将距离阈值设置为0，然后在计算的过程中，记录下每个点计算的距离，最后将距离从大到小排序，保留前n个点即可。
注意，由于没有找到严格的按点数抽稀的算法，因此本算法仅供参考，只能大致保存折线特征

import math

class Point(object):
    def __init__(self, id, x, y):
        self.id = id
        self.x = x
        self.y = y
        self.distance = float("inf")
    def __str__(self):
        return '{} {} {} {}'.format(self.id,self.x,self.y, self.distance)

class DPCompress(object):
    def __init__(self, pointList):
        self.Compressed = list()
        self.pointList = pointList
        self.runDP(pointList)

    def calc_height(self, point1, point2, point):
        """
        计算point到[point1, point2所在直线]的距离
        点到直线距离：
        A = point2.y - point1.y;
        B = point1.x - point2.x;
        C = point2.x * point1.y - point1.x * point2.y
        Dist = abs((A * point3.X + B * point3.Y + C) / sqrt(A * A + B * B))
        """
        tops = abs(point1.x * point.y + point2.x * point1.y + point.x * point2.y
                - point1.x * point2.y - point2.x * point.y - point.x * point1.y
                )

        bottom = math.sqrt(
            math.pow(point2.y - point1.y, 2) + math.pow(point2.x - point1.x, 2)
        )

        height = tops / bottom
        return height

    def DouglasPeucker(self, pointList, firsPoint, lastPoint):
        """
        计算通过的内容
        DP算法
        :param pointList: 点列表
        :param firsPoint: 第一个点
        :param lastPoint: 最后一个点
        :param tolerance: 容差
        :return:
        """
        maxDistance = 0.0
        indexFarthest = firsPoint
        for i in range(firsPoint+1, lastPoint):
            distance = self.calc_height(pointList[firsPoint], pointList[lastPoint], pointList[i])
            if distance < pointList[i].distance:
                pointList[i].distance = distance
            if (distance >= maxDistance):
                maxDistance = distance
                indexFarthest = i

        if indexFarthest != firsPoint:

            # pointList[indexFarthest].distance = maxDistance
            self.Compressed.append(pointList[indexFarthest])
            self.DouglasPeucker(pointList, firsPoint, indexFarthest)
            self.DouglasPeucker(pointList, indexFarthest, lastPoint)

    def runDP(self, pointList):
        """
        主要运行结果
        :param pointList: Point 列表
        :param tolerance: 值越小，压缩后剩余的越多
        :return:
        """
        if pointList == None or pointList.__len__() < 3:
            return pointList

        firspoint = 0
        lastPoint = len(pointList) - 1

        self.Compressed.append(pointList[firspoint])
        self.Compressed.append(pointList[lastPoint])

        self.DouglasPeucker(pointList, firspoint, lastPoint)

    def getCompressed(self,pointCount):
        
        self.Compressed.sort(key=lambda point: point.distance,reverse=True)
        result = self.Compressed[:pointCount]
        result.sort(key=lambda point: point.id)
        
        return result

结果

coords = [[1, 2], [2, 4], [3, 3], [4, 6],[7,9],[8,5],[9,2]]
pointList = []
for i,coord in enumerate(coords):
    pointList.append(Point(i,coord[0],coord[1]))
dp = DPCompress(pointList)
simplify2 = []
for p in dp.getCompressed(4):
    simplify2.append([p.x,p.y])
plt.plot([i[0] for i in coords],[i[1] for i in coords],color='r',marker='*',markersize=10)
plt.plot([i[0] for i in simplify2],[i[1] for i in simplify2],color='b',marker='*',markersize=10)
plt.show()