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

Commit87e21f6

Browse files
committed
add function to compute (signed) area of path
1 parentd3eaffd commit87e21f6

File tree

6 files changed

+237
-15
lines changed

6 files changed

+237
-15
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 & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -12,14 +12,14 @@
1212

1313

1414
# same algorithm as 3.8's math.comb
15-
@np.vectorize
1615
@lru_cache(maxsize=128)
1716
def_comb(n,k):
1817
ifk>n:
1918
return0
2019
k=min(k,n-k)
2120
i=np.arange(1,k+1)
2221
returnnp.prod((n+1-i)/i).astype(int)
22+
_comb=np.vectorize(_comb,otypes=[np.int])
2323

2424

2525
classNonIntersectingPathException(ValueError):
@@ -220,6 +220,114 @@ def point_at_t(self, t):
220220
"""Evaluate curve at a single point *t*. Returns a Tuple[float*d]."""
221221
returntuple(self(t))
222222

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+
223331
@property
224332
defcontrol_points(self):
225333
"""The control points of the curve."""

‎lib/matplotlib/path.py‎

Lines changed: 56 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -626,6 +626,62 @@ def intersects_bbox(self, bbox, filled=True):
626626
return_path.path_intersects_rectangle(
627627
self,bbox.x0,bbox.y0,bbox.x1,bbox.y1,filled)
628628

