correctly compute bounding box for path

brunobeltran · brunobeltran · commit b29f403cf5ab · 2020-04-08T17:32:24.000-07:00
diff --git a/doc/users/next_whats_new/2020-03-31-path-size-methods.rst b/doc/users/next_whats_new/2020-03-31-path-size-methods.rst
@@ -0,0 +1,12 @@
+
+Path extents now computed correctly
+-----------------------------------
+
+Historically, `~.path.Path.get_extents` has always simply returned the Bbox of
+its control points. While this is a correct upper bound for the path's extents,
+it can differ dramatically from the Path's actual extents for non-linear Bezier
+curves.
+
+In addition, `~.bezier.BezierSegment` has gained more documentation and
+usability improvements, including properties that contain its dimension, degree,
+control_points, and more.
diff --git a/lib/matplotlib/bezier.py b/lib/matplotlib/bezier.py
@@ -3,11 +3,19 @@
 """
 
 import math
+import warnings
 
 import numpy as np
 
 import matplotlib.cbook as cbook
 
+# same algorithm as 3.8's math.comb
+@np.vectorize
+def _comb(n, k):
+    k = min(k, n - k)
+    i = np.arange(1, k + 1)
+    return np.prod((n + 1 - i)/i).astype(int)
+
 
 class NonIntersectingPathException(ValueError):
     pass
@@ -168,27 +176,113 @@ def find_bezier_t_intersecting_with_closedpath(
 
 class BezierSegment:
     """
-    A D-dimensional Bezier segment.
+    A d-dimensional Bezier segment.
 
     Parameters
     ----------
