10000 add function to compute (signed) area of path · matplotlib/matplotlib@eb92405 · GitHub
[go: up one dir, main page]
Skip to content

Commit eb92405

Browse files
brunobeltrangreglucas
brunobeltran
authored and
greglucas
committed
add function to compute (signed) area of path
1 parent bda958a commit eb92405

File tree

7 files changed

+253
-24
lines changed

7 files changed

+253
-24
lines changed

‎doc/users/next_whats_new/path-size-methods.rst

Lines changed: 6 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,6 @@
1+
Path area
2+
~~~~~~~~~
3+
4+
A `~.path.Path.signed_area` method was added to compute the signed filled area
5+
of a Path object analytically (i.e. without integration). This should be useful
6+
for constructing Paths of a desired area.

‎lib/matplotlib/bezier.py

Lines changed: 116 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -2,7 +2,6 @@
22
A module providing some utility functions regarding Bézier path manipulation.
33
"""
44

5-
from functools import lru_cache
65
import math
76
import warnings
87

@@ -11,15 +10,7 @@
1110
from matplotlib import _api
1211

1312

14-
# same algorithm as 3.8's math.comb
15-
@np.vectorize
16-
@lru_cache(maxsize=128)
17-
def _comb(n, k):
18-
if k > n:
19-
return 0
20-
k = min(k, n - k)
21-
i = np.arange(1, k + 1)
22-
return np.prod((n + 1 - i)/i).astype(int)
13+
_comb = np.vectorize(math.comb, otypes=[int])
2314

2415

2516
class NonIntersectingPathException(ValueError):
@@ -229,6 +220,121 @@ def point_at_t(self, t):
229220
"""
230221
return tuple(self(t))
231222

