matplotlib
diff --git a/‎lib/matplotlib/bezier.py
Lines changed: 71 additions & 8 deletions b/‎lib/matplotlib/bezier.py
Lines changed: 71 additions & 8 deletions
diff --git a/‎lib/matplotlib/tests/test_path.py
Lines changed: 16 additions & 0 deletions b/‎lib/matplotlib/tests/test_path.py
Lines changed: 16 additions & 0 deletions
@@ -6,7 +6,6 @@
 import warnings
 
 import numpy as np
-import scipy.integrate
 
 import matplotlib.cbook as cbook
 
@@ -256,13 +255,77 @@ def point_at_t(self, t):
         return tuple(self(t))
 
     def arc_area(self):
-        """(Signed) area between curve an origin."""
-        dB = self.differentiate(self)
-        def integrand(t):
-            dr = dB(t).T
-            n = np.array([dr[1], -dr[0]])
-            return (self(t).T @ n) / 2
-        return scipy.integrate.quad(integrand, 0, 1)[0]
+        r"""(Signed) area swept out by ray from origin to curve.
+
+        Notes
+        -----
+        Exact formula used for Bezier curves.
+
+        Now for an arbitrary bezier curve B(t), the area of the arc swept
+        out by the ray from the origin to the curve, we simply need to compute
+        :math:`\frac{1}{2}\int_0^1 B(t) \cdot n(t) dt`, where :math:`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`.
+
+        The control points of the curve are just its coefficients in a
+        Bernstein expansion, so if we let :math:`P_i = [x_i, y_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)
+                      - \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)` is the :math:`\nu`'th Bernstein polynomial
+        of degree :math:`n`, :math:`b_{\nu, n}(t) = {n \choose \nu} t^\nu {(1 -
+        t)}^{n-\nu}`.
+
+        grouping :math:`t^l(1-t)^m` terms together for each :math:`l`,
+        :math:`m`, we get that the integrand becomes
+
+        .. math::
+
+            \frac{n}{2}\int_0^1 \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) dt
+
+        Interchanging sum and integral, we notice that :math:`\int_0^1 b_{\nu,
+        n}(t) dt = \frac{1}{n + 1}`, so we conclude that
+
+        .. 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)})]
+        """
+        n = self.degree
+        area = 0
+        P = self.control_points
+        dP = np.diff(P, axis=0)
+        for j in range(n + 1):
+            for k in range(n):
+                area += _comb(n, j)*_comb(n-1, k)/_comb(2*n - 1, j + k) \
+                        * (P[j, 0]*dP[k, 1] - P[j, 1]*dP[k, 0])
+        return (1/4)*area
 
     @classmethod
     def differentiate(cls, B):
 
@@ -55,6 +55,22 @@ def test_exact_extents_cubic():
     np.testing.assert_equal(hard_curve.get_exact_extents(), [0., 0., 0.75, 1.])
 
 
+def test_signed_area_unit_circle():
+    circ = Path.unit_circle()
+    # not quite pi...since it's not quite a circle!
+    assert(np.isclose(circ.signed_area(), 3.1415935732517166))
+    # now counter-clockwise
+    rverts = circ.vertices[-2::-1]
+    rverts = np.append(rverts, np.atleast_2d(circ.vertices[0]), axis=0)
+    rcirc = Path(rverts, circ.codes)
+    assert(np.isclose(rcirc.signed_area(), -3.1415935732517166))
+
+
+def test_signed_area_unit_rectangle():
+    rect = Path.unit_rectangle()
+    assert(np.isclose(rect.signed_area(), 1))
+
+
 def test_point_in_path_nan():
     box = np.array([[0, 0], [1, 0], [1, 1], [0, 1], [0, 0]])
     p = Path(box)