6

I'm looking for a function to compute the euclidian distance between a numpy array of points with two coordinates (x, y) and a line segment. My goal is to have a result in under 0.01 sec for a line segment and 10k points.

I already found a function for a single point. But running a for loop is very inefficient.

I also found this function that calculates the distance to the infinite line:

def line_dists(points, start, end):
    if np.all(start == end):
        return np.linalg.norm(points - start, axis=1)

    vec = end - start
    cross = np.cross(vec, start - points)
    return np.divide(abs(cross), np.linalg.norm(vec))

It is very efficient and I would like to have a similar approach for a bounded line.

Thank you for your help.

meowgoesthedog
  • 14,670
  • 4
  • 27
  • 40
Jérôme
  • 142
  • 2
  • 8
  • 2
    Please share the input data so its easier to answer – yatu Jan 30 '19 at 13:49
  • As a side not, for 2-dimensional vector, you will earn a tremendous amount of time, rewriting a small function that computes the norm, directly from the coordinates of the vector. `np.linalg.norm` is very slow compared to direct computation for vectors of low dimensions. The same holds for `np.cross` – Eskapp Jan 30 '19 at 14:17

1 Answers1

10

Setup – test point P, endpoints A and B:

enter image description here

  • Take the dot-product of P - A with normalize(A - B) to obtain the signed parallel distance component s from A. Likewise with B and t.

  • Take the maximum of these two numbers and zero to get the clamped parallel distance component. This will only be non-zero if the point is outside the "boundary" (Voronoi region?) of the segment.

  • Calculate the perpendicular distance component as before, using the cross-product.

  • Use Pythagoras to compute the required closest distance (gray line from P to A).

The above is branchless and thus easy to vectorize with numpy:

def lineseg_dists(p, a, b):
    # Handle case where p is a single point, i.e. 1d array.
    p = np.atleast_2d(p)

    # TODO for you: consider implementing @Eskapp's suggestions
    if np.all(a == b):
        return np.linalg.norm(p - a, axis=1)

    # normalized tangent vector
    d = np.divide(b - a, np.linalg.norm(b - a))

    # signed parallel distance components
    s = np.dot(a - p, d)
    t = np.dot(p - b, d)

    # clamped parallel distance
    h = np.maximum.reduce([s, t, np.zeros(len(p))])

    # perpendicular distance component, as before
    # note that for the 3D case these will be vectors
    c = np.cross(p - a, d)

    # use hypot for Pythagoras to improve accuracy
    return np.hypot(h, c)
Kel Solaar
  • 3,660
  • 1
  • 23
  • 29
meowgoesthedog
  • 14,670
  • 4
  • 27
  • 40
  • 3
    Thank you so much for your quick and detailed solution! This seems to work as expected and in around 1ms for 10k points. – Jérôme Jan 30 '19 at 14:42