matplotlib · brunobeltran · Mar 20, 2020 · Mar 23, 2020
diff --git a/doc/users/next_whats_new/path-size-methods.rst b/doc/users/next_whats_new/path-size-methods.rst
@@ -0,0 +1,6 @@
+Path area
+~~~~~~~~~
+
+A `~.path.Path.signed_area` method was added to compute the signed filled area
+of a Path object analytically (i.e. without integration). This should be useful
+for constructing Paths of a desired area.
diff --git a/lib/matplotlib/bezier.py b/lib/matplotlib/bezier.py
@@ -2,24 +2,16 @@
 A module providing some utility functions regarding Bézier path manipulation.
 """
 
-from functools import lru_cache
 import math
 import warnings
+from collections import deque
 
 import numpy as np
 
 from matplotlib import _api
 
 
-# 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)
+_comb = np.vectorize(math.comb, otypes=[int])
 
 
 class NonIntersectingPathException(ValueError):
@@ -229,6 +221,208 @@
         """
         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=1e-4, atol=1e-6):
+        """
+        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
+        estimate 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.
+        """
+        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
+
+    @property
+    def arc_area(self):
+        r"""
+        Signed area swept out by ray from origin to curve.
+
+        Counterclockwise area is counted as positive, and clockwise area as
+        negative.
+
+        The sum of this function for each Bezier curve in a Path will give the
+        signed area enclosed by the Path.
+
+        Returns
+        -------
+        float
+            The signed area of the arc swept out by the curve.
+
+        Notes
+        -----
+        An analytical formula is possible for arbitrary bezier curves.
+        The formulas can be found in computer graphics references [1]_ and
+        an example derivation is given below.
+
+        For a bezier curve :math:`\vec{B}(t)`, to calculate the area of the arc
+        swept out by the ray from the origin to the curve, we need to compute
+        :math:`\frac{1}{2}\int_0^1 \vec{B}(t) \cdot \vec{n}(t) dt`, where
+        :math:`\vec{n}(t) = u^{(1)}(t)\hat{x}^{(0)} - u^{(0)}(t)\hat{x}^{(1)}`
+        is the normal vector oriented away from the origin and
+        :math:`u^{(i)}(t) = \frac{d}{dt} B^{(i)}(t)` is the :math:`i`\th
+        component of the curve's tangent vector.  (This formula can be found by
+        applying the divergence theorem to :math:`F(x,y) = [x, y]/2`, and
+        calculates the *signed* area for a counter-clockwise curve, by the
+        right hand rule).
+
+        The control points of the curve are its coefficients in a Bernstein
+        expansion, so if we let :math:`P_i = [P^{(0)}_i, P^{(1)}_i]` be the
+        :math:`i`\th control point, then
+
+        .. math::
+
+            \frac{1}{2}\int_0^1 B(t) \cdot n(t) dt
+                    &= \frac{1}{2}\int_0^1 B^{(0)}(t) \frac{d}{dt} B^{(1)}(t)
+                                        - B^{(1)}(t) \frac{d}{dt} B^{(0)}(t)
+                                        dt \\
+                    &= \frac{1}{2}\int_0^1
+                        \left( \sum_{j=0}^n P_j^{(0)} b_{j,n} \right)
+                        \left( n \sum_{k=0}^{n-1} (P_{k+1}^{(1)} -
+                                P_{k}^{(1)}) b_{j,n} \right)
+                        \\
+                    &\hspace{1em} - \left( \sum_{j=0}^n P_j^{(1)} b_{j,n}
+                        \right) \left( n \sum_{k=0}^{n-1} (P_{k+1}^{(0)}
+                                        - P_{k}^{(0)}) b_{j,n} \right)
+                        dt,
+
+        where :math:`b_{\nu, n}(t) = {n \choose \nu} t^\nu {(1 - t)}^{n-\nu}`
+        is the :math:`\nu`\th Bernstein polynomial of degree :math:`n`.
+
+        Grouping :math:`t^l(1-t)^m` terms together for each :math:`l`,
+        :math:`m`, we get that the integrand becomes
+
+        .. math::
+
+            \sum_{j=0}^n \sum_{k=0}^{n-1}
+                {n \choose j} {{n - 1} \choose k}
+                &\left[P_j^{(0)} (P_{k+1}^{(1)} - P_{k}^{(1)})
+                    - P_j^{(1)} (P_{k+1}^{(0)} - P_{k}^{(0)})\right] \\
+                &\hspace{1em}\times{}t^{j + k} {(1 - t)}^{2n - 1 - j - k}
+
+        or more compactly,
+
+        .. math::
+
+            \sum_{j=0}^n \sum_{k=0}^{n-1}
+                \frac{{n \choose j} {{n - 1} \choose k}}
+                        {{{2n - 1} \choose {j+k}}}
+                [P_j^{(0)} (P_{k+1}^{(1)} - P_{k}^{(1)})
+                    - P_j^{(1)} (P_{k+1}^{(0)} - P_{k}^{(0)})]
+                b_{j+k,2n-1}(t).
+
+        Interchanging sum and integral, and using the fact that :math:`\int_0^1
+        b_{\nu, n}(t) dt = \frac{1}{n + 1}`, we conclude that the original
+        integral can be written as
+
+        .. math::
+
+            \frac{1}{2}&\int_0^1 B(t) \cdot n(t) dt
+            \\
+            &= \frac{1}{4}\sum_{j=0}^n \sum_{k=0}^{n-1}
+                \frac{{n \choose j} {{n - 1} \choose k}}
+                    {{{2n - 1} \choose {j+k}}}
+                [P_j^{(0)} (P_{k+1}^{(1)} - P_{k}^{(1)})
+                - P_j^{(1)} (P_{k+1}^{(0)} - P_{k}^{(0)})]
+
+        References
+        ----------
+        .. [1] Sederberg, Thomas W., "Computer Aided Geometric Design" (2012).
+          Faculty Publications. 1. https://scholarsarchive.byu.edu/facpub/1
+        """
+        n = self.degree
+        P = self.control_points
+        dP = np.diff(P, axis=0)
+        j = np.arange(n + 1)
+        k = np.arange(n)
+        return (1/4)*np.sum(
+            np.multiply.outer(_comb(n, j), _comb(n - 1, k))
+            / _comb(2*n - 1, np.add.outer(j, k))
+            * (np.multiply.outer(P[j, 0], dP[k, 1]) -
+               np.multiply.outer(P[j, 1], dP[k, 0]))
+        )
+
+    @classmethod
+    def differentiate(cls, B):
+        """Return the derivative of a BezierSegment, itself a BezierSegment"""
+        dcontrol_points = B.degree*np.diff(B.control_points, axis=0)
+        return cls(dcontrol_points)
+
     @property
     def control_points(self):
         """The control points of the curve."""