-    control_points : (N, D) array
+    control_points : (N, d) array
         Location of the *N* control points.
     """
 
     def __init__(self, control_points):
-        n = len(control_points)
-        self._orders = np.arange(n)
-        coeff = [math.factorial(n - 1)
-                 // (math.factorial(i) * math.factorial(n - 1 - i))
-                 for i in range(n)]
-        self._px = np.asarray(control_points).T * coeff
+        self._cpoints = np.asarray(control_points)
+        self._N, self._d = self._cpoints.shape
+        self._orders = np.arange(self._N)
+        coeff = [math.factorial(self._N - 1)
+                 // (math.factorial(i) * math.factorial(self._N - 1 - i))
+                 for i in range(self._N)]
+        self._px = self._cpoints.T * coeff
 
     def point_at_t(self, t):
         """Return the point on the Bezier curve for parameter *t*."""
         return tuple(
             self._px @ (((1 - t) ** self._orders)[::-1] * t ** self._orders))
 
+    @property
+    def control_points(self):
+        """The control points of the curve."""
+        return self._cpoints
+
+    @property
+    def dimension(self):
+        """The dimension of the curve."""
+        return self._d
+
+    @property
+    def degree(self):
+        """The number of control points in the curve."""
+        return self._N - 1
+
+    @property
+    def polynomial_coefficients(self):
+        r"""
+        The polynomial coefficients of the Bezier curve.
+
+        Returns
+        -------
+        float, (n+1, d) array_like
+            Coefficients after expanding in polynomial basis, where :math:`n`
+            is the degree of the bezier curve and :math:`d` its dimension.
+            These are the numbers (:math:`C_j`) such that the curve can be
+            written :math:`\sum_{j=0}^n C_j t^j`.
+
+        Notes
+        -----
+        The coefficients are calculated as
+
+        .. math::
+
+            {n \choose j} \sum_{i=0}^j (-1)^{i+j} {j \choose i} P_i
+
+        where :math:`P_i` are the control points of the curve.
+        """
+        n = self.degree
+        # matplotlib uses n <= 4. overflow plausible starting around n = 15.
+        if n > 10:
+            warnings.warn("Polynomial coefficients formula unstable for high "
+                          "order Bezier curves!", RuntimeWarning)
+        d = self.dimension
+        P = self.control_points
+        coefs = np.zeros((n+1, d))
+        for j in range(n+1):
+            i = np.arange(j+1)
+            prefactor = np.power(-1, i + j) * _comb(j, i)
+            coefs[j] = _comb(n, j) * np.sum(prefactor[:, None]*P[i], axis=0)
+        return coefs
+
+    @property
+    def axis_aligned_extrema(self):
+        """
+        Return the dimension and location of the curve's interior extrema.
+
+        The extrema are the points along the curve where one of its partial
+        derivatives is zero.
+
+        Returns
+        -------
+        dims : int, array_like
+            Index :math:`i` of the partial derivative which is zero at each
+            interior extrema.
+        dzeros : float, array_like
+            Of same size as dims. The :math:`t` such that :math:`d/dx_i B(t) =
+            0`
+        """
+        n = self.degree
+        Cj = self.polynomial_coefficients
+        dCj = np.arange(1, n+1)[:, None] * Cj[1:]
+        if len(dCj) == 0:
+            return np.array([]), np.array([])
+        dims = []
+        roots = []
+        for i, pi in enumerate(dCj.T):
+            r = np.roots(pi[::-1])
+            roots.append(r)
+            dims.append(np.full_like(r, i))
+        roots = np.concatenate(roots)
+        dims = np.concatenate(dims)
+        in_range = np.isreal(roots) & (roots >= 0) & (roots <= 1)
+        return dims[in_range], np.real(roots)[in_range]
+
 
 def split_bezier_intersecting_with_closedpath(
         bezier, inside_closedpath, tolerance=0.01):
diff --git a/lib/matplotlib/path.py b/lib/matplotlib/path.py
@@ -17,6 +17,18 @@
 import matplotlib as mpl
 from . import _path, cbook
 from .cbook import _to_unmasked_float_array, simple_linear_interpolation
+from .bezier import BezierSegment
+
+
+def _update_extents(extents, point):
+    dim = len(point)
+    for i, xi in enumerate(point):
+        if xi < extents[i]:
+            extents[i] = xi
+        # elif here would fail to correctly update from "null" extents of
+        # np.array([np.inf, np.inf, -np.inf, -np.inf])
+        if extents[i+dim] < xi:
+            extents[i+dim] = xi
 
 
 class Path:
@@ -420,6 +432,53 @@ def iter_segments(self, transform=None, remove_nans=True, clip=None,
                     curr_vertices = np.append(curr_vertices, next(vertices))
             yield curr_vertices, code
 
+    def iter_bezier(self, **kwargs):
+        """
+        Iterate over each bezier curve (lines included) in a Path.
+
+        Parameters
+        ----------
+        **kwargs
+            Forwarded to `.iter_segments`.
+
+        Yields
+        ------
+        B : matplotlib.bezier.BezierSegment
+            The bezier curves that make up the current path. Note in particular
+            that freestanding points are bezier curves of order 0, and lines
+            are bezier curves of order 1 (with two control points).
+        code : Path.code_type
+            The code describing what kind of curve is being returned.
+            Path.MOVETO, Path.LINETO, Path.CURVE3, Path.CURVE4 correspond to
+            bezier curves with 1, 2, 3, and 4 control points (respectively).
+            Path.CLOSEPOLY is a Path.LINETO with the control points correctly
+            chosen based on the start/end points of the current stroke.
+        """
+        first_vert = None
+        prev_vert = None
+        for verts, code in self.iter_segments(**kwargs):
+            if first_vert is None:
+                if code != Path.MOVETO:
+                    raise ValueError("Malformed path, must start with MOVETO.")
+            if code == Path.MOVETO:  # a point is like "CURVE1"
+                first_vert = verts
+                yield BezierSegment(np.array([first_vert])), code
+            elif code == Path.LINETO:  # "CURVE2"
+                yield BezierSegment(np.array([prev_vert, verts])), code
+            elif code == Path.CURVE3:
+                yield BezierSegment(np.array([prev_vert, verts[:2],
+                                              verts[2:]])), code
+            elif code == Path.CURVE4:
+                yield BezierSegment(np.array([prev_vert, verts[:2],
+                                              verts[2:4], verts[4:]])), code
+            elif code == Path.CLOSEPOLY:
+                yield BezierSegment(np.array([prev_vert, first_vert])), code
+            elif code == Path.STOP:
+                return
+            else:
+                raise ValueError("Invalid Path.code_type: " + str(code))
+            prev_vert = verts[-2:]
+
     @cbook._delete_parameter("3.3", "quantize")
     def cleaned(self, transform=None, remove_nans=False, clip=None,
                 quantize=False, simplify=False, curves=False,
@@ -528,22 +587,39 @@ def contains_path(self, path, transform=None):
             transform = transform.frozen()
         return _path.path_in_path(self, None, path, transform)
 
-    def get_extents(self, transform=None):
+    def get_extents(self, transform=None, **kwargs):
         """