223+
@property
224+
def arc_area(self):
225+
r"""
226+
Signed area swept out by ray from origin to curve.
227+
228+
Counterclockwise area is counted as positive, and clockwise area as
229+
negative.
230+
231+
The sum of this function for each Bezier curve in a Path will give the
232+
signed area enclosed by the Path.
233+
234+
Returns
235+
-------
236+
float
237+
The signed area of the arc swept out by the curve.
238+
239+
Notes
240+
-----
241+
An analytical formula is possible for arbitrary bezier curves.
242+
The formulas can be found in computer graphics references [1]_ and
243+
an example derivation is given below.
244+
245+
For a bezier curve :math:`\vec{B}(t)`, to calculate the area of the arc
246+
swept out by the ray from the origin to the curve, we need to compute
247+
:math:`\frac{1}{2}\int_0^1 \vec{B}(t) \cdot \vec{n}(t) dt`, where
248+
:math:`\vec{n}(t) = u^{(1)}(t)\hat{x}^{(0)} - u^{(0)}(t)\hat{x}^{(1)}`
249+
is the normal vector oriented away from the origin and
250+
:math:`u^{(i)}(t) = \frac{d}{dt} B^{(i)}(t)` is the :math:`i`\th
251+
component of the curve's tangent vector. (This formula can be found by
252+
applying the divergence theorem to :math:`F(x,y) = [x, y]/2`, and
253+
calculates the *signed* area for a counter-clockwise curve, by the
254+
right hand rule).
255+
256+
The control points of the curve are its coefficients in a Bernstein
257+
expansion, so if we let :math:`P_i = [P^{(0)}_i, P^{(1)}_i]` be the
258+
:math:`i`\th control point, then
259+
260+
.. math::
261+
262+
\frac{1}{2}\int_0^1 B(t) \cdot n(t) dt
263+
&= \frac{1}{2}\int_0^1 B^{(0)}(t) \frac{d}{dt} B^{(1)}(t)
264+
- B^{(1)}(t) \frac{d}{dt} B^{(0)}(t)
265+
dt \\
266+
&= \frac{1}{2}\int_0^1
267+
\left( \sum_{j=0}^n P_j^{(0)} b_{j,n} \right)
268+
\left( n \sum_{k=0}^{n-1} (P_{k+1}^{(1)} -
269+
P_{k}^{(1)}) b_{j,n} \right)
270+
\\
271+
&\hspace{1em} - \left( \sum_{j=0}^n P_j^{(1)} b_{j,n}
272+
\right) \left( n \sum_{k=0}^{n-1} (P_{k+1}^{(0)}
273+
- P_{k}^{(0)}) b_{j,n} \right)
274+
dt,
275+
276+
where :math:`b_{\nu, n}(t) = {n \choose \nu} t^\nu {(1 - t)}^{n-\nu}`
277+
is the :math:`\nu`\th Bernstein polynomial of degree :math:`n`.
278+
279+
Grouping :math:`t^l(1-t)^m` terms together for each :math:`l`,
280+
:math:`m`, we get that the integrand becomes
281+
282+
.. math::
283+
284+
\sum_{j=0}^n \sum_{k=0}^{n-1}
285+
{n \choose j} {{n - 1} \choose k}
286+
&\left[P_j^{(0)} (P_{k+1}^{(1)} - P_{k}^{(1)})
287+
- P_j^{(1)} (P_{k+1}^{(0)} - P_{k}^{(0)})\right] \\
288+
&\hspace{1em}\times{}t^{j + k} {(1 - t)}^{2n - 1 - j - k}
289+
290+
or more compactly,
291+
292+
.. math::
293+
294+
\sum_{j=0}^n \sum_{k=0}^{n-1}
295+
\frac{{n \choose j} {{n - 1} \choose k}}
296+
{{{2n - 1} \choose {j+k}}}
297+
[P_j^{(0)} (P_{k+1}^{(1)} - P_{k}^{(1)})
298+
- P_j^{(1)} (P_{k+1}^{(0)} - P_{k}^{(0)})]
299+
b_{j+k,2n-1}(t).
300+
301+
Interchanging sum and integral, and using the fact that :math:`\int_0^1
302+
b_{\nu, n}(t) dt = \frac{1}{n + 1}`, we conclude that the original
303+
integral can be written as
304+
305+
.. math::
306+
307+
\frac{1}{2}&\int_0^1 B(t) \cdot n(t) dt
308+
\\
309+
&= \frac{1}{4}\sum_{j=0}^n \sum_{k=0}^{n-1}
310+
\frac{{n \choose j} {{n - 1} \choose k}}
311+
{{{2n - 1} \choose {j+k}}}
312+
[P_j^{(0)} (P_{k+1}^{(1)} - P_{k}^{(1)})
313+
- P_j^{(1)} (P_{k+1}^{(0)} - P_{k}^{(0)})]
314+
315+
References
316+
----------
317+
.. [1] Sederberg, Thomas W., "Computer Aided Geometric Design" (2012).
318+
Faculty Publications. 1. https://scholarsarchive.byu.edu/facpub/1
319+
"""
320+
n = self.degree
321+
P = self.control_points
322+
dP = np.diff(P, axis=0)
323+
j = np.arange(n + 1)
324+
k = np.arange(n)
325+
return (1/4)*np.sum(
326+
np.multiply.outer(_comb(n, j), _comb(n - 1, k))
327+
/ _comb(2*n - 1, np.add.outer(j, k))
328+
* (np.multiply.outer(P[j, 0], dP[k, 1]) -
329+
np.multiply.outer(P[j, 1], dP[k, 0]))
330+
)
331+
332+
@classmethod
333+
def differentiate(cls, B):
334+
"""Return the derivative of a BezierSegment, itself a BezierSegment"""
335+
dcontrol_points = B.degree*np.diff(B.control_points, axis=0)
336+
return cls(dcontrol_points)
337+
232338
@property
233339
def control_points(self):
234340
"""The control points of the curve."""

‎lib/matplotlib/bezier.pyi

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -36,6 +36,10 @@ class BezierSegment:
3636
def __call__(self, t: ArrayLike) -> np.ndarray: ...
3737
def point_at_t(self, t: float) -> tuple[float, ...]: ...
3838
@property
39+
def arc_area(self) -> float: ...
40+
@classmethod
41+
def differentiate(cls, B: BezierSegment) -> BezierSegment: ...
42+
@property
3943
def control_points(self) -> np.ndarray: ...
4044
@property
4145
def dimension(self) -> int: ...

‎lib/matplotlib/path.py

Lines changed: 56 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -666,6 +666,62 @@ def intersects_bbox(self, bbox, filled=True):
666666
return _path.path_intersects_rectangle(
667667
self, bbox.x0, bbox.y0, bbox.x1, bbox.y1, filled)
668668

