32

I'm trying to come up with an algorithm that will determine turning points in a trajectory of x/y coordinates. The following figures illustrates what I mean: green indicates the starting point and red the final point of the trajectory (the entire trajectory consists of ~ 1500 points): trajectory

In the following figure, I added by hand the possible (global) turning points that an algorithm could return:

trajectory with possible turning points

Obviously, the true turning point is always debatable and will depend on the angle that one specifies that has to lie between points. Furthermore a turning point can be defined on a global scale (what I tried to do with the black circles), but could also be defined on a high-resolution local scale. I'm interested in the global (overall) direction changes, but I'd love to see a discussion on the different approaches that one would use to tease apart global vs local solutions.

What I've tried so far:

  • calculate distance between subsequent points
  • calculate angle between subsequent points
  • look how distance / angle changes between subsequent points

Unfortunately this doesn't give me any robust results. I probably have too calculate the curvature along multiple points, but that's just an idea. I'd really appreciate any algorithms / ideas that might help me here. The code can be in any programming language, matlab or python are preferred.

EDIT here's the raw data (in case somebody want's to play with it):

2
  • Very interesting problem, but I am not sure if this forum is the right place to ask it. I see many subjective ways to define a turning point in the trajectory, so for example on what scale you see it. When you look very closely I can see many different turning points. The way to proceed would be maybe some kind of smoothing of the points on either side of each point (or just draw a line using n points) and make a decision on the angle between those two straight lines. Then you would have 'only' two parameters (n and min. angle), despite the straightening algorithms. Maybe this helps anyway?
    – Alex
    Jan 31, 2013 at 17:51
  • @Alex I'm aware of the subjectivity of this problem. I still think that this might be a problem of general interest and I'd love to see people discuss the different approaches one might use to tease apart local turning points vs global ones.
    – memyself
    Jan 31, 2013 at 17:54

5 Answers 5

31

You could use the Ramer-Douglas-Peucker (RDP) algorithm to simplify the path. Then you could compute the change in directions along each segment of the simplified path. The points corresponding to the greatest change in direction could be called the turning points:

A Python implementation of the RDP algorithm can be found 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()

enter image description here

Two parameters were used above:

  1. The RDP algorithm takes one parameter, the tolerance, which represents the maximum distance the simplified path can stray from the original path. The larger the tolerance, the cruder the simplified path.
  2. The other parameter is the min_angle which defines what is considered a turning point. (I'm taking a turning point to be any point on the original path, whose angle between the entering and exiting vectors on the simplified path is greater than min_angle).
2
  • This seems like what I was trying to achieve with the convex hull approach, as I kinda had quick-hull in mind. +1 and will have to remember this one. My only issue is that it should probably be based on a minimum angle to be considered a turn instead of number of points.
    – Nuclearman
    Jan 31, 2013 at 22:22
  • @MC: Excellent idea. Thank you. (The change has been made.)
    – unutbu
    Jan 31, 2013 at 22:37
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.

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

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.

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.

Not the answer you're looking for? Browse other questions tagged or ask your own question.