correctly compute bounding box for path

brunobeltran · brunobeltran · commit e0d3da149d98 · 2020-04-15T02:10:09.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,27 @@
+
+Functions to compute a Path's size
+----------------------------------
+
+Various functions were added to `~.bezier.BezierSegment` and `~.path.Path` to
+allow computation of the shape/size of a `~.path.Path` and its composite Bezier
+curves.
+
+In addition to the fixes below, `~.bezier.BezierSegment` has gained more
+documentation and usability improvements, including properties that contain its
+dimension, degree, control_points, and more.
+
+Better interface for Path segment iteration
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+`~.path.Path.iter_bezier` iterates through the `~.bezier.BezierSegment`'s that
+make up the Path. This is much more useful typically than the existing
+`~.path.Path.iter_segments` function, which returns the absolute minimum amount
+of information possible to reconstruct the Path.
+
+Fixed bug that computed a Path's Bbox incorrectly
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+Historically, `~.path.Path.get_extents` has always simply returned the Bbox of
+a curve's control points, instead of the Bbox of the curve itself. 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.
diff --git a/lib/matplotlib/bezier.py b/lib/matplotlib/bezier.py
@@ -2,12 +2,24 @@
 A module providing some utility functions regarding Bezier path manipulation.
 """
 
+from functools import lru_cache
 import math
+import warnings
 
 import numpy as np
 
 import matplotlib.cbook as cbook
 
+# same algorithm as 3.8's math.comb
+@np.vectorize
+@lru_cache(maxsize=128)
+def _comb(n, k):
+    if k > n:
+        return 0
+    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,26 +180,127 @@ 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).T
+
+    def __call__(self, t):
+        """
+        Evaluate the Bezier curve at point(s) t in [0, 1].
+
+        Parameters
+        ----------
+        t : float (k,), array_like
+            Points at which to evaluate the curve.
+
+        Returns
+        -------
+        float (k, d), array_like
+            Value of the curve for each point in *t*.
+        """
+        t = np.asarray(t)
+        return (np.power.outer(1 - t, self._orders[::-1])
+                * np.power.outer(t, self._orders)) @ self._px
 
     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))
+        """Evaluate curve at a single point *t*. Returns a Tuple[float*d]."""
+        return tuple(self(t))
+
+    @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.
+
+        .. warning:: Follows opposite convention from `numpy.polyval`.
+
+        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)
+        P = self.control_points
+        j = np.arange(n+1)[:, None]
+        i = np.arange(n+1)[None, :]  # _comb is non-zero for i <= j
+        prefactor = (-1)**(i + j) * _comb(j, i)  # j on axis 0, i on axis 1
+        return _comb(n, j) * prefactor @ P  # j on axis 0, self.dimension on 1
+
+    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(
diff --git a/lib/matplotlib/path.py b/lib/matplotlib/path.py
@@ -17,6 +17,7 @@
 import matplotlib as mpl
 from . import _path, cbook
 from .cbook import _to_unmasked_float_array, simple_linear_interpolation
+from .bezier import BezierSegment
 
 
 class Path:
@@ -420,6 +421,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 +576,32 @@ 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)
+        bbox = Bbox.null()
+        for curve, code in self.iter_bezier(**kwargs):
+            # places where the derivative is zero can be extrema
+            _, dzeros = curve.axis_aligned_extrema()
+            # as can the ends of the curve
+            bbox.update_from_data_xy(curve([0, *dzeros, 1]), ignore=False)
+        return bbox
 
     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)
diff --git a/lib/matplotlib/transforms.py b/lib/matplotlib/transforms.py
@@ -845,8 +845,8 @@ def ignore(self, value):
 
     def update_from_path(self, path, ignore=None, updatex=True, updatey=True):
         """
-        Update the bounds of the `Bbox` based on the passed in
-        data.  After updating, the bounds will have positive *width*
+        Update the bounds of the `Bbox` to contain the vertices of the
+        provided path. After updating, the bounds will have positive *width*
         and *height*; *x0* and *y0* will be the minimal values.
 
         Parameters