669+
def signed_area(self):
670+
"""
671+
Get signed area of the filled path.
672+
673+
Area of a filled region is treated as positive if the path encloses it
674+
in a counter-clockwise direction, but negative if the path encloses it
675+
moving clockwise.
676+
677+
All sub paths are treated as if they had been closed. That is, if there
678+
is a MOVETO without a preceding CLOSEPOLY, one is added.
679+
680+
If the path is made up of multiple components that overlap, the
681+
overlapping area is multiply counted.
682+
683+
Returns
684+
-------
685+
float
686+
The signed area enclosed by the path.
687+
688+
Notes
689+
-----
690+
If the Path is not self-intersecting and has no overlapping components,
691+
then the absolute value of the signed area is equal to the actual
692+
filled area when the Path is drawn (e.g. as a PathPatch).
693+
694+
Examples
695+
--------
696+
A symmetric figure eight, (where one loop is clockwise and
697+
the other counterclockwise) would have a total *signed_area* of zero,
698+
since the two loops would cancel each other out.
699+
"""
700+
area = 0
701+
prev_point = None
702+
prev_code = None
703+
start_point = None
704+
for B, code in self.iter_bezier():
705+
# MOVETO signals the start of a new connected component of the path
706+
if code == Path.MOVETO:
707+
# if the previous segment exists and it doesn't close the
708+
# previous connected component of the path, do so manually to
709+
# match Agg's convention of filling unclosed path segments
710+
if prev_code not in (None, Path.CLOSEPOLY):
711+
Bclose = BezierSegment([prev_point, start_point])
712+
area += Bclose.arc_area
713+
# to allow us to manually close this connected component later
714+
start_point = B.control_points[0]
715+
area += B.arc_area
716+
prev_point = B.control_points[-1]
717+
prev_code = code
718+
# add final implied CLOSEPOLY, if necessary
719+
if start_point is not None \
720+
and not np.all(np.isclose(start_point, prev_point)):
721+
B = BezierSegment([prev_point, start_point])
722+
area += B.arc_area
723+
return area
724+
669725
def interpolated(self, steps):
670726
"""
671727
Return a new path resampled to length N x *steps*.

‎lib/matplotlib/path.pyi

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -90,6 +90,7 @@ class Path:
9090
def get_extents(self, transform: Transform | None = ..., **kwargs) -> Bbox: ...
9191
def intersects_path(self, other: Path, filled: bool = ...) -> bool: ...
9292
def intersects_bbox(self, bbox: Bbox, filled: bool = ...) -> bool: ...
93+
def signed_area(self) -> float: ...
9394
def interpolated(self, steps: int) -> Path: ...
9495
def to_polygons(
9596
self,

‎lib/matplotlib/tests/test_bezier.py

Lines changed: 29 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -3,6 +3,10 @@
33
"""
44

55
from matplotlib.bezier import inside_circle, split_bezier_intersecting_with_closedpath
6+
from matplotlib.tests.test_path import _test_curves
7+
8+
import numpy as np
9+
import pytest
610

711

812
def test_split_bezier_with_large_values():
@@ -15,3 +19,28 @@ def test_split_bezier_with_large_values():
1519
# All we are testing is that this completes
1620
# The failure case is an infinite loop resulting from floating point precision
1721
# pytest will timeout if that occurs
22+
23+
24+
# get several curves to test our code on by borrowing the tests cases used in
25+
# `~.tests.test_path`. get last path element ([-1]) and curve, not code ([0])
26+
_test_curves = [list(tc.path.iter_bezier())[-1][0] for tc in _test_curves]
27+
# np2+ uses trapezoid, but we need to fallback to trapz for np<2 since it isn't there
28+
_trapezoid = getattr(np, "trapezoid", np.trapz) # type: ignore[attr-defined]
29+
30+
31+
def _integral_arc_area(B):
32+
"""(Signed) area swept out by ray from origin to curve."""
33+
dB = B.differentiate(B)
34+
def integrand(t):
35+
x = B(t)
36+
y = dB(t)
37+
# np.cross for 2d input
38+
return (x[:, 0] * y[:, 1] - x[:, 1] * y[:, 0]) / 2
39+
x = np.linspace(0, 1, 1000)
40+
y = integrand(x)
41+
return _trapezoid(y, x)
42+
43+
44+
@pytest.mark.parametrize("B", _test_curves)
45+
def test_area_formula(B):
46+
assert np.isclose(_integral_arc_area(B), B.arc_area)

‎lib/matplotlib/tests/test_path.py

Lines changed: 41 additions & 14 deletions
Original file line numberDiff line numberDiff line change
@@ -1,8 +1,8 @@
11
import platform
22
import re
3+
from collections import namedtuple
34

45
import numpy as np
5-
66
from numpy.testing import assert_array_equal
77
import pytest
88

