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

Commit7054e98

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

File tree

5 files changed

+234
-24
lines changed

5 files changed

+234
-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: 109 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,114 @@ 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+
243+
For a bezier curve :math:`\vec{B}(t)`, to calculate the area of the arc
244+
swept out by the ray from the origin to the curve, we need to compute
245+
:math:`\frac{1}{2}\int_0^1 \vec{B}(t) \cdot \vec{n}(t) dt`, where
246+
:math:`\vec{n}(t) = u^{(1)}(t)\hat{x}^{(0)} - u^{(0)}(t)\hat{x}^{(1)}`
247+
is the normal vector oriented away from the origin and
248+
:math:`u^{(i)}(t) = \frac{d}{dt} B^{(i)}(t)` is the :math:`i`\th
249+
component of the curve's tangent vector. (This formula can be found by
250+
applying the divergence theorem to :math:`F(x,y) = [x, y]/2`, and
251+
calculates the *signed* area for a counter-clockwise curve, by the
252+
right hand rule).
253+
254+
The control points of the curve are its coefficients in a Bernstein
255+
expansion, so if we let :math:`P_i = [P^{(0)}_i, P^{(1)}_i]` be the
256+
:math:`i`\th control point, then
257+
258+
.. math::
259+
260+
\frac{1}{2}\int_0^1 B(t) \cdot n(t) dt
261+
&= \frac{1}{2}\int_0^1 B^{(0)}(t) \frac{d}{dt} B^{(1)}(t)
262+
- B^{(1)}(t) \frac{d}{dt} B^{(0)}(t)
263+
dt \\
264+
&= \frac{1}{2}\int_0^1
265+
\left( \sum_{j=0}^n P_j^{(0)} b_{j,n} \right)
266+
\left( n \sum_{k=0}^{n-1} (P_{k+1}^{(1)} -
267+
P_{k}^{(1)}) b_{j,n} \right)
268+
\\
269+
&\hspace{1em} - \left( \sum_{j=0}^n P_j^{(1)} b_{j,n}
270+
\right) \left( n \sum_{k=0}^{n-1} (P_{k+1}^{(0)}
271+
- P_{k}^{(0)}) b_{j,n} \right)
272+
dt,
273+
274+
where :math:`b_{\nu, n}(t) = {n \choose \nu} t^\nu {(1 - t)}^{n-\nu}`
275+
is the :math:`\nu`\th Bernstein polynomial of degree :math:`n`.
276+
277+
Grouping :math:`t^l(1-t)^m` terms together for each :math:`l`,
278+
:math:`m`, we get that the integrand becomes
279+
280+
.. math::
281+
282+
\sum_{j=0}^n \sum_{k=0}^{n-1}
283+
{n \choose j} {{n - 1} \choose k}
284+
&\left[P_j^{(0)} (P_{k+1}^{(1)} - P_{k}^{(1)})
285+
- P_j^{(1)} (P_{k+1}^{(0)} - P_{k}^{(0)})\right] \\
286+
&\hspace{1em}\times{}t^{j + k} {(1 - t)}^{2n - 1 - j - k}
287+
288+
or more compactly,
289+
290+
.. math::
291+
292+
\sum_{j=0}^n \sum_{k=0}^{n-1}
293+
\frac{{n \choose j} {{n - 1} \choose k}}
294+
{{{2n - 1} \choose {j+k}}}
295+
[P_j^{(0)} (P_{k+1}^{(1)} - P_{k}^{(1)})
296+
- P_j^{(1)} (P_{k+1}^{(0)} - P_{k}^{(0)})]
297+
b_{j+k,2n-1}(t).
298+
299+
Interchanging sum and integral, and using the fact that :math:`\int_0^1
300+
b_{\nu, n}(t) dt = \frac{1}{n + 1}`, we conclude that the original
301+
integral can be written as
302+
303+
.. math::
304+
305+
\frac{1}{2}&\int_0^1 B(t) \cdot n(t) dt
306+
\\
307+
&= \frac{1}{4}\sum_{j=0}^n \sum_{k=0}^{n-1}
308+
\frac{{n \choose j} {{n - 1} \choose k}}
309+
{{{2n - 1} \choose {j+k}}}
310+
[P_j^{(0)} (P_{k+1}^{(1)} - P_{k}^{(1)})
311+
- P_j^{(1)} (P_{k+1}^{(0)} - P_{k}^{(0)})]
312+
"""
313+
n=self.degree
314+
P=self.control_points
315+
dP=np.diff(P,axis=0)
316+
j=np.arange(n+1)
317+
k=np.arange(n)
318+
return (1/4)*np.sum(
319+
np.multiply.outer(_comb(n,j),_comb(n-1,k))
320+
/_comb(2*n-1,np.add.outer(j,k))
321+
* (np.multiply.outer(P[j,0],dP[k,1])-
322+
np.multiply.outer(P[j,1],dP[k,0]))
323+
)
324+
325+
@classmethod
326+
defdifferentiate(cls,B):
327+
"""Return the derivative of a BezierSegment, itself a BezierSegment"""
328+
dcontrol_points=B.degree*np.diff(B.control_points,axis=0)
329+
returncls(dcontrol_points)
330+
232331
@property
233332
defcontrol_points(self):
234333
"""The control points of the curve."""

‎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/tests/test_bezier.py

Lines changed: 24 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,23 @@ 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+
returnnp.cross(B(t),dB(t))/2
34+
x=np.linspace(0,1,1000)
35+
y=integrand(x)
36+
returnnp.trapz(y,x)
37+
38+
39+
@pytest.mark.parametrize("B",_test_curves)
40+
deftest_area_formula(B):
41+
assertnp.isclose(_integral_arc_area(B),B.arc_area)

‎lib/matplotlib/tests/test_path.py

Lines changed: 39 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,27 @@ 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.),
101107
]
102108

103109

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):
110+
@pytest.mark.parametrize('precomputed_curve',_test_curves)
111+
deftest_exact_extents(precomputed_curve):
110112
# notice that if we just looked at the control points to get the bounding
111113
# box of each curve, we would get the wrong answers. For example, for
112114
# hard_curve = Path([[0, 0], [1, 0], [1, 1], [0, 1]],
@@ -116,6 +118,7 @@ def test_exact_extents(path, extents):
116118
# the way out to the control points.
117119
# Note that counterintuitively, path.get_extents() returns a Bbox, so we
118120
# have to get that Bbox's `.extents`.
121+
path,extents=precomputed_curve.path,precomputed_curve.extents
119122
assertnp.all(path.get_extents().extents==extents)
120123

121124

@@ -129,6 +132,28 @@ def test_extents_with_ignored_codes(ignored_code):
129132
assertnp.all(path.get_extents().extents== (0.,0.,1.,1.))
130133

131134

135+
@pytest.mark.parametrize('precomputed_curve',_test_curves)
136+
deftest_signed_area(precomputed_curve):
137+
path,area=precomputed_curve.path,precomputed_curve.area
138+
assertnp.isclose(path.signed_area(),area)
139+
# now flip direction, sign of *signed_area* should flip
140+
rcurve=Path(path.vertices[::-1],path.codes)
141+
assertnp.isclose(rcurve.signed_area(),-area)
142+
143+
144+
deftest_signed_area_unit_rectangle():
145+
rect=Path.unit_rectangle()
146+
assertnp.isclose(rect.signed_area(),1)
147+
148+
149+
deftest_signed_area_unit_circle():
150+
circ=Path.unit_circle()
151+
# Not a "real" circle, just an approximation of a circle made out of bezier
152+
# curves. The actual value is 3.1415935732517166, which is close enough to
153+
# pass here.
154+
assertnp.isclose(circ.signed_area(),np.pi)
155+
156+
132157
deftest_point_in_path_nan():
133158
box=np.array([[0,0], [1,0], [1,1], [0,1], [0,0]])
134159
p=Path(box)

0 commit comments

Comments
 (0)
[8]ページ先頭
©2009-2025 Movatter.jp