道格拉斯-普克算法 python实现

 道格拉斯-普克算法（Douglas–Peucker algorithm）

将曲线近似表示为一系列点，并减少点的数量的一种算法。它的优点是具有平移和旋转不变性，给定曲线与阈值后，抽样结果一定。可用于矢量地图数字水印算法

算法的基本思路是:对每一条曲线的首末点虚连一条直线,求所有点与直线的距离,并找出最大距离值dmax ,用dmax与限差D相比:若dmax＜D,这条曲线上的中间点全部舍去;若dmax ≥D,保留dmax 对应的做标点,并以该点为界,把曲线分为两部分,对这两部分重复使用该方法

 ​

 

以下为Python代码实现

输入：坐标数组以及阈值

返回：一个要素中的特征点以及对应的索引

import math
from decimal import Decimal


# 道格拉斯普克算法 输出一条直线上的特征点以及对应的索引
def point_to_line_Distance(point_a, point_b, point_c):
    """
    计算点 a 到点 b 和点 c 所在直线的垂直距离
    :param point_a: 要计算距离的点
    :param point_b: 直线的一端点
    :param point_c: 直线的另一端点
    :return: 垂直距离
    """
    if point_b.all() == point_c.all():  # 若直线两端点为同一点,则范围点与直线其中一端点的距离
        distance = Decimal(math.sqrt((point_b[0] - point_a[0]) ** 2 + (point_b[1] - point_a[1]) ** 2))
    elif point_b[0] == point_c[0]:  # 若直线为垂线
        distance = Decimal(abs(point_a[0] - point_b[0]))
    else:
        slope = Decimal(point_b[1] - point_c[1]) / Decimal(point_b[0] - point_c[0])
        intercept = Decimal(point_b[1]) - slope * Decimal(point_b[0])
        distance = abs(slope * Decimal(point_a[0]) - Decimal(point_a[1]) + intercept) / Decimal(
            math.sqrt(1 + slope ** 2))
    return distance


def douglas_peuker(point_list, threshold, lowerLimit=2, ceiling=None):
    """
    道格拉斯-普克算法用于线简化
    :param point_list: 点列表 形似[[0, 0], [1, 2], [2, 1], [3, 5], [4, 0], [0, 0]]
    :param threshold: 简化的距离阈值
    :param lowerLimit: 简化后点的最小数量
    :param ceiling: 简化后点的最大数量
    :return: 简化后的点列表及其索引
    """

    def diluting(points, current_disqualify_indices, threshold, qualify_points, disqualify_points, qualify_indices,
                 disqualify_indices):
        if len(points) < 3:
            qualify_points.extend(points[::-1])
            qualify_indices.extend(current_disqualify_indices[::-1])
        else:
            max_distance_index, max_distance = 0,0
            for index, point in enumerate(points):
                if index in [0, len(points) - 1]:
                    continue
                distance = point_to_line_Distance(point, points[0], points[-1])
                if distance > max_distance:
                    max_distance_index = index
                    max_distance = distance
            if max_distance < threshold:
                qualify_points.append(points[-1])
                qualify_points.append(points[0])
                qualify_indices.append(current_disqualify_indices[-1])
                qualify_indices.append(current_disqualify_indices[0])
            else:
                sequence_a = points[:max_distance_index]
                sequence_b = points[max_distance_index:]
                indices_a = current_disqualify_indices[:max_distance_index]
                indices_b = current_disqualify_indices[max_distance_index:]

                disqualify_points.append(sequence_a)
                disqualify_points.append(sequence_b)
                disqualify_indices.append(indices_a)
                disqualify_indices.append(indices_b)

        return qualify_points, disqualify_points, qualify_indices, disqualify_indices

    def get_qualify_list(points, threshold):
        qualify_points = []
        disqualify_points = [points]
        qualify_indices = []
        disqualify_indices = [list(range(len(points)))]

        while len(disqualify_indices) > 0:
            current_disqualify_points = disqualify_points.pop()
            current_disqualify_indices = disqualify_indices.pop()
            if not isinstance(current_disqualify_indices, list):
                current_disqualify_indices = [current_disqualify_indices]
            qualify_points, disqualify_points, qualify_indices, disqualify_indices = diluting(
                current_disqualify_points, current_disqualify_indices, threshold, qualify_points, disqualify_points,
                qualify_indices,
                disqualify_indices)

        return qualify_points, qualify_indices

    if len(point_list) < 5:
        return point_list, list(range(len(point_list)))

    result, indices = get_qualify_list(point_list, threshold)

    if len(result) < lowerLimit:
        while len(result) < lowerLimit:
            threshold *= 0.9
            result, indices = get_qualify_list(point_list, threshold)

    if ceiling and len(result) > ceiling:
        while len(result) > ceiling:
            threshold *= 1.1
            result, indices = get_qualify_list(point_list, threshold)

    if len(result) > len(point_list):
        return point_list, list(range(len(point_list)))

    return result, indices


if __name__ == '__main__':
    # 测试示例
    point_list = [[0, 0], [1, 2], [2, 1], [3, 5], [4, 0], [0, 0]]
    threshold = 4
    result, indices = douglas_peuker(point_list, threshold)
    print("简化后的点:", result)
    print("索引:", indices)