Movatterモバイル変換

[0]ホーム
URL:

Skip to content

Navigation Menu

Sign in
Appearance settings

Search code, repositories, users, issues, pull requests...

Provide feedback

We read every piece of feedback, and take your input very seriously.

Saved searches

Use saved searches to filter your results more quickly

 Sign up
Appearance settings

Commitcdf05fe

Browse files
brunobeltrangreglucas
authored andcommitted
add function to compute (signed) area of path
1 parentbda958a commitcdf05fe

File tree

7 files changed

+251
-24
lines changed

7 files changed

+251
-24
lines changed
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-
fromfunctoolsimportlru_cache
65
importmath
76
importwarnings
87

@@ -11,15 +10,7 @@
1110
frommatplotlibimport_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-
ifk>n:
19-
return0
20-
k=min(k,n-k)
21-
i=np.arange(1,k+1)
22-
returnnp.prod((n+1-i)/i).astype(int)
13+
_comb=np.vectorize(math.comb,otypes=[int])
2314

2415

2516
classNonIntersectingPathException(ValueError):
@@ -229,6 +220,121 @@ def point_at_t(self, t):
229220
"""
230221
returntuple(self(t))
231222

223+
@property
224+
defarc_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+
defdifferentiate(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+
returncls(dcontrol_points)
337+
232338
@property
233339
defcontrol_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
defpoint_at_t(self,t:float)->tuple[float, ...]: ...
3838
@property
39+
defarc_area(self)->float: ...
40+
@classmethod
41+
defdifferentiate(cls,B:BezierSegment)->BezierSegment: ...
42+
@property
3943
defcontrol_points(self)->np.ndarray: ...
4044
@property
4145
defdimension(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+
defsigned_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+
forB,codeinself.iter_bezier():
705+
# MOVETO signals the start of a new connected component of the path
706+
ifcode==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+
ifprev_codenotin (None,Path.CLOSEPOLY):
711+
Bclose=BezierSegment(np.array([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+
ifstart_pointisnotNone \
720+
andnotnp.all(np.isclose(start_point,prev_point)):
721+
B=BezierSegment(np.array([prev_point,start_point]))
722+
area+=B.arc_area
723+
returnarea
724+
669725
definterpolated(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
defget_extents(self,transform:Transform|None= ...,**kwargs)->Bbox: ...
9191
defintersects_path(self,other:Path,filled:bool= ...)->bool: ...
9292
defintersects_bbox(self,bbox:Bbox,filled:bool= ...)->bool: ...
93+
defsigned_area(self)->float: ...
9394
definterpolated(self,steps:int)->Path: ...
9495
defto_polygons(
9596
self,

‎lib/matplotlib/tests/test_bezier.py

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

55
frommatplotlib.bezierimportinside_circle,split_bezier_intersecting_with_closedpath
6+
frommatplotlib.tests.test_pathimport_test_curves
7+
8+
importnumpyasnp
9+
importpytest
610

711

812
deftest_split_bezier_with_large_values():
@@ -15,3 +19,26 @@ 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]fortcin_test_curves]
27+
28+
29+
def_integral_arc_area(B):
30+
"""(Signed) area swept out by ray from origin to curve."""
31+
dB=B.differentiate(B)
32+
defintegrand(t):
33+
x=B(t)
34+
y=dB(t)
35+
# np.cross for 2d input
36+
return (x[:,0]*y[:,1]-x[:,1]*y[:,0])/2
37+
x=np.linspace(0,1,1000)
38+
y=integrand(x)
39+
returnnp.trapezoid(y,x)
40+
41+
42+
@pytest.mark.parametrize("B",_test_curves)
43+
deftest_area_formula(B):
44+
assertnp.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
importplatform
22
importre
3+
fromcollectionsimportnamedtuple
34

45
importnumpyasnp
5-
66
fromnumpy.testingimportassert_array_equal
77
importpytest
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-
deftest_exact_extents(path,extents):
112+
@pytest.mark.parametrize('precomputed_curve',_test_curves)
113+
deftest_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
assertnp.all(path.get_extents().extents==extents)
120125

121126

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

131136

137+
@pytest.mark.parametrize('precomputed_curve',_test_curves)
138+
deftest_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+
deftest_signed_area_unit_rectangle():
147+
rect=Path.unit_rectangle()
148+
assertrect.signed_area()==1
149+
150+
151+
deftest_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+
assertnp.isclose(circ.signed_area(),np.pi)
157+
158+
132159
deftest_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)
[8]ページ先頭
©2009-2025 Movatter.jp