Find minimum distance from point to complicated curve

Question

I have a complicated curve defined as a set of points in a table like so (the full table is here):

#  x   y
1.0577  12.0914
1.0501  11.9946
1.0465  11.9338
...

If I plot this table with the commands:

plt.plot(x_data, y_data, c='b',lw=1.)
plt.scatter(x_data, y_data, marker='o', color='k', s=10, lw=0.2)

I get the following:

where I've added the red points and segments manually. What I need is a way to calculate those segments for each of those points, that is: a way to find the minimum distance from a given point in this 2D space to the interpolated curve.

I can't use the distance to the data points themselves (the black dots that generate the blue curve) since they are not located at equal intervals, sometimes they are close and sometimes they are far apart and this deeply affects my results further down the line.

Since this is not a well behaved curve I'm not really sure what I could do. I've tried interpolating it with a UnivariateSpline but it returns a very poor fit:

# Sort data according to x.
temp_data = zip(x_data, y_data)
temp_data.sort()
# Unpack sorted data.
x_sorted, y_sorted = zip(*temp_data)

# Generate univariate spline.
s = UnivariateSpline(x_sorted, y_sorted, k=5)
xspl = np.linspace(0.8, 1.1, 100)
yspl = s(xspl)

# Plot.
plt.scatter(xspl, yspl, marker='o', color='r', s=10, lw=0.2)

I also tried increasing the number of interpolating points but got a mess:

# Sort data according to x.
temp_data = zip(x_data, y_data)
temp_data.sort()
# Unpack sorted data.
x_sorted, y_sorted = zip(*temp_data)

t = np.linspace(0, 1, len(x_sorted))
t2 = np.linspace(0, 1, 100)    
# One-dimensional linear interpolation.
x2 = np.interp(t2, t, x_sorted)
y2 = np.interp(t2, t, y_sorted)
plt.scatter(x2, y2, marker='o', color='r', s=10, lw=0.2)

Any ideas/pointers will be greatly appreciated.

Jamnis

This problem is not convex (en.wikipedia.org/wiki/Convex_optimization)... this means that optimization techniques are generally not guaranteed to give the global minimum. That said, simulated annealing is one option for trying to find global optima in an imperfect world. – Gouda

Answer 1

If you're open to using a library for this, have a look at shapely: https://github.com/Toblerity/Shapely

As a quick example (points.txt contains the data you linked to in your question):

import shapely.geometry as geom
import numpy as np

coords = np.loadtxt('points.txt')

line = geom.LineString(coords)
point = geom.Point(0.8, 10.5)

# Note that "line.distance(point)" would be identical
print(point.distance(line))

As an interactive example (this also draws the line segments you wanted):

import numpy as np
import shapely.geometry as geom
import matplotlib.pyplot as plt

class NearestPoint(object):
    def __init__(self, line, ax):
        self.line = line
        self.ax = ax
        ax.figure.canvas.mpl_connect('button_press_event', self)

    def __call__(self, event):
        x, y = event.xdata, event.ydata
        point = geom.Point(x, y)
        distance = self.line.distance(point)
        self.draw_segment(point)
        print 'Distance to line:', distance

    def draw_segment(self, point):
        point_on_line = line.interpolate(line.project(point))
        self.ax.plot([point.x, point_on_line.x], [point.y, point_on_line.y], 
                     color='red', marker='o', scalex=False, scaley=False)
        fig.canvas.draw()

if __name__ == '__main__':
    coords = np.loadtxt('points.txt')

    line = geom.LineString(coords)

    fig, ax = plt.subplots()
    ax.plot(*coords.T)
    ax.axis('equal')
    NearestPoint(line, ax)
    plt.show()

Note that I've added ax.axis('equal'). shapely operates in the coordinate system that the data is in. Without the equal axis plot, the view will be distorted, and while shapely will still find the nearest point, it won't look quite right in the display:

Osanna

I don't know how I missed this answer until now. It is the most amazing answer I've gotten so far in StackOverflow. Not only did you answer my question, you showed me how to make an interactive plot. I can't thank you enough Joe. – Jamnis

If I might ask, is this going to return the minimum distance to the curve (any point in the curve) or the minimum distance to an existent point (it seems that is the latter)? – Godly

How to do it in 3D (calculate distance from 3D point to 3D spline/curve/line)? – Waterborne

Answer 2

The curve is by nature parametric, i.e. for each x there isn't necessary a unique y and vice versa. So you shouldn't interpolate a function of the form y(x) or x(y). Instead, you should do two interpolations, x(t) and y(t) where t is, say, the index of the corresponding point.

Then you use scipy.optimize.fminbound to find the optimal t such that (x(t) - x0)^2 + (y(t) - y0)^2 is the smallest, where (x0, y0) are the red dots in your first figure. For fminsearch, you could specify the min/max bound for t to be 1 and len(x_data)

Deel

Would you mind clarifying what fminsearch is? Also what you say about making two interpolations, one for x and one for y, wouldn't that be what I tried last in my question which gave me a mess? – Jamnis

don't sort, the initial orders of x and y already in the right sequence. – Deel

also it's 'fminbound', 'fminsearch' is a matlab equivalence. it finds the minimum of a scalar function between two specified bounds. see docs.scipy.org/doc/scipy/reference/generated/… – Deel

The curve quite clearly presents discontinuities, I wouldn't interpolate at all, unless you are going to interpolate only within the non-discontinuous parts, and that's non-trivial to automate. – Geronimo

Answer 3

You could try implementing a calculation of distance from point to line on incremental pairs of points on the curve and finding that minimum. This will introduce a small bit of error from the curve as drawn, but it should be very small, as the points are relatively close together.

http://en.wikipedia.org/wiki/Distance_from_a_point_to_a_line

Schmuck

Moreover, the curve is only defined by the points, not by a specific interpolant of the points. So, unless we assume a specific interpolant function is the right one, we can't even discuss the error made. – Geronimo

Oh, good catch! Yeah, there is no particular error in this method, then. – Schmuck

I had the same idea, but actually it won't work: you must check that the projection of your point belongs to the line. But this may not happen in certain circumstances. Example: your point is just above a "hat" angled curve. Then the minimum distance would be between your point and the upper point of the hat, but you won't be able to find it by orthogonal distance to a line. – Diplegia

If the point on the line is closer than either of the two points in question, you should reject the distance. – Schmuck

Answer 4

You can easily use the package trjtrypy in PyPI: https://pypi.org/project/trjtrypy/

All needed computations and visualizations are available in this package. You can get your answer within a line of code like:

to get the minimum distance use: trjtrypy.basedists.distance(points, curve)

to visualize the curve and points use: trjtrypy.visualizations.draw_landmarks_trajectory(points, curve)

Draftee