629+
defsigned_area(self):
630+
"""
631+
Get signed area of the filled path.
632+
633+
Area of a filled region is treated as positive if the path encloses it
634+
in a counter-clockwise direction, but negative if the path encloses it
635+
moving clockwise.
636+
637+
All sub paths are treated as if they had been closed. That is, if there
638+
is a MOVETO without a preceding CLOSEPOLY, one is added.
639+
640+
If the path is made up of multiple components that overlap, the
641+
overlapping area is multiply counted.
642+
643+
Returns
644+
-------
645+
float
646+
The signed area enclosed by the path.
647+
648+
Notes
649+
-----
650+
If the Path is not self-intersecting and has no overlapping components,
651+
then the absolute value of the signed area is equal to the actual
652+
filled area when the Path is drawn (e.g. as a PathPatch).
653+
654+
Examples
655+
--------
656+
A symmetric figure eight, (where one loop is clockwise and
657+
the other counterclockwise) would have a total *signed_area* of zero,
658+
since the two loops would cancel each other out.
659+
"""
660+
area=0
661+
prev_point=None
662+
prev_code=None
663+
start_point=None
664+
forB,codeinself.iter_bezier():
665+
# MOVETO signals the start of a new connected component of the path
666+
ifcode==Path.MOVETO:
667+
# if the previous segment exists and it doesn't close the
668+
# previous connected component of the path, do so manually to
669+
# match Agg's convention of filling unclosed path segments
670+
ifprev_codenotin (None,Path.CLOSEPOLY):
671+
Bclose=BezierSegment(np.array([prev_point,start_point]))
672+
area+=Bclose.arc_area
673+
# to allow us to manually close this connected component later
674+
start_point=B.control_points[0]
675+
area+=B.arc_area
676+
prev_point=B.control_points[-1]
677+
prev_code=code
678+
# add final implied CLOSEPOLY, if necessary
679+
ifstart_pointisnotNone \
680+
andnotnp.all(np.isclose(start_point,prev_point)):
681+
B=BezierSegment(np.array([prev_point,start_point]))
682+
area+=B.arc_area
683+
returnarea
684+
629685
definterpolated(self,steps):
630686
"""
631687
Return a new path resampled to length N x steps.
Lines changed: 26 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,26 @@
1+
frommatplotlib.tests.test_pathimport_test_curves
2+
3+
importnumpyasnp
4+
importpytest
5+
6+
7+
# all tests here are currently comparing against integrals
8+
integrate=pytest.importorskip('scipy.integrate')
9+
10+
11+
# get several curves to test our code on by borrowing the tests cases used in
12+
# `~.tests.test_path`. get last path element ([-1]) and curve, not code ([0])
13+
_test_curves= [list(tc.path.iter_bezier())[-1][0]fortcin_test_curves]
14+
15+
16+
def_integral_arc_area(B):
17+
"""(Signed) area swept out by ray from origin to curve."""
18+
dB=B.differentiate(B)
19+
defintegrand(t):
20+
returnnp.cross(B(t),dB(t))/2
21+
returnintegrate.quad(integrand,0,1)[0]
22+
23+
24+
@pytest.mark.parametrize("B",_test_curves)
25+
deftest_area_formula(B):
26+
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
importcopy
22
importre
3+
fromcollectionsimportnamedtuple
34

45
importnumpyasnp
5-
66
fromnumpy.testingimportassert_array_equal
77
importpytest
88

@@ -71,25 +71,27 @@ def test_contains_points_negative_radius():
7171
np.testing.assert_equal(result, [True,False,False])
7272

7373

74-
_test_paths= [
74+
_ExampleCurve=namedtuple('ExampleCurve', ['path','extents','area'])
75+
_test_curves= [
7576
# interior extrema determine extents and degenerate derivative
76-
Path([[0,0], [1,0], [1,1], [0,1]],
77-
[Path.MOVETO,Path.CURVE4,Path.CURVE4,Path.CURVE4]),
78-
# a quadratic curve
79-
Path([[0,0], [0,1], [1,0]], [Path.MOVETO,Path.CURVE3,Path.CURVE3]),
77+
_ExampleCurve(Path([[0,0], [1,0], [1,1], [0,1]],
78+
[Path.MOVETO,Path.CURVE4,Path.CURVE4,Path.CURVE4]),
79+
extents=(0.,0.,0.75,1.),area=0.6),
80+
# a quadratic curve, clockwise
81+
_ExampleCurve(Path([[0,0], [0,1], [1,0]],
82+
[Path.MOVETO,Path.CURVE3,Path.CURVE3]),
83+
extents=(0.,0.,1.,0.5),area=-1/3),
8084
# a linear curve, degenerate vertically
81-
Path([[0,1], [1,1]], [Path.MOVETO,Path.LINETO]),
85+
_ExampleCurve(Path([[0,1], [1,1]], [Path.MOVETO,Path.LINETO]),
86+
extents=(0.,1.,1.,1.),area=0.),
8287
# a point
83-
Path([[1,2]], [Path.MOVETO]),
88+
_ExampleCurve(Path([[1,2]], [Path.MOVETO]),extents=(1.,2.,1.,2.),
89+
area=0.),
8490
]
8591

8692

87-
_test_path_extents= [(0.,0.,0.75,1.), (0.,0.,1.,0.5), (0.,1.,1.,1.),
88-
(1.,2.,1.,2.)]
89-
90-
91-
@pytest.mark.parametrize('path, extents',zip(_test_paths,_test_path_extents))
92-
deftest_exact_extents(path,extents):
93+
@pytest.mark.parametrize('precomputed_curve',_test_curves)
94+
deftest_exact_extents(precomputed_curve):
9395
# notice that if we just looked at the control points to get the bounding
9496
# box of each curve, we would get the wrong answers. For example, for
9597
# hard_curve = Path([[0, 0], [1, 0], [1, 1], [0, 1]],
@@ -99,9 +101,32 @@ def test_exact_extents(path, extents):
99101
# the way out to the control points.
100102
# Note that counterintuitively, path.get_extents() returns a Bbox, so we
101103
# have to get that Bbox's `.extents`.
104+
path,extents=precomputed_curve.path,precomputed_curve.extents
102105
assertnp.all(path.get_extents().extents==extents)
103106

104107

108+
@pytest.mark.parametrize('precomputed_curve',_test_curves)
109+
deftest_signed_area(precomputed_curve):
110+
path,area=precomputed_curve.path,precomputed_curve.area
111+
assertnp.isclose(path.signed_area(),area)
112+
# now flip direction, sign of *signed_area* should flip
113+
rcurve=Path(path.vertices[::-1],path.codes)
114+
assertnp.isclose(rcurve.signed_area(),-area)
115+
116+
117+
deftest_signed_area_unit_rectangle():
118+
rect=Path.unit_rectangle()
119+
assertnp.isclose(rect.signed_area(),1)
120+
121+
122+
deftest_signed_area_unit_circle():
123+
circ=Path.unit_circle()
124+
# Not a "real" circle, just an approximation of a circle made out of bezier
125+
# curves. The actual value is 3.1415935732517166, which is close enough to
126+
# pass here.
127+
assertnp.isclose(circ.signed_area(),np.pi)
128+
129+
105130
deftest_point_in_path_nan():
106131
box=np.array([[0,0], [1,0], [1,1], [0,1], [0,0]])
107132
p=Path(box)

‎requirements/testing/travis_all.txt‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -9,3 +9,4 @@ pytest-timeout
99
pytest-xdist
1010
python-dateutil
1111
tornado
12+
scipy

0 commit comments

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