diff --git a/lib/matplotlib/bezier.pyi b/lib/matplotlib/bezier.pyi
@@ -35,6 +35,13 @@ class BezierSegment:
     def __init__(self, control_points: ArrayLike) -> None: ...
     def __call__(self, t: ArrayLike) -> np.ndarray: ...
     def point_at_t(self, t: float) -> tuple[float, ...]: ...
+    def split_at_t(self, t: float) -> tuple[BezierSegment, BezierSegment]: ...
+    def control_net_length(self) -> float: ...
+    def arc_length(self, rtol: float=1e-4, atol: float=1e-6) -> float: ...
+    @property
+    def arc_area(self) -> float: ...
+    @classmethod
+    def differentiate(cls, B: BezierSegment) -> BezierSegment: ...
     @property
     def control_points(self) -> np.ndarray: ...
     @property

diff --git a/lib/matplotlib/path.py b/lib/matplotlib/path.py
@@ -666,6 +666,93 @@
         return _path.path_intersects_rectangle(
             self, bbox.x0, bbox.y0, bbox.x1, bbox.y1, filled)
 
+    def length(self, rtol=1e-4, atol=1e-6, **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.
+
+        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.
+
+        Returns
+        -------
+        length : float
+            The path length.
+        """
+        return np.sum([B.arc_length(rtol=rtol, atol=atol)
+                       for B, code in self.iter_bezier(**kwargs)])
+
+    def signed_area(self):
+        """
+        Get signed area of the filled path.
+
+        Area of a filled region is treated as positive if the path encloses it
+        in a counter-clockwise direction, but negative if the path encloses it
+        moving clockwise.
+
+        All sub paths are treated as if they had been closed. That is, if there
+        is a MOVETO without a preceding CLOSEPOLY, one is added.
+
+        If the path is made up of multiple components that overlap, the
+        overlapping area is multiply counted.
+
+        Returns
+        -------
+        float
+            The signed area enclosed by the path.
+
+        Notes
+        -----
+        If the Path is not self-intersecting and has no overlapping components,
+        then the absolute value of the signed area is equal to the actual
+        filled area when the Path is drawn (e.g. as a PathPatch).
+
+        Examples
+        --------
+        A symmetric figure eight, (where one loop is clockwise and
+        the other counterclockwise) would have a total *signed_area* of zero,
+        since the two loops would cancel each other out.
+        """
+        area = 0
+        prev_point = None
+        prev_code = None
+        start_point = None
+        for B, code in self.iter_bezier():
+            # MOVETO signals the start of a new connected component of the path
+            if code == Path.MOVETO:
+                # if the previous segment exists and it doesn't close the
+                # previous connected component of the path, do so manually to
+                # match Agg's convention of filling unclosed path segments
+                if prev_code not in (None, Path.CLOSEPOLY):
+                    Bclose = BezierSegment([prev_point, start_point])
+                    area += Bclose.arc_area
+                # to allow us to manually close this connected component later
+                start_point = B.control_points[0]
+            area += B.arc_area
+            prev_point = B.control_points[-1]
+            prev_code = code
+        # add final implied CLOSEPOLY, if necessary
+        if start_point is not None \
+                and not np.all(np.isclose(start_point, prev_point)):
+            B = BezierSegment([prev_point, start_point])
+            area += B.arc_area
+        return area
+
     def interpolated(self, steps):
         """
         Return a new path resampled to length N x *steps*.

diff --git a/lib/matplotlib/path.pyi b/lib/matplotlib/path.pyi
@@ -90,6 +90,8 @@ class Path:
     def get_extents(self, transform: Transform | None = ..., **kwargs) -> Bbox: ...
     def intersects_path(self, other: Path, filled: bool = ...) -> bool: ...
     def intersects_bbox(self, bbox: Bbox, filled: bool = ...) -> bool: ...
+    def length(self, rtol: float=1e-4, atol: float=1e-6, **kwargs) -> float: ...
+    def signed_area(self) -> float: ...
     def interpolated(self, steps: int) -> Path: ...
     def to_polygons(
         self,