修改 numpy 中的 Ramer–Douglas–Peucker (RDP) 算法以返回掩码而不是值答案

【问题标题】：Modify the Ramer–Douglas–Peucker (RDP) algorithm in numpy to return a mask instead of the values修改 numpy 中的 Ramer–Douglas–Peucker (RDP) 算法以返回掩码而不是值
【发布时间】：2015-03-08 00:30:48
【问题描述】：

我正在尝试修改numpy implementation of the Ramer–Douglas–Peucker (RDP) algorithm，以返回过滤值的掩码，而不是带有值的过滤数组。

问题是如何递归返回正确的掩码。

【问题讨论】：

RDP 算法需要对单个 2D 点进行大量计算。例如，pldist 对三个 (2,1) 形的 numpy 数组进行操作。当在小数组上进行大量计算时，与纯 Python 代码相比，NumPy 相当慢。例如，Sebleier's non-NumPy implementation 的速度提高了约 42 倍。
这是有道理的。也许可以修改pldist 函数以对多个点进行操作，而不是循环遍历它们。

标签： python algorithm numpy

【解决方案1】：

这有助于正确理解算法。如果每个递归调用中第一个和最后一个段之间的所有点都在 epsilon 内，则这些点应标记为 False，而 start 和 end 应为 True。这是我的解决方案：

import numpy as np

def pldist(x0, x1, x2):
    return np.divide(np.linalg.norm(np.linalg.det([x2 - x1, x1 - x0])),
                     np.linalg.norm(x2 - x1))

def rdp_index(M, epsilon=0, dist=pldist):

    dmax = 0.0
    index = -1

    for i in xrange(1, M.shape[0]):
        d = dist(M[i], M[0], M[-1])

        if d > dmax:
            index = i
            dmax = d

    if dmax > epsilon:
        r1 = rdp_index(M[:index + 1], epsilon, dist)
        r2 = rdp_index(M[index:], epsilon, dist)

        return np.vstack((r1[:-1], r2))
    else:
        part_mask = np.empty_like(M, dtype = bool)
        part_mask.fill(False)
        part_mask[0] = True
        part_mask[-1] = True
        return part_mask

请注意，此实现使用递归，并且此解决方案在大型数据集https://stackoverflow.com/a/2401520/2082968 上存在问题。 IE。超出了最大递归调用量。下面是一个使用堆栈和 while 循环而不是递归调用的解决方案。此外，距离的计算效率更高。

def dsquared_line_points(P1, P2, points):
    '''
    Calculate only squared distance, only needed for comparison
    http://en.wikipedia.org/wiki/Distance_from_a_point_to_a_line
    '''
    xdiff = P2[0] - P1[0]
    ydiff = P2[1] - P1[1]
    nom  = (
        ydiff*points[:,0] - \
        xdiff*points[:,1] + \
        P2[0]*P1[1] - \
        P2[1]*P1[0]
    )**2
    denom = ydiff**2 + xdiff**2
    return np.divide(nom, denom)

def rdp_numpy(M, epsilon = 0):

    # initiate mask array
    # same amount of points
    mask = np.empty(M.shape[0], dtype = bool)

    # Assume all points are valid and falsify those which are found
    mask.fill(True)

    # The stack to select start and end index
    stack = [(0 , M.shape[0]-1)]

    while (len(stack) > 0):
        # Pop the last item
        (start, end) = stack.pop()

        # nothing to calculate if no points in between
        if end - start <= 1:
            continue

        # Calculate distance to points
        P1 = M[start]
        P2 = M[end]
        points = M[start + 1:end]
        dsq = dsquared_line_points(P1, P2, points)

        mask_eps = dsq > epsilon**2

        if mask_eps.any():
            # max point outside eps
            # Include index that was sliced out
            # Also include the start index to get absolute index
            # And not relative 
            mid = np.argmax(dsq) + 1 + start
            stack.append((start, mid))
            stack.append((mid, end))

        else:
            # Points in between are redundant
            mask[start + 1:end] = False

    return mask

【讨论】：