code to compute bezier segment / path lengths

brunobeltran · brunobeltran · commit 23ffc9f7099f · 2020-04-06T12:29:29.000-07:00
diff --git a/lib/matplotlib/bezier.py b/lib/matplotlib/bezier.py
@@ -4,6 +4,7 @@
 
 import math
 import warnings
+from collections import deque
 
 import numpy as np
 
@@ -219,6 +220,96 @@ def point_at_t(self, t):
         """Evaluate curve at a single point *t*. Returns a Tuple[float*d]."""
         return tuple(self(t))
 
+    def split_at_t(self, t):
+        """Split into two Bezier curves using de casteljau's algorithm.
+
+        Parameters
+        ----------
+        t : float
+            Point in [0,1] at which to split into two curves
+
+        Returns
+        -------
+        B1, B2 : BezierSegment
+            The two sub-curves.
+        """
+        new_cpoints = split_de_casteljau(self._cpoints, t)
+        return BezierSegment(new_cpoints[0]), BezierSegment(new_cpoints[1])
+
+    def control_net_length(self):
+        """Sum of lengths between control points"""
+        L = 0
+        N, d = self._cpoints.shape
+        for i in range(N - 1):
+            L += np.linalg.norm(self._cpoints[i+1] - self._cpoints[i])
+        return L
+
+    def arc_length(self, rtol=None, atol=None):
+        """Estimate the length using iterative refinement.
+
+        Our estimate is just the average between the length of the chord and
+        the length of the control net.
+
+        Since the chord length and control net give lower and upper bounds
+        (respectively) on the length, this maximum possible error is tested
+        against an absolute tolerance threshold at each subdivision.
+
+        However, sometimes this estimator converges much faster than this error
+        esimate would suggest. Therefore, the relative change in the length
+        estimate between subdivisions is compared to a relative error tolerance
+        after each set of subdivisions.
+
+        Parameters
+        ----------
+        rtol : float, default 1e-4
+            If :code:`abs(est[i+1] - est[i]) <= rtol * est[i+1]`, we return
+            :code:`est[i+1]`.
+        atol : float, default 1e-6
+            If the distance between chord length and control length at any
+            point falls below this number, iteration is terminated.
+        """
+        if rtol is None:
+            rtol = 1e-4
+        if atol is None:
+            atol = 1e-6
+
+        chord = np.linalg.norm(self._cpoints[-1] - self._cpoints[0])
+        net = self.control_net_length()
+        max_err = (net - chord)/2
+        curr_est = chord + max_err
+        # early exit so we don't try to "split" paths of zero length
+        if max_err < atol:
+            return curr_est
+
+        prev_est = np.inf
+        curves = deque([self])
+        errs = deque([max_err])
+        lengths = deque([curr_est])
+        while np.abs(curr_est - prev_est) > rtol * curr_est:
+            # subdivide the *whole* curve before checking relative convergence
+            # again
+            prev_est = curr_est
+            num_curves = len(curves)
+            for i in range(num_curves):
+                curve = curves.popleft()
+                new_curves = curve.split_at_t(0.5)
+                max_err -= errs.popleft()
+                curr_est -= lengths.popleft()
+                for ncurve in new_curves:
+                    chord = np.linalg.norm(
+                            ncurve._cpoints[-1] - ncurve._cpoints[0])
+                    net = ncurve.control_net_length()
+                    nerr = (net - chord)/2
+                    nlength = chord + nerr
+                    max_err += nerr
+                    curr_est += nlength
+                    curves.append(ncurve)
+                    errs.append(nerr)
+                    lengths.append(nlength)
+                if max_err < atol:
+                    return curr_est
+        return curr_est
+
     def arc_area(self):
         r"""
         (Signed) area swept out by ray from origin to curve.
diff --git a/lib/matplotlib/path.py b/lib/matplotlib/path.py
@@ -642,6 +642,27 @@ def intersects_bbox(self, bbox, filled=True):
         return _path.path_intersects_rectangle(
             self, bbox.x0, bbox.y0, bbox.x1, bbox.y1, filled)
 
+    def length(self, rtol=None, atol=None, **kwargs):
+        r"""Get length of Path.
+
+        Equivalent to (but not computed as)
+
+        .. math::
+
+            \sum_{j=1}^N \int_0^1 ||B'_j(t)|| dt
+
+        where the sum is over the :math:`N` Bezier curves that comprise the
+        Path. Notice that this measure of length will assign zero weight to all
+        isolated points on the Path.
+
+        Returns
+        -------
+        length : float
+            The path length.
+        """
+        return np.sum([B.arc_length(rtol, atol)
+                       for B, code in self.iter_bezier(**kwargs)])
+
     def signed_area(self, **kwargs):
         """
         Get signed area of the filled path.
