在轨迹计算转向点/枢转点（路径）(calculate turning points / pivot

2019-08-16 16:16发布

站内文章 / 后端开发

387 0

爷、活的狠高调

女 | 书童

私信

我试图想出一个算法，将决定X / Y坐标的轨迹的转折点。下图说明了我的意思是：绿色表示的出发点和红色轨迹的终点（整个轨迹由〜1500分）：

在下图中，我加手工的可能（全球）转折点，一个算法可以返回：

显然，真正的转折点总是有争议的，而且将取决于角度有一个指定的点之间。此外的一个转折点可以在全球范围内（我试图与黑色圆圈做）来定义，但也可以在高分辨率局部范围定义。我感兴趣的全局（整体）的方向变化，但我很乐意看到的，一个会用梳理出全球VS本地解决方案的不同方法的讨论。

我试过到目前为止：

计算随后的点之间的距离
计算后续点之间角度
看距离/角度后续点之间如何变化

不幸的是这并没有给我任何稳定的结果。我大概有太多计算沿多个点的曲率，但是这只是一个想法。我真的很感激任何算法/想法，可以帮助我在这里。代码可以是任何编程语言，MATLAB或Python是优选的。

这里编辑的原始数据（如果有人想对玩它）：

垫档
文本文件（x坐标第一，Y在第二线坐标）

Answer 1:

您可以使用拉默-道格拉斯-普克（RDP）算法，简化了路径。然后，你可以计算沿简化路径的每个段的方向变化。对应于方向变化最大的点可以被称为转折点：

在RDP算法的Python实现，可以发现在github上。

import matplotlib.pyplot as plt
import numpy as np
import os
import rdp

def angle(dir):
    """
    Returns the angles between vectors.

    Parameters:
    dir is a 2D-array of shape (N,M) representing N vectors in M-dimensional space.

    The return value is a 1D-array of values of shape (N-1,), with each value
    between 0 and pi.

    0 implies the vectors point in the same direction
    pi/2 implies the vectors are orthogonal
    pi implies the vectors point in opposite directions
    """
    dir2 = dir[1:]
    dir1 = dir[:-1]
    return np.arccos((dir1*dir2).sum(axis=1)/(
        np.sqrt((dir1**2).sum(axis=1)*(dir2**2).sum(axis=1))))

tolerance = 70
min_angle = np.pi*0.22
filename = os.path.expanduser('~/tmp/bla.data')
points = np.genfromtxt(filename).T
print(len(points))
x, y = points.T

# Use the Ramer-Douglas-Peucker algorithm to simplify the path
# http://en.wikipedia.org/wiki/Ramer-Douglas-Peucker_algorithm
# Python implementation: https://github.com/sebleier/RDP/
simplified = np.array(rdp.rdp(points.tolist(), tolerance))

print(len(simplified))
sx, sy = simplified.T

# compute the direction vectors on the simplified curve
directions = np.diff(simplified, axis=0)
theta = angle(directions)
# Select the index of the points with the greatest theta
# Large theta is associated with greatest change in direction.
idx = np.where(theta>min_angle)[0]+1

fig = plt.figure()
ax =fig.add_subplot(111)

ax.plot(x, y, 'b-', label='original path')
ax.plot(sx, sy, 'g--', label='simplified path')
ax.plot(sx[idx], sy[idx], 'ro', markersize = 10, label='turning points')
ax.invert_yaxis()
plt.legend(loc='best')
plt.show()

两个参数被用于以上：

的RDP算法需要一个参数，该tolerance ，它代表了简化路径可以从原始路径偏离的最大距离。在较大的tolerance ，简化路径的较粗略的。
另一个参数是min_angle定义什么被认为是一个转折点。（我正在转折点是原始路径中，进入和离开矢量简化路径上之间，其角度大于上的任何点min_angle ）。

Answer 2:

我将在下文中给出numpy的/ SciPy的代码，因为我几乎没有任何经验的Matlab的。

如果你的曲线是不够顺畅，你能确定你的转折点那些最高的曲率。取点索引号作为曲线参数，和一个中央差异方案，可以计算用下面的代码曲率

import numpy as np
import matplotlib.pyplot as plt
import scipy.ndimage

def first_derivative(x) :
    return x[2:] - x[0:-2]

def second_derivative(x) :
    return x[2:] - 2 * x[1:-1] + x[:-2]

def curvature(x, y) :
    x_1 = first_derivative(x)
    x_2 = second_derivative(x)
    y_1 = first_derivative(y)
    y_2 = second_derivative(y)
    return np.abs(x_1 * y_2 - y_1 * x_2) / np.sqrt((x_1**2 + y_1**2)**3)

你可能会想先用手整理的曲线出来，然后计算曲率，然后确定最高曲率点。下面的功能做到了这一点：

def plot_turning_points(x, y, turning_points=10, smoothing_radius=3,
                        cluster_radius=10) :
    if smoothing_radius :
        weights = np.ones(2 * smoothing_radius + 1)
        new_x = scipy.ndimage.convolve1d(x, weights, mode='constant', cval=0.0)
        new_x = new_x[smoothing_radius:-smoothing_radius] / np.sum(weights)
        new_y = scipy.ndimage.convolve1d(y, weights, mode='constant', cval=0.0)
        new_y = new_y[smoothing_radius:-smoothing_radius] / np.sum(weights)
    else :
        new_x, new_y = x, y
    k = curvature(new_x, new_y)
    turn_point_idx = np.argsort(k)[::-1]
    t_points = []
    while len(t_points) < turning_points and len(turn_point_idx) > 0:
        t_points += [turn_point_idx[0]]
        idx = np.abs(turn_point_idx - turn_point_idx[0]) > cluster_radius
        turn_point_idx = turn_point_idx[idx]
    t_points = np.array(t_points)
    t_points += smoothing_radius + 1
    plt.plot(x,y, 'k-')
    plt.plot(new_x, new_y, 'r-')
    plt.plot(x[t_points], y[t_points], 'o')
    plt.show()

有的解释是为了：

turning_points是要识别点数
smoothing_radius是平滑卷积的半径计算的曲率之前必须施加到数据
cluster_radius是选自作为转折点，其中没有其他点应被视为候选高曲率的点的距离。

您可能不得不玩的参数了一点，但我得到了这样的事情：

>>> x, y = np.genfromtxt('bla.data')
>>> plot_turning_points(x, y, turning_points=20, smoothing_radius=15,
...                     cluster_radius=75)

也许不够好完全自动化的检测，但它是非常接近你想要的。

Answer 3:

一个非常有趣的问题。这里是我的解决方案，它允许可变分辨率。虽然，微调也未必是简单的，因为它的主要目的是缩小

每k个点，计算所述凸包并将其存储为一组。通过至多k个点去，并删除不在凸包，在这样的点不会失去原来的顺序的方式任意点。

这里的目的是凸包将作为一个过滤器，清除所有的“不重要点”只留下了极值点。当然，如果k值太高，你会用什么，而不是你真正想要的东西太接近实际的凸包，结束了。

这应该有一个小k，至少4开始，然后增加它，直到你得到你所追求的。你也应该可能只针对每3个点的中间点的角度小于一定金额，d。这将确保所有的匝的至少d度（在下面的代码未实现）。然而，这也许应该逐步进行，以避免信息的损失，同样增加了K值。另一种可能的改进是实际上被删除的点重新运行，并且，仅删除在不同时凸包点，尽管这需要的至少8个较高的最小k值。

下面的代码似乎工作得相当好，但仍然可以使用效率和噪声去除的改进。这也是在确定何时应该停止，而不雅，因此代码真的只能（因为它代表）从大约K = 4到k = 14。

def convex_filter(points,k):
    new_points = []
    for pts in (points[i:i + k] for i in xrange(0, len(points), k)):
        hull = set(convex_hull(pts))
        for point in pts:
            if point in hull:
                new_points.append(point)
    return new_points

# How the points are obtained is a minor point, but they need to be in the right order.
x_coords = [float(x) for x in x.split()]
y_coords = [float(y) for y in y.split()]
points = zip(x_coords,y_coords)

k = 10

prev_length = 0
new_points = points

# Filter using the convex hull until no more points are removed
while len(new_points) != prev_length:
    prev_length = len(new_points)
    new_points = convex_filter(new_points,k)

下面是其中k = 14上面的代码的屏幕截图。 61个红点是依然存在的过滤后的人。