@@ -88,25 +88,29 @@ def test_contains_points_negative_radius():
8888
np.testing.assert_equal(result, [True, False, False])
8989

9090

91-
_test_paths = [
91+
_ExampleCurve = namedtuple('_ExampleCurve', ['path', 'extents', 'area'])
92+
_test_curves = [
9293
# interior extrema determine extents and degenerate derivative
93-
Path([[0, 0], [1, 0], [1, 1], [0, 1]],
94-
[Path.MOVETO, Path.CURVE4, Path.CURVE4, Path.CURVE4]),
95-
# a quadratic curve
96-
Path([[0, 0], [0, 1], [1, 0]], [Path.MOVETO, Path.CURVE3, Path.CURVE3]),
94+
_ExampleCurve(Path([[0, 0], [1, 0], [1, 1], [0, 1]],
95+
[Path.MOVETO, Path.CURVE4, Path.CURVE4, Path.CURVE4]),
96+
extents=(0., 0., 0.75, 1.), area=0.6),
97+
# a quadratic curve, clockwise
98+
_ExampleCurve(Path([[0, 0], [0, 1], [1, 0]],
99+
[Path.MOVETO, Path.CURVE3, Path.CURVE3]),
100+
extents=(0., 0., 1., 0.5), area=-1/3),
97101
# a linear curve, degenerate vertically
98-
Path([[0, 1], [1, 1]], [Path.MOVETO, Path.LINETO]),
102+
_ExampleCurve(Path([[0, 1], [1, 1]], [Path.MOVETO, Path.LINETO]),
103+
extents=(0., 1., 1., 1.), area=0.),
99104
# a point
100-
Path([[1, 2]], [Path.MOVETO]),
105+
_ExampleCurve(Path([[1, 2]], [Path.MOVETO]), extents=(1., 2., 1., 2.),
106+
area=0.),
107+
# non-curved triangle
108+
_ExampleCurve(Path([(1, 1), (2, 1), (1.5, 2)]), extents=(1, 1, 2, 2), area=0.5),
101109
]
102110

103111

104-
_test_path_extents = [(0., 0., 0.75, 1.), (0., 0., 1., 0.5), (0., 1., 1., 1.),
105-
(1., 2., 1., 2.)]
106-
107-
108-
@pytest.mark.parametrize('path, extents', zip(_test_paths, _test_path_extents))
109-
def test_exact_extents(path, extents):
112+
@pytest.mark.parametrize('precomputed_curve', _test_curves)
113+
def test_exact_extents(precomputed_curve):
110114
# notice that if we just looked at the control points to get the bounding
111115
# box of each curve, we would get the wrong answers. For example, for
112116
# hard_curve = Path([[0, 0], [1, 0], [1, 1], [0, 1]],
@@ -116,6 +120,7 @@ def test_exact_extents(path, extents):
116120
# the way out to the control points.
117121
# Note that counterintuitively, path.get_extents() returns a Bbox, so we
118122
# have to get that Bbox's `.extents`.
123+
path, extents = precomputed_curve.path, precomputed_curve.extents
119124
assert np.all(path.get_extents().extents == extents)
120125

121126

@@ -129,6 +134,28 @@ def test_extents_with_ignored_codes(ignored_code):
129134
assert np.all(path.get_extents().extents == (0., 0., 1., 1.))
130135

131136

137+
@pytest.mark.parametrize('precomputed_curve', _test_curves)
138+
def test_signed_area(precomputed_curve):
139+
path, area = precomputed_curve.path, precomputed_curve.area
140+
np.testing.assert_allclose(path.signed_area(), area)
141+
# now flip direction, sign of *signed_area* should flip
142+
rcurve = Path(path.vertices[::-1], path.codes)
143+
np.testing.assert_allclose(rcurve.signed_area(), -area)
144+
145+
146+
def test_signed_area_unit_rectangle():
147+
rect = Path.unit_rectangle()
148+
assert rect.signed_area() == 1
149+
150+
151+
def test_signed_area_unit_circle():
152+
circ = Path.unit_circle()
153+
# Not a "real" circle, just an approximation of a circle made out of bezier
154+
# curves. The actual value is 3.1415935732517166, which is close enough to
155+
# pass here.
156+
assert np.isclose(circ.signed_area(), np.pi)
157+
158+
132159
def test_point_in_path_nan():
133160
box = np.array([[0, 0], [1, 0], [1, 1], [0, 1], [0, 0]])
134161
p = Path(box)

0 commit comments

Comments
 (0)
0