-        Return the extents (*xmin*, *ymin*, *xmax*, *ymax*) of the path.
+        Get Bbox of the path.
 
-        Unlike computing the extents on the *vertices* alone, this
-        algorithm will take into account the curves and deal with
-        control points appropriately.
+        Parameters
+        ----------
+        transform : matplotlib.transforms.Transform, optional
+            Transform to apply to path before computing extents, if any.
+        **kwargs
+            Forwarded to `.iter_bezier`.
+
+        Returns
+        -------
+        matplotlib.transforms.Bbox
+            The extents of the path Bbox([[xmin, ymin], [xmax, ymax]])
         """
         from .transforms import Bbox
-        path = self
         if transform is not None:
-            transform = transform.frozen()
-            if not transform.is_affine:
-                path = self.transformed(transform)
-                transform = None
-        return Bbox(_path.get_path_extents(path, transform))
+            self = transform.transform_path(self)
+        # return value for empty paths to match what used to be done in _path.h
+        extents = np.array([np.inf, np.inf, -np.inf, -np.inf])
+        for curve, code in self.iter_bezier(**kwargs):
+            # start and endpoints can be extrema of the curve
+            _update_extents(extents, curve.point_at_t(0))  # start point
+            _update_extents(extents, curve.point_at_t(1))  # end point
+            # interior extrema where d/ds B(s) == 0
+            _, dzeros = curve.axis_aligned_extrema
+            if len(dzeros) == 0:
+                continue
+            for ti in dzeros:
+                potential_extrema = curve.point_at_t(ti)
+                _update_extents(extents, potential_extrema)
+        return Bbox.from_extents(extents)
 
     def intersects_path(self, other, filled=True):
         """
diff --git a/lib/matplotlib/tests/test_path.py b/lib/matplotlib/tests/test_path.py
@@ -49,6 +49,37 @@ def test_contains_points_negative_radius():
     np.testing.assert_equal(result, [True, False, False])
 
 
+_test_paths = [
+    # interior extrema determine extents and degenerate derivative
+    Path([[0, 0], [1, 0], [1, 1], [0, 1]],
+           [Path.MOVETO, Path.CURVE4, Path.CURVE4, Path.CURVE4]),
+    # a quadratic curve
+    Path([[0, 0], [0, 1], [1, 0]], [Path.MOVETO, Path.CURVE3, Path.CURVE3]),
+    # a linear curve, degenerate vertically
+    Path([[0, 1], [1, 1]], [Path.MOVETO, Path.LINETO]),
+    # a point
+    Path([[1, 2]], [Path.MOVETO]),
+]
+
+
+_test_path_extents = [(0., 0., 0.75, 1.), (0., 0., 1., 0.5), (0., 1., 1., 1.),
+                      (1., 2., 1., 2.)]
+
+
+@pytest.mark.parametrize('path, extents', zip(_test_paths, _test_path_extents))
+def test_exact_extents(path, extents):
+    # notice that if we just looked at the control points to get the bounding
+    # box of each curve, we would get the wrong answers. For example, for
+    # hard_curve = Path([[0, 0], [1, 0], [1, 1], [0, 1]],
+    #                   [Path.MOVETO, Path.CURVE4, Path.CURVE4, Path.CURVE4])
+    # we would get that the extents area (0, 0, 1, 1). This code takes into
+    # account the curved part of the path, which does not typically extend all
+    # the way out to the control points.
+    # Note that counterintuitively, path.get_extents() returns a Bbox, so we
+    # have to get that Bbox's `.extents`.
+    assert np.all(path.get_extents().extents == extents)
+
+
 def test_point_in_path_nan():
     box = np.array([[0, 0], [1, 0], [1, 1], [0, 1], [0, 0]])
     p = Path(box)
