The Python Oracle

calculate turning points / pivot points in trajectory (path)

This video explains
calculate turning points / pivot points in trajectory (path)

--

Become part of the top 3% of the developers by applying to Toptal
https://topt.al/25cXVn

--

Music by Eric Matyas
https://www.soundimage.org
Track title: Puzzle Game 3 Looping

--

Chapters
00:00 Question
02:31 Accepted answer (Score 29)
04:35 Answer 2 (Score 9)
06:52 Answer 3 (Score 4)
09:39 Answer 4 (Score 1)
10:29 Thank you

--

Full question
https://stackoverflow.com/questions/1463...

Question links:
[mat file]: http://www.speedyshare.com/8CgN2/bla.mat
[text file]: http://www.speedyshare.com/nx8Er/bla.dat...

Accepted answer links:
[Ramer-Douglas-Peucker (RDP) algorithm]: http://en.wikipedia.org/wiki/Ramer-Dougl...
[on github]: https://github.com/sebleier/RDP/

Answer 2 links:
[curvature]: http://en.wikipedia.org/wiki/Curvature#L...
[central differences scheme]: http://en.wikipedia.org/wiki/Finite_diff...

--

Content licensed under CC BY-SA
https://meta.stackexchange.com/help/lice...

--

Tags
#python #algorithm #matlab #signalprocessing

#avk47



ANSWER 1

Score 11


I will be giving numpy/scipy code below, as I have almost no Matlab experience.

If your curve is smooth enough, you could identify your turning points as those of highest curvature. Taking the point index number as the curve parameter, and a central differences scheme, you can compute the curvature with the following code

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)

You will probably want to smooth your curve out first, then calculate the curvature, then identify the highest curvature points. The following function does just that:

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

Some explaining is in order:

  • turning_points is the number of points you want to identify
  • smoothing_radius is the radius of a smoothing convolution to be applied to your data before computing the curvature
  • cluster_radius is the distance from a point of high curvature selected as a turning point where no other point should be considered as a candidate.

You may have to play around with the parameters a little, but I got something like this:

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

enter image description here

Probably not good enough for a fully automated detection, but it's pretty close to what you wanted.




ANSWER 2

Score 4


A very interesting question. Here is my solution, that allows for variable resolution. Although, fine-tuning it may not be simple, as it's mostly intended to narrow down

Every k points, calculate the convex hull and store it as a set. Go through the at most k points and remove any points that are not in the convex hull, in such a way that the points don't lose their original order.

The purpose here is that the convex hull will act as a filter, removing all of "unimportant points" leaving only the extreme points. Of course, if the k-value is too high, you'll end up with something too close to the actual convex hull, instead of what you actually want.

This should start with a small k, at least 4, then increase it until you get what you seek. You should also probably only include the middle point for every 3 points where the angle is below a certain amount, d. This would ensure that all of the turns are at least d degrees (not implemented in code below). However, this should probably be done incrementally to avoid loss of information, same as increasing the k-value. Another possible improvement would be to actually re-run with points that were removed, and and only remove points that were not in both convex hulls, though this requires a higher minimum k-value of at least 8.

The following code seems to work fairly well, but could still use improvements for efficiency and noise removal. It's also rather inelegant in determining when it should stop, thus the code really only works (as it stands) from around k=4 to 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)

Here is a screen shot of the above code with k=14. The 61 red dots are the ones that remain after the filter.

Convex Filter Example




ANSWER 3

Score 1


The approach you took sounds promising but your data is heavily oversampled. You could filter the x and y coordinates first, for example with a wide Gaussian and then downsample.

In MATLAB, you could use x = conv(x, normpdf(-10 : 10, 0, 5)) and then x = x(1 : 5 : end). You will have to tweak those numbers depending on the intrinsic persistence of the objects you are tracking and the average distance between points.

Then, you will be able to detect changes in direction very reliably, using the same approach you tried before, based on the scalar product, I imagine.




ANSWER 4

Score 0


Another idea is to examine the left and the right surroundings at every point. This may be done by creating a linear regression of N points before and after each point. If the intersecting angle between the points is below some threshold, then you have an corner.

This may be done efficiently by keeping a queue of the points currently in the linear regression and replacing old points with new points, similar to a running average.

You finally have to merge adjacent corners to a single corner. E.g. choosing the point with the strongest corner property.