diff --git a/lib/matplotlib/tests/test_bezier.py b/lib/matplotlib/tests/test_bezier.py
@@ -9,6 +9,13 @@
 _test_curves = [list(tc.path.iter_bezier())[-1][0] for tc in _test_curves]
 
 
+def _integral_arc_length(B):
+    dB = B.differentiate(B)
+    def integrand(t):
+        return np.linalg.norm(dB(t))
+    return scipy.integrate.quad(integrand, 0, 1)[0]
+
+
 def _integral_arc_area(B):
     """(Signed) area swept out by ray from origin to curve."""
     dB = B.differentiate(B)
@@ -22,3 +29,9 @@ def integrand(t):
 def test_area_formula():
     for B in _test_curves:
         assert(np.isclose(_integral_arc_area(B), B.arc_area()))
+
+
+def test_length_iteration():
+    for B in _test_curves:
+        assert(np.isclose(_integral_arc_length(B), B.arc_length(
+                rtol=1e-5, atol=1e-8), rtol=1e-5, atol=1e-8))
diff --git a/lib/matplotlib/tests/test_path.py b/lib/matplotlib/tests/test_path.py
@@ -14,28 +14,27 @@
 from matplotlib.backend_bases import MouseEvent
 
 
-ExampleCurve = namedtuple('ExampleCurve', ['path', 'extents', 'area'])
+ExampleCurve = namedtuple('ExampleCurve',
+                          ['path', 'extents', 'area', 'length'])
 _test_curves = [
     # interior extrema determine extents and degenerate derivative
     ExampleCurve(Path([[0, 0], [1, 0], [1, 1], [0, 1]],
                        [Path.MOVETO, Path.CURVE4, Path.CURVE4, Path.CURVE4]),
-                 (0., 0., 0.75, 1.), 0.6),
+                 (0., 0., 0.75, 1.), area=0.6, length=2.0),
     # a quadratic curve, clockwise
     ExampleCurve(Path([[0, 0], [0, 1], [1, 0]], [Path.MOVETO, Path.CURVE3,
-                 Path.CURVE3]), (0., 0., 1., 0.5), -1/3),
+                 Path.CURVE3]), (0., 0., 1., 0.5), area=-1/3, length=(1/25)*(10
+                 + 15*np.sqrt(2) + np.sqrt(5)*(np.arcsinh(2) + np.arcsinh(3)))
+                 ),
     # a linear curve, degenerate vertically
     ExampleCurve(Path([[0, 1], [1, 1]], [Path.MOVETO, Path.LINETO]),
-                 (0., 1., 1., 1.), 0.),
+                 (0., 1., 1., 1.), area=0., length=1.),
     # a point
-    ExampleCurve(Path([[1, 2]], [Path.MOVETO]), (1., 2., 1., 2.), 0.),
+    ExampleCurve(Path([[1, 2]], [Path.MOVETO]), (1., 2., 1., 2.), area=0.,
+                 length=0.),
 ]
 
 
-# convenient curve (xmax=0.75, area=0.6, length=2.0)
-_hard_curve = Path([[0, 0], [1, 0], [1, 1], [0, 1]],
-                    [Path.MOVETO, Path.CURVE4, Path.CURVE4, Path.CURVE4])
-
-
 def test_empty_closed_path():
     path = Path(np.zeros((0, 2)), closed=True)
     assert path.vertices.shape == (0, 2)
@@ -101,6 +100,12 @@ def test_signed_area_unit_rectangle():
     assert np.isclose(rect.signed_area(), 1)
 
 
+def test_length_curve():
+    for test_curve in _test_curves:
+        length, path = test_curve.length, test_curve.path
+    assert(np.isclose(path.length(), length))
+
+
 def test_point_in_path_nan():
     box = np.array([[0, 0], [1, 0], [1, 1], [0, 1], [0, 0]])
     p = Path(box)
