вычислять точки поворота / точки поворота в траектории (пути)

Я пытаюсь придумать алгоритм, который определит поворотные точки в траектории координат x / y. На следующих рисунках показано, что я имею в виду: зеленый указывает начальную точку и красную конечную точку траектории (вся траектория состоит из ~ 1500 точек): траектория

На следующем рисунке я добавил вручную возможные (глобальные) точки поворота, которые мог бы вернуть алгоритм:

траектория с возможными точками поворота

Очевидно, что истинный поворотный момент всегда дискуссионный и будет зависеть от угла, который указывается, который должен лежать между точками. Кроме того, поворотный момент можно определить в глобальном масштабе (что я пытался сделать с черными кругами), но также можно определить в локальном масштабе высокого разрешения. Меня интересуют глобальные (общие) изменения направления, но мне бы хотелось обсудить обсуждение различных подходов, которые можно было бы использовать, чтобы разделить глобальные и локальные решения.

То, что я пробовал до сих пор:

  • рассчитать расстояние между последующими точками
  • вычислить угол между последующими точками
  • посмотрите, как расстояние / угол изменяется между последующими точками

К сожалению, это не дает мне никаких надежных результатов. Я, вероятно, слишком вычислил кривизну вдоль нескольких точек, но это всего лишь идея. Я бы очень признателен за любые алгоритмы / идеи, которые могут помочь мне здесь. Код может быть на любом языке программирования, предпочтительнее использовать matlab или python.

EDIT – это исходные данные (в случае, если кто-то хочет сыграть с ним):

  • мат-файл
  • текстовый файл (сначала координата x, координата y во второй строке)

Вы можете использовать алгоритм Ramer-Douglas-Peucker (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() 

введите описание изображения здесь

Два параметра были использованы выше:

  1. Алгоритм RDP принимает один параметр, tolerance , который представляет максимальное расстояние, которое упрощенное путь может отклоняться от исходного пути. Чем больше tolerance , тем грубее упрощен путь.
  2. Другим параметром является min_angle который определяет то, что считается поворотной точкой. (Я беру точку поворота, чтобы быть любой точкой исходного пути, угол между входящим и выходящим векторами на упрощенном пути больше min_angle ).

Я буду давать 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) 

введите описание изображения здесь

Наверное, недостаточно хорошо для полностью автоматизированного обнаружения, но это довольно близко к тому, что вы хотели.

Очень интересный вопрос. Вот мое решение, которое допускает переменное разрешение. Хотя, тонкая настройка может быть не простой, так как она в основном предназначена для сужения

Каждый k точек вычисляет выпуклую оболочку и сохраняет ее как множество. Пройдите через не более k точек и удалите любые точки, которые не находятся в выпуклой оболочке, таким образом, чтобы точки не теряли свой первоначальный порядок.

Цель здесь состоит в том, что выпуклая оболочка будет действовать как фильтр, удаляя все «несущественные точки», оставляя только крайние точки. Конечно, если значение k слишком велико, вы получите что-то слишком близко к фактической выпуклой оболочке, а не то, что вы действительно хотите.

Это должно начинаться с маленького k, по крайней мере 4, а затем увеличивать его, пока вы не получите то, что ищете. Вероятно, вы также должны включать только среднюю точку для каждых 3 точек, где угол меньше определенной величины, d. Это обеспечит, чтобы все повороты были не менее d степеней (не реализованы в коде ниже). Однако это, вероятно, следует делать постепенно, чтобы избежать потери информации, так же как и увеличение k-значения. Другим возможным улучшением было бы фактическое повторное выполнение с удалением точек и удаление только точек, которые не были в обоих выпуклых оболочках, хотя для этого требуется более высокое минимальное значение k, равное по меньшей мере 8.

Следующий код, кажется, работает достаточно хорошо, но может по-прежнему использовать улучшения для эффективности и удаления шума. Это также довольно нелогично в определении того, когда он должен остановиться, поэтому код действительно работает (только стоит) от 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 красные точки – это те, которые остаются после фильтра.

Пример выпуклого фильтра

Подход, который вы сделали, звучит многообещающе, но ваши данные сильно передискретизированы. Сначала вы можете отфильтровать координаты x и y, например, с широким гауссовским, а затем с понижением.

В MATLAB вы можете использовать x = conv(x, normpdf(-10 : 10, 0, 5)) а затем x = x(1 : 5 : end) . Вам придется настроить эти числа в зависимости от внутренней сохранности объектов, которые вы отслеживаете, и среднего расстояния между точками.

Тогда вы сможете надежно обнаружить изменения в направлении, используя тот же подход, который вы пробовали раньше, на основе скалярного продукта.

Другая идея – исследовать левое и правое окружение в каждой точке. Это может быть сделано путем создания линейной регрессии N точек до и после каждой точки. Если угол пересечения между точками ниже некоторого порога, то у вас есть угол.

Это можно сделать эффективно, сохраняя очередь точек, находящихся в настоящее время в линейной регрессии, и заменяя старые точки новыми точками, аналогичными среднему среднему.

Вы, наконец, должны объединить смежные углы в один угол. Например, выбирая точку с самым сильным угловым свойством.