From eb92405d21db3e81bf5148356b0c9a014c28b7c1 Mon Sep 17 00:00:00 2001 From: Bruno Beltran Date: Fri, 20 Mar 2020 01:59:28 -0700 Subject: [PATCH] add function to compute (signed) area of path --- .../next_whats_new/path-size-methods.rst | 6 + lib/matplotlib/bezier.py | 126 ++++++++++++++++-- lib/matplotlib/bezier.pyi | 4 + lib/matplotlib/path.py | 56 ++++++++ lib/matplotlib/path.pyi | 1 + lib/matplotlib/tests/test_bezier.py | 29 ++++ lib/matplotlib/tests/test_path.py | 55 ++++++-- 7 files changed, 253 insertions(+), 24 deletions(-) create mode 100644 doc/users/next_whats_new/path-size-methods.rst diff --git a/doc/users/next_whats_new/path-size-methods.rst b/doc/users/next_whats_new/path-size-methods.rst new file mode 100644 index 000000000000..1ab7c7452586 --- /dev/null +++ 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 index 42a6b478d729..bdba09bf1e00 100644 --- a/lib/matplotlib/bezier.py +++ b/lib/matplotlib/bezier.py @@ -2,7 +2,6 @@ A module providing some utility functions regarding Bézier path manipulation. """ -from functools import lru_cache import math import warnings @@ -11,15 +10,7 @@ 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 +220,121 @@ def point_at_t(self, t): """ return tuple(self(t)) + @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 index ad82b873affd..5900287a318e 100644 --- a/lib/matplotlib/bezier.pyi +++ b/lib/matplotlib/bezier.pyi @@ -36,6 +36,10 @@ class BezierSegment: def __call__(self, t: ArrayLike) -> np.ndarray: ... def point_at_t(self, t: float) -> tuple[float, ...]: ... @property + def arc_area(self) -> float: ... + @classmethod + def differentiate(cls, B: BezierSegment) -> BezierSegment: ... + @property def control_points(self) -> np.ndarray: ... @property def dimension(self) -> int: ... diff --git a/lib/matplotlib/path.py b/lib/matplotlib/path.py index 5f5a0f3de423..ee5757677259 100644 --- a/lib/matplotlib/path.py +++ b/lib/matplotlib/path.py @@ -666,6 +666,62 @@ def intersects_bbox(self, bbox, filled=True): return _path.path_intersects_rectangle( self, bbox.x0, bbox.y0, bbox.x1, bbox.y1, filled) + 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 index 464fc6d9a912..3bd0689e79a2 100644 --- a/lib/matplotlib/path.pyi +++ b/lib/matplotlib/path.pyi @@ -90,6 +90,7 @@ 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 signed_area(self) -> float: ... def interpolated(self, steps: int) -> Path: ... def to_polygons( self, diff --git a/lib/matplotlib/tests/test_bezier.py b/lib/matplotlib/tests/test_bezier.py index 65e2c616e738..0cbef9578962 100644 --- a/lib/matplotlib/tests/test_bezier.py +++ b/lib/matplotlib/tests/test_bezier.py @@ -3,6 +3,10 @@ """ from matplotlib.bezier import inside_circle, split_bezier_intersecting_with_closedpath +from matplotlib.tests.test_path import _test_curves + +import numpy as np +import pytest def test_split_bezier_with_large_values(): @@ -15,3 +19,28 @@ def test_split_bezier_with_large_values(): # All we are testing is that this completes # The failure case is an infinite loop resulting from floating point precision # pytest will timeout if that occurs + + +# get several curves to test our code on by borrowing the tests cases used in +# `~.tests.test_path`. get last path element ([-1]) and curve, not code ([0]) +_test_curves = [list(tc.path.iter_bezier())[-1][0] for tc in _test_curves] +# np2+ uses trapezoid, but we need to fallback to trapz for np<2 since it isn't there +_trapezoid = getattr(np, "trapezoid", np.trapz) # type: ignore[attr-defined] + + +def _integral_arc_area(B): + """(Signed) area swept out by ray from origin to curve.""" + dB = B.differentiate(B) + def integrand(t): + x = B(t) + y = dB(t) + # np.cross for 2d input + return (x[:, 0] * y[:, 1] - x[:, 1] * y[:, 0]) / 2 + x = np.linspace(0, 1, 1000) + y = integrand(x) + return _trapezoid(y, x) + + +@pytest.mark.parametrize("B", _test_curves) +def test_area_formula(B): + assert np.isclose(_integral_arc_area(B), B.arc_area) diff --git a/lib/matplotlib/tests/test_path.py b/lib/matplotlib/tests/test_path.py index 2c4df6ea3b39..305007665046 100644 --- a/lib/matplotlib/tests/test_path.py +++ b/lib/matplotlib/tests/test_path.py @@ -1,8 +1,8 @@ import platform import re +from collections import namedtuple import numpy as np - from numpy.testing import assert_array_equal import pytest @@ -88,25 +88,29 @@ def test_contains_points_negative_radius(): np.testing.assert_equal(result, [True, False, False]) -_test_paths = [ +_ExampleCurve = namedtuple('_ExampleCurve', ['path', 'extents', 'area']) +_test_curves = [ # 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]), + _ExampleCurve(Path([[0, 0], [1, 0], [1, 1], [0, 1]], + [Path.MOVETO, Path.CURVE4, Path.CURVE4, Path.CURVE4]), + extents=(0., 0., 0.75, 1.), area=0.6), + # a quadratic curve, clockwise + _ExampleCurve(Path([[0, 0], [0, 1], [1, 0]], + [Path.MOVETO, Path.CURVE3, Path.CURVE3]), + extents=(0., 0., 1., 0.5), area=-1/3), # a linear curve, degenerate vertically - Path([[0, 1], [1, 1]], [Path.MOVETO, Path.LINETO]), + _ExampleCurve(Path([[0, 1], [1, 1]], [Path.MOVETO, Path.LINETO]), + extents=(0., 1., 1., 1.), area=0.), # a point - Path([[1, 2]], [Path.MOVETO]), + _ExampleCurve(Path([[1, 2]], [Path.MOVETO]), extents=(1., 2., 1., 2.), + area=0.), + # non-curved triangle + _ExampleCurve(Path([(1, 1), (2, 1), (1.5, 2)]), extents=(1, 1, 2, 2), area=0.5), ] -_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): +@pytest.mark.parametrize('precomputed_curve', _test_curves) +def test_exact_extents(precomputed_curve): # 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]], @@ -116,6 +120,7 @@ def test_exact_extents(path, extents): # 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`. + path, extents = precomputed_curve.path, precomputed_curve.extents assert np.all(path.get_extents().extents == extents) @@ -129,6 +134,28 @@ def test_extents_with_ignored_codes(ignored_code): assert np.all(path.get_extents().extents == (0., 0., 1., 1.)) +@pytest.mark.parametrize('precomputed_curve', _test_curves) +def test_signed_area(precomputed_curve): + path, area = precomputed_curve.path, precomputed_curve.area + np.testing.assert_allclose(path.signed_area(), area) + # now flip direction, sign of *signed_area* should flip + rcurve = Path(path.vertices[::-1], path.codes) + np.testing.assert_allclose(rcurve.signed_area(), -area) + + +def test_signed_area_unit_rectangle(): + rect = Path.unit_rectangle() + assert rect.signed_area() == 1 + + +def test_signed_area_unit_circle(): + circ = Path.unit_circle() + # Not a "real" circle, just an approximation of a circle made out of bezier + # curves. The actual value is 3.1415935732517166, which is close enough to + # pass here. + assert np.isclose(circ.signed_area(), np.pi) + + def test_point_in_path_nan(): box = np.array([[0, 0], [1, 0], [1, 1], [0, 1], [0, 0]]) p = Path(box)