Simplify/robustify segment-point distance calculation.

The version in poly_editor is relatively simple because it only supports
array inputs and doesn't vectorize over any input.

The version in proj3d is private so the API can be changed, but it needs
(currently) to at least support non-array inputs and to vectorize over
`p`.

- Rename the parameters to make the difference between the "segment
  ends" (`s0, s1`) and the "point" (`p`) parameters clearer.
- Switch `p` to support (N, ndim) inputs instead of (ndim, N)
  (consistently with most other APIs); adjust test_lines_dists
  accordingly.
- Use vectorized ops everywhere, which also caught the fact that
  previously, entries beyond the third in (what was) `p1, p2` would be
  silently ignored (because access was via `p1[0]`, `p1[1]`, `p2[0]`,
  `p2[1]`).  Instead now the vectorized version naturally extends to any
  number of dimensions.  Adjust format_coord and
  test_lines_dists_nowarning accordingly.
- Also support vectorizing over `s0`, `s1`, if they have the same length
  as `p` (this comes basically for free).

Author: anntzer

examples/event_handling/poly_editor.py

@@ -14,37 +14,24 @@
     You can copy and paste individual parts, or download the entire example
     using the link at the bottom of the page.
 """
+
 import numpy as np
 from matplotlib.lines import Line2D
 from matplotlib.artist import Artist
 
 
-def dist(x, y):
-    """
-    Return the distance between two points.
-    """
-    d = x - y
-    return np.sqrt(np.dot(d, d))
-
-
 def dist_point_to_segment(p, s0, s1):
     """
-    Get the distance of a point to a segment.
-      *p*, *s0*, *s1* are *xy* sequences
-    This algorithm from
-    http://www.geomalgorithms.com/algorithms.html
+    Get the distance from the point *p* to the segment (*s0*, *s1*), where
+    *p*, *s0*, *s1* are ``[x, y]`` arrays.
     """
-    v = s1 - s0
-    w = p - s0
-    c1 = np.dot(w, v)
-    if c1 <= 0:
-        return dist(p, s0)
-    c2 = np.dot(v, v)
-    if c2 <= c1:
-        return dist(p, s1)
-    b = c1 / c2
-    pb = s0 + b * v
-    return dist(p, pb)
+    s01 = s1 - s0
+    s0p = p - s0
+    if (s01 == 0).all():
+        return np.hypot(*s0p)
+    # Project onto segment, without going past segment ends.
+    p1 = s0 + np.clip((s0p @ s01) / (s01 @ s01), 0, 1) * s01
+    return np.hypot(*(p - p1))
 
 
 class PolygonInteractor:

lib/mpl_toolkits/mplot3d/axes3d.py

@@ -1066,7 +1066,7 @@ def format_coord(self, xd, yd):
         # nearest edge
         p0, p1 = min(self._tunit_edges(),
                      key=lambda edge: proj3d._line2d_seg_dist(
-                         edge[0], edge[1], (xd, yd)))
+                         (xd, yd), edge[0][:2], edge[1][:2]))
 
         # scale the z value to match
         x0, y0, z0 = p0

lib/mpl_toolkits/mplot3d/proj3d.py

@@ -6,31 +6,27 @@
 import numpy.linalg as linalg
 
 
-def _line2d_seg_dist(p1, p2, p0):
+def _line2d_seg_dist(p, s0, s1):
     """
-    Return the distance(s) from line defined by p1 - p2 to point(s) p0.
+    Return the distance(s) from point(s) *p* to segment(s) (*s0*, *s1*).
 
-    p0[0] = x(s)
-    p0[1] = y(s)
-
-    intersection point p = p1 + u*(p2-p1)
-    and intersection point lies within segment if u is between 0 and 1.
-
-    If p1 and p2 are identical, the distance between them and p0 is returned.
-    """
-
-    x01 = np.asarray(p0[0]) - p1[0]
-    y01 = np.asarray(p0[1]) - p1[1]
-    if np.all(p1[0:2] == p2[0:2]):
-        return np.hypot(x01, y01)
-
-    x21 = p2[0] - p1[0]
-    y21 = p2[1] - p1[1]
-    u = (x01*x21 + y01*y21) / (x21**2 + y21**2)
-    u = np.clip(u, 0, 1)
-    d = np.hypot(x01 - u*x21, y01 - u*y21)
-
-    return d
+    Parameters
+    ----------
+    p : (ndim,) or (N, ndim) array-like
+        The points from which the distances are computed.
+    s0, s1 : (ndim,) or (N, ndim) array-like
+        The xy(z...) coordinates of the segment endpoints.
+    """
+    s0 = np.asarray(s0)
+    s01 = s1 - s0  # shape (ndim,) or (N, ndim)
+    s0p = p - s0  # shape (ndim,) or (N, ndim)
+    l2 = s01 @ s01  # squared segment length
+    # Avoid div. by zero for degenerate segments (for them, s01 = (0, 0, ...)
+    # so the value of l2 doesn't matter; this just replaces 0/0 by 0/1).
+    l2 = np.where(l2, l2, 1)
+    # Project onto segment, without going past segment ends.
+    p1 = s0 + np.multiply.outer(np.clip(s0p @ s01 / l2, 0, 1), s01)
+    return ((p - p1) ** 2).sum(axis=-1) ** (1/2)
 
 
 def world_transformation(xmin, xmax,

lib/mpl_toolkits/mplot3d/tests/test_axes3d.py

@@ -1120,8 +1120,8 @@ def test_lines_dists():
     ys = (100, 150, 30, 200)
     ax.scatter(xs, ys)
 
-    dist0 = proj3d._line2d_seg_dist(p0, p1, (xs[0], ys[0]))
-    dist = proj3d._line2d_seg_dist(p0, p1, np.array((xs, ys)))
+    dist0 = proj3d._line2d_seg_dist((xs[0], ys[0]), p0, p1)
+    dist = proj3d._line2d_seg_dist(np.array((xs, ys)).T, p0, p1)
     assert dist0 == dist[0]
 
     for x, y, d in zip(xs, ys, dist):
@@ -1133,15 +1133,11 @@ def test_lines_dists():
 
 
 def test_lines_dists_nowarning():
-    # Smoke test to see that no RuntimeWarning is emitted when two first
-    # arguments are the same, see GH#22624
-    p0 = (10, 30, 50)
-    p1 = (10, 30, 20)
-    p2 = (20, 150)
-    proj3d._line2d_seg_dist(p0, p0, p2)
-    proj3d._line2d_seg_dist(p0, p1, p2)
-    p0 = np.array(p0)
-    proj3d._line2d_seg_dist(p0, p0, p2)
+    # No RuntimeWarning must be emitted for degenerate segments, see GH#22624.
+    s0 = (10, 30, 50)
+    p = (20, 150, 180)
+    proj3d._line2d_seg_dist(p, s0, s0)
+    proj3d._line2d_seg_dist(np.array(p), s0, s0)
 
 
 def test_autoscale():