我编写了这段代码,以删除绘制时会形成直线的点。

例如,在一条直线上有许多点:


import matplotlib.pyplot as plt
import numpy as np

pattern = np.array(
    [ 9, 10, 11, 10,  9,  8,  6,  5,  4,  3,  4,  5,  7,  6,  5,  4,  2,
      1,  2,  3,  4,  5,  7,  8,  9, 10,  8,  7,  8,  9, 10, 11, 10, 11,
     12, 13, 12, 11, 10,  8,  7,  6,  5,  4,  6,  5,  3,  2,  3,  2,  0,
      1,  2,  3,  4,  5,  6,  8,  9, 10, 11, 12, 13, 14, 15, 16, 17, 19,
     21, 19, 18, 16, 15, 14, 12, 11, 12, 11, 10, 11, 12, 13, 14, 15, 13,
     12, 11, 10,  9,  8,  7,  6,  5,  6,  5,  4,  6])

plt.plot(pattern)
plt.plot(pattern, 'bo')
plt.show()





使用degrees函数,我可以删除在两个相邻点之间形成180°夹角的线。

from numpy.lib.stride_tricks import as_strided

def moving_slice(a, k):
    a = a.ravel()
    return as_strided(a, (a.size - k + 1, k), 2 * a.strides)

def degrees(x, y):
    x_window, y_window = moving_slice(x, 3), moving_slice(y, 3)
    x_window = np.concatenate((x_window, x_window[:, None, 0]), axis=1)
    y_window = np.concatenate((y_window, y_window[:, None, 0]), axis=1)

    triangle_sides = np.sqrt(np.square(np.diff(x_window)) + np.square(np.diff(y_window)))
    squared_sides = np.square(triangle_sides)

    cos_num = np.sum(squared_sides[:, :2], axis=1) - squared_sides[:, 2]
    cos_den = 2 * np.prod(triangle_sides[:, :2], axis=1)
    return np.degrees(np.arccos(cos_num / cos_den))


这是我要实现的目标:

# My review function
angles = degrees(np.arange(len(pattern)), pattern)

pattern = pattern[1:-1] # <- No angles between the end points
less_than_180 = angles < 180
plt.plot(pattern[less_than_180])
plt.plot(pattern[less_than_180], 'bo')
plt.show()




您可以看到很好。遗漏了一些要删除的点(由于十进制精度),否则该功能可以满足我的要求。我想知道是否有一种更简单的方法可以在不计算所有角度的情况下达到期望的结果。你觉得呢?

评论

考虑如果两个相邻点之间的dx相同,则它们具有相同的斜率iff dy也相同。如果您的dx总是固定的(如您的测试示例中所示),则此事实可能有助于简化操作。

@DanBryant好观察。是的,在我的用例中,x始终是连续的整数。这样也可以。

#1 楼



看一下帖子中的第二个情节,发现出了点问题。一行中有四个点:



这应该是不可能的,因为已删除了“绘制时会形成直线的点”。这怎么发生的?看一下\ $ x \ $轴就可以找到原因。第一个图的\ $ x \ $坐标从0到96,但是第二个图的\ $ x \ $坐标从0到43,因为省略点的\ $ x \ $值也已经

要获得正确的绘图,您必须记住原始\ $ x \ $坐标并在遮罩\ $ y \ $值的同时遮罩它们,如下所示:

x = np.arange(len(pat))[less_than_180]
y = pat[less_than_180]
plt.plot(x, y)
plt.plot(x, y, 'bo')


请参见下面的更正图。在NumPy中,读者可能会想象您的函数执行了类似的操作。像degrees这样的名称会更好(请参阅Wikipedia)。

该行:

pattern = pattern[1:-1] # <- No angles between the end points


丢弃数据的端点。但是确定要在曲线图中保留端点吗?


给出三个点\ $ A,B,C \ $,代码找到顶点角\ $∠ABC\ $和将其与180°进行比较。另一种方法是在端点之间进行线性插值,找到\ $ B'=½(A + C)\ $,然后查看\ $ B'\ $是否足够接近\ $ B \ $。像这样:

# Boolean mask for the points we are going to plot.
mask = np.ones_like(pattern, dtype=bool)

if len(pattern) > 2:
    # Linearly interpolate between points that are two apart.
    interpolants = (pattern[:-2] + pattern[2:]) / 2

    # Don't plot points that are equal to their interpolant.
    mask[1:-1] = interpolants != pattern[1:-1]

x = np.arange(len(pattern))[mask]
y = pattern[mask]
plt.plot(x, y)
plt.plot(x, y, 'bo')
plt.show()


结果图如下:



如果数据是浮点则在将点与内插值进行比较时,您将需要一些公差。例如,您可以使用math.degrees

# Don't plot points that are close their interpolant.
mask[1:-1] = ~np.isclose(interpolants, pattern[1:-1])




评论


\ $ \ begingroup \ $
这是一个了不起的答案!谢谢你的时间。我不了解术语“顶角”或“线性插值”,我需要做一些阅读工作!是的,我确实想要终点:)
\ $ \ endgroup \ $
–詹姆斯·辛纳(James Schinner)
18-2-8在12:56



\ $ \ begingroup \ $
作为线性插值的一种替代方法,也可以使用以下事实:如果中间点的二阶导数为零,则3个点在同一行上:mask [1:-1] = np.diff(pattern,n = 2)!= 0
\ $ \ endgroup \ $
–乔治
18年2月8日在20:26

\ $ \ begingroup \ $
嗯,这也是一个有趣的想法。哇,我真的错了。
\ $ \ endgroup \ $
–詹姆斯·辛纳(James Schinner)
18年2月8日在23:09

#2 楼

您可以使用Ramer Douglas Peuker算法。 RDP绘制曲线并消除接近直线的点。它是基于距离的。从两个端点开始,它形成一条线,并选择距离其最远的点。该点用于形成线段并进行迭代,直到所有其余点都靠近该线为止。

Python具有可使用的RDP包。您需要做的就是定义线条的直线度。

结果:



评论


\ $ \ begingroup \ $
好,那很好。
\ $ \ endgroup \ $
–詹姆斯·辛纳(James Schinner)
18年2月9日在2:12

#3 楼

这实际上不是代码审查,而是解决问题的另一种方法。

我的博士论文是关于折线对噪声数据序列的双标准拟合。两个标准是段数和均方估计误差。权重参数允许在两者之间进行权衡。以下是使用您的数据和权衡参数的不同值得出的一些结果:

0.1



1.0



4.0



20.0



重要减少段数时,选择结果点以最大程度地减少误差,并且不限于等于任何输入点(因为假定它们被噪声破坏了)。

我没有可以与您共享的python版本代码,但是该算法(这是一种动态编程方法)已在https://repository.upenn.edu/edissertations/1158/在线描述,我建议查看显示错误性能和计算复杂度的图表,而不是尝试重现代码。如果这是您想要的,可以说服我在python绑定上进行协作。

评论


\ $ \ begingroup \ $
哇,谢谢你的分享!我正在读你的论文。
\ $ \ endgroup \ $
–詹姆斯·辛纳(James Schinner)
18年2月9日在21:19