calcolare punti di svolta/punti di articolazione in traiettoria (percorso)

Question

Sto provando a elaborare un algoritmo che determinerà i punti di svolta in una traiettoria di coordinate x/y. Le seguenti figure illustrano cosa intendo: verde indica il punto di partenza e rosso l'ultimo punto della traiettoria (l'intera traiettoria consiste di ~ 1500 punti):

Nella figura seguente, sono aggiunti dai mano possibile () svolte globali che un algoritmo potrebbe tornare:

Ovviamente, la vera svolta è sempre discutibile e dipenderà l'angolo che si specifica che deve trovano tra i punti. Inoltre un punto di svolta può essere definito su scala globale (cosa ho cercato di fare con i cerchi neri), ma potrebbe anche essere definito su una scala locale ad alta risoluzione. Sono interessato ai cambiamenti globali (generali) della direzione, ma mi piacerebbe vedere una discussione sui diversi approcci che si potrebbero utilizzare per separare le soluzioni globali e locali.

Quello che ho provato finora: distanza

• calcolare tra i punti successivi
• angolo calcolare tra i punti successivi
• sguardo come i cambiamenti distanza/angolo tra i punti successivi

Purtroppo questo non mi dà risultati robusti. Probabilmente ho calcolato anche la curvatura su più punti, ma è solo un'idea. Apprezzerei davvero qualsiasi algoritmo/idea che potrebbe aiutarmi qui. Il codice può essere in qualsiasi linguaggio di programmazione, matlab o python sono preferiti.

EDIT ecco i dati grezzi (nel caso qualcuno voglia di giocare con lui):

• mat file
• text file (coordinata x prima, coordinata y in seconda linea)

Answer 1

È possibile utilizzare Ramer-Douglas-Peucker (RDP) algorithm per semplificare il percorso. Quindi è possibile calcolare la modifica delle direzioni lungo ciascun segmento del percorso semplificato. I punti corrispondenti al più grande cambiamento di direzione potrebbero essere chiamati i punti di svolta:

Un'implementazione Python dell'algoritmo RDP può essere trovata 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() 

Due parametri sono stati usati in precedenza:

1. L'algoritmo RDP accetta un parametro, il tolerance, che rappresenta la distanza massima tracciato semplificato può deviare dal percorso originale . Più grande è il tolerance, il più semplice è il percorso semplificato.
2. L'altro parametro è min_angle che definisce ciò che è considerato un punto di svolta. (Sto prendendo un punto di svolta per essere qualsiasi punto sul percorso originale, il cui angolo tra i vettori in entrata e in uscita sul percorso semplificato è maggiore di min_angle).

Answer 2

L'approccio che hai preso sembra promettente, ma i tuoi dati sono pesantemente sovracampionati. È possibile filtrare prima le coordinate xey, ad esempio con un gaussiano ampio e quindi un downsampling.

In MATLAB, è possibile utilizzare x = conv(x, normpdf(-10 : 10, 0, 5)) e quindi x = x(1 : 5 : end). Dovrai modificare questi numeri a seconda della persistenza intrinseca degli oggetti che stai monitorando e della distanza media tra i punti.

Quindi, sarete in grado di rilevare i cambiamenti di direzione in modo molto affidabile, usando lo stesso approccio che avete provato prima, basato sul prodotto scalare, immagino.

Answer 3

Una domanda molto interessante. Ecco la mia soluzione, che consente la risoluzione variabile. Sebbene, la regolazione fine potrebbe non essere semplice, in quanto è principalmente intesa a restringere

Ogni punto k, calcolare lo scafo convesso e memorizzarlo come un set. Passare attraverso al massimo k punti e rimuovere tutti i punti che non sono nello scafo convesso, in modo tale che i punti non perdano il loro ordine originale.

Lo scopo qui è che lo scafo convesso fungerà da filtro, rimuovendo tutti i "punti non importanti" lasciando solo i punti estremi. Ovviamente, se il valore k è troppo alto, ti ritroverai con qualcosa di troppo vicino allo scafo convesso reale, invece di ciò che effettivamente desideri.

Questo dovrebbe iniziare con un piccolo k, almeno 4, quindi aumentarlo fino a ottenere quello che cerchi. Probabilmente dovresti anche includere solo il punto medio per ogni 3 punti in cui l'angolo è inferiore ad una certa quantità, d. Ciò garantirebbe che tutti i turni siano almeno d gradi (non implementati nel codice di seguito). Tuttavia, questo dovrebbe probabilmente essere fatto in modo incrementale per evitare la perdita di informazioni, come aumentare il valore k. Un altro possibile miglioramento potrebbe essere quello di eseguire nuovamente i punti rimossi e rimuovere solo i punti che non erano in entrambi gli scafi convessi, sebbene ciò richieda un valore minimo minimo di almeno 8.

Il seguente codice sembra funzionare abbastanza bene, ma potrebbe ancora utilizzare miglioramenti per l'efficienza e la rimozione del rumore. È anche piuttosto inelegante nel determinare quando dovrebbe fermarsi, quindi il codice funziona davvero solo (così com'è) da circa k = 4 a 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) 

Ecco una schermata del codice precedente con k = 14. I 61 punti rossi sono quelli che rimangono dopo il filtro.

Answer 4

Un'altra idea è quella di esaminare il circostante giusti in ogni punto e sinistro. Questo può essere fatto creando una regressione lineare di N punti prima e dopo ogni punto. Se l'angolo di intersezione tra i punti è inferiore a qualche soglia, allora hai un angolo.

Questo può essere fatto in modo efficiente mantenendo una coda dei punti correntemente nella regressione lineare e sostituendo i vecchi punti con nuovi punti, simili a una media corrente.

Infine, è necessario unire gli angoli adiacenti in un singolo angolo. Per esempio. scegliendo il punto con la proprietà angolo più forte.

Answer 5

Fornirò codice numpy/scipy di seguito, poiché non ho quasi nessuna esperienza di Matlab.

Se la curva è abbastanza regolare, è possibile identificare i punti di svolta come quelli di massimo curvature. Prendendo il numero di indice punto come parametro della curva, e un central differences scheme, si può calcolare la curvatura con il seguente codice

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) 

Probabilmente si desidera per spianare la curva per primo, quindi calcolare la curvatura, quindi identificare il più alto punti di curvatura. La seguente funzione fa proprio questo:

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

qualche spiegazione è in ordine:

• turning_points è il numero di punti che si desidera identificare
• smoothing_radius è il raggio di un convoluzione lisciatura da applicare ai tuoi dati prima di calcolare la curvatura
• cluster_radius è la distanza da un punto di alta curvatura selezionato come un punto di svolta in cui nessun altro punto deve essere considerato come candidato.

Potrebbe essere necessario giocare con i parametri un po ', ma ho avuto qualcosa di simile:

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

Probabilmente non è abbastanza buono per una rilevazione completamente automatizzato, ma è abbastanza vicino a quello che volevi.