obliczyć punkty zwrotne/punkty obrotu w trajektorii (ścieżka)

Question

Próbuję wymyślić algorytm, który określi punkty zwrotne w trajektorii współrzędnych x/y. Poniższe dane przedstawia to, co oznacza, zielonym wskazuje punkt startowy i czerwony punkt końcowy trajektorii (cały tor składa się z ~ 1500 punktów)

Na rysunku I dodaje ręcznie Ewentualny (globalne) punkty zwrotne, że algorytm może wrócić:

Oczywiście, prawdziwy punkt zwrotny jest zawsze dyskusyjna i zależy od kąta, który określa, że ​​jeden musi znajdować się pomiędzy punktami. Ponadto punkt zwrotny można zdefiniować w skali globalnej (co próbowałem zrobić z czarnymi okręgami), ale można go również zdefiniować w lokalnej skali o wysokiej rozdzielczości. Interesują mnie globalne (ogólne) zmiany kierunku, ale chciałbym zobaczyć dyskusję na temat różnych podejść, które można zastosować, by rozdzielić globalne i lokalne rozwiązania.

Co próbowałem dotąd: odległość

• oblicz pomiędzy kolejnymi punktami
• oblicz kąt pomiędzy kolejnymi punktami
• patrzeć jak zmienia się odległość/kąt pomiędzy kolejnymi punktami

Niestety to nie daje żadnych solidnych wyników. Prawdopodobnie również obliczyłem krzywiznę wzdłuż wielu punktów, ale to tylko pomysł. Naprawdę doceniam wszelkie algorytmy/pomysły, które mogą mi w tym pomóc. Kod może być w dowolnym języku programowania, preferowany jest matlab lub python.

EDIT oto dane surowe (w przypadku ktoś chce, aby grać z nim):

• mat file
• text file (współrzędna x pierwszy, współrzędna y w drugiej linii)

Answer 1

Możesz użyć Ramer-Douglas-Peucker (RDP) algorithm, aby uprościć ścieżkę. Następnie możesz obliczyć zmianę kierunków wzdłuż każdego segmentu uproszczonej ścieżki. Punkty odpowiadające największej zmianie kierunku można nazwać punktami zwrotnymi:

Implementację algorytmu RDP w języku Python można znaleźć pod adresem on 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() 

dwa parametry zastosowano powyżej:

1. Algorytm RDP przyjmuje jeden parametr, tolerance, który reprezentuje maksymalną odległość uproszczony ścieżka może zboczyć z oryginalnej ścieżki . Im większy numer tolerance, tym bardziej uproszczona jest ścieżka. Kolejnym parametrem jest min_angle, który określa, co jest punktem zwrotnym. (Punktem zwrotnym jest punkt na oryginalnej ścieżce, którego kąt pomiędzy wektorami wejściowymi i wyjściowymi na uproszczonej ścieżce jest większy niż min_angle).

Answer 2

Podejście, które podjąłeś, brzmi obiecująco, ale twoje dane są mocno nadpróbkowane. Możesz najpierw filtrować współrzędne X i Y, na przykład szerokim krzywa Gaussa, a następnie w dół.

W programie MATLAB można użyć x = conv(x, normpdf(-10 : 10, 0, 5)), a następnie x = x(1 : 5 : end). Będziesz musiał poprawić te liczby w zależności od wewnętrznej trwałości obiektów, które śledzisz i średniej odległości między punktami.

Wówczas będziesz w stanie bardzo dokładnie wykryć zmiany kierunku, używając tej samej metody, którą próbowałem wcześniej, na podstawie produktu skalarnego, jak sobie wyobrażam.

Answer 3

Bardzo interesujące pytanie. Oto moje rozwiązanie, które pozwala na różną rozdzielczość. Mimo że precyzyjne dostrojenie może nie być proste, ponieważ głównie ma na celu zawężenie każdego k, obliczyć wypukły kadłub i zapisać go jako zestaw. Przejdź przez co najwyżej k punktów i usuń punkty, które nie znajdują się w wypukłym kadłubie, w taki sposób, aby punkty nie utraciły pierwotnej kolejności.

Celem jest tutaj, aby wypukły kadłub działał jak filtr, usuwając wszystkie "nieistotne punkty" pozostawiając tylko skrajne punkty. Oczywiście, jeśli wartość k jest zbyt wysoka, otrzymasz coś zbyt zbliżonego do rzeczywistego wypukłego kadłuba, zamiast tego, czego naprawdę chcesz.

Powinno się zaczynać od małego k, co najmniej 4, a następnie zwiększać, aż otrzymasz to, czego szukasz. Powinieneś również prawdopodobnie uwzględniać tylko środkowy punkt za każde 3 punkty, w których kąt jest poniżej określonej wartości, d. Zapewni to, że wszystkie zwoje są co najmniej d stopni (nie zaimplementowane w kodzie poniżej). Jednak powinno to być wykonywane stopniowo, aby uniknąć utraty informacji, tak samo jak zwiększenie wartości k. Inną możliwą poprawą byłoby ponowne uruchomienie z usuniętymi punktami i tylko usunięcie punktów, które nie były w obu wypukłych kadłubach, jednak wymaga to wyższej minimalnej wartości k co najmniej 8.

Poniższy kod Wydaje się, że działa całkiem dobrze, ale nadal może używać ulepszeń w zakresie wydajności i usuwania hałasu. Jest również dość nieelegancka w określaniu, kiedy powinna się zatrzymać, tak więc kod działa tylko (tak jak stoi) od około k = 4 do 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) 

Oto zrzut ekranu powyższego kodu z k = 14. 61 czerwonych kropek to te, które pozostają po filtrze.

Answer 4

Innym pomysłem jest zbadanie lewy i prawy otoczenie w każdym punkcie. Można tego dokonać, tworząc liniową regresję N punktów przed i po każdym punkcie. Jeśli kąt przecięcia między punktami jest poniżej pewnego progu, to masz róg.

Można to zrobić skutecznie, utrzymując kolejkę punktów będących obecnie w regresji liniowej i zastępując stare punkty nowymi punktami, podobnymi do średniej bieżącej.

Musisz w końcu połączyć sąsiadujące narożniki z jednym rogiem. Na przykład. wybierając punkt z najsilniejszą właściwością narożną.

Answer 5

Podam poniżej kod numpy/scipy, ponieważ nie mam prawie żadnego doświadczenia z Matlabą.

Jeśli twoja krzywa jest wystarczająco gładka, możesz zidentyfikować swoje punkty zwrotne jako najwyższe wartości curvature. Biorąc numer indeksu punktu jako parametr krzywej, a central differences scheme można obliczyć krzywiznę z następującego kodu

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) 

Prawdopodobnie będzie chciał wygładzić krzywą w pierwszej kolejności, a następnie obliczyć krzywiznę, a następnie zidentyfikować najwyższa punkty krzywizny. Poniższa funkcja właśnie to robi:

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() 

Niektóre tłumacząc to w kolejności:

• turning_points jest liczba punktów, które chcesz zidentyfikować
• smoothing_radius jest promień splotu wygładzającym być stosowane do twoich danych przed obliczeniem krzywizny
• cluster_radius to odległość od punktu wysokiej krzywizny wybranego jako punkt zwrotny, w którym żaden inny punkt nie powinien być uważany za kandydata.

Być może trzeba się bawić z parametrami trochę, ale mam coś takiego:

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

prawdopodobnie nie wystarczająco dobry, w pełni zautomatyzowane wykrywanie, ale to dość blisko tego, co chciałeś.