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

Commit4a41a2f

Browse files
committed
add function to compute (signed) area of path
1 parent405aba6 commit4a41a2f

File tree

4 files changed

+192
-7
lines changed

4 files changed

+192
-7
lines changed

‎lib/matplotlib/bezier.py

Lines changed: 104 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -191,15 +191,114 @@ def __init__(self, control_points):
191191
coeff= [math.factorial(self._N-1)
192192
// (math.factorial(i)*math.factorial(self._N-1-i))
193193
foriinrange(self._N)]
194-
self._px=self._cpoints.T*coeff
194+
self._px=(self._cpoints.T*coeff).T
195195

196196
def__call__(self,t):
197-
returnself.point_at_t(t)
197+
t=np.array(t)
198+
orders_shape= (1,)*t.ndim+self._orders.shape
199+
t_shape=t.shape+ (1,)# self._orders.ndim == 1
200+
orders=np.reshape(self._orders,orders_shape)
201+
rev_orders=np.reshape(self._orders[::-1],orders_shape)
202+
t=np.reshape(t,t_shape)
203+
return ((1-t)**rev_orders*t**orders) @self._px
198204

199205
defpoint_at_t(self,t):
200206
"""Return the point on the Bezier curve for parameter *t*."""
201-
returntuple(
202-
self._px @ (((1-t)**self._orders)[::-1]*t**self._orders))
207+
returntuple(self(t))
208+
209+
defarc_area(self):
210+
r"""
211+
(Signed) area swept out by ray from origin to curve.
212+
213+
Notes
214+
-----
215+
A simple, analytical formula is possible for arbitrary bezier curves.
216+
217+
Given a bezier curve B(t), in order to calculate the area of the arc
218+
swept out by the ray from the origin to the curve, we simply need to
219+
compute :math:`\frac{1}{2}\int_0^1 B(t) \cdot n(t) dt`, where
220+
:math:`n(t) = u^{(1)}(t) \hat{x}_0 - u{(0)}(t) \hat{x}_1` is the normal
221+
vector oriented away from the origin and :math:`u^{(i)}(t) =
222+
\frac{d}{dt} B^{(i)}(t)` is the :math:`i`th component of the curve's
223+
tangent vector. (This formula can be found by applying the divergence
224+
theorem to :math:`F(x,y) = [x, y]/2`, and calculates the *signed* area
225+
for a counter-clockwise curve, by the right hand rule).
226+
227+
The control points of the curve are just its coefficients in a
228+
Bernstein expansion, so if we let :math:`P_i = [P^{(0)}_i, P^{(1)}_i]`
229+
be the :math:`i`'th control point, then
230+
231+
.. math::
232+
233+
\frac{1}{2}\int_0^1 B(t) \cdot n(t) dt
234+
&= \frac{1}{2}\int_0^1 B^{(0)}(t) \frac{d}{dt} B^{(1)}(t)
235+
- B^{(1)}(t) \frac{d}{dt} B^{(0)}(t)
236+
dt \\
237+
&= \frac{1}{2}\int_0^1
238+
\left( \sum_{j=0}^n P_j^{(0)} b_{j,n} \right)
239+
\left( n \sum_{k=0}^{n-1} (P_{k+1}^{(1)} -
240+
P_{k}^{(1)}) b_{j,n} \right)
241+
\\
242+
&\hspace{1em} - \left( \sum_{j=0}^n P_j^{(1)} b_{j,n}
243+
\right) \left( n \sum_{k=0}^{n-1} (P_{k+1}^{(0)}
244+
- P_{k}^{(0)}) b_{j,n} \right)
245+
dt,
246+
247+
where :math:`b_{\nu, n}(t) = {n \choose \nu} t^\nu {(1 - t)}^{n-\nu}`
248+
is the :math:`\nu`'th Bernstein polynomial of degree :math:`n`.
249+
250+
Grouping :math:`t^l(1-t)^m` terms together for each :math:`l`,
251+
:math:`m`, we get that the integrand becomes
252+
253+
.. math::
254+
255+
\sum_{j=0}^n \sum_{k=0}^{n-1}
256+
{n \choose j} {{n - 1} \choose k}
257+
&\left[P_j^{(0)} (P_{k+1}^{(1)} - P_{k}^{(1)})
258+
- P_j^{(1)} (P_{k+1}^{(0)} - P_{k}^{(0)})\right] \\
259+
&\hspace{1em}\times{}t^{j + k} {(1 - t)}^{2n - 1 - j - k}
260+
261+
or just
262+
263+
.. math::
264+
265+
\sum_{j=0}^n \sum_{k=0}^{n-1}
266+
\frac{{n \choose j} {{n - 1} \choose k}}
267+
{{{2n - 1} \choose {j+k}}}
268+
[P_j^{(0)} (P_{k+1}^{(1)} - P_{k}^{(1)})
269+
- P_j^{(1)} (P_{k+1}^{(0)} - P_{k}^{(0)})]
270+
b_{j+k,2n-1}(t).
271+
272+
Interchanging sum and integral, and using the fact that :math:`\int_0^1
273+
b_{\nu, n}(t) dt = \frac{1}{n + 1}`, we conclude that the
274+
original integral can
275+
simply be written as
276+
277+
.. math::
278+
279+
\frac{1}{2}&\int_0^1 B(t) \cdot n(t) dt
280+
\\
281+
&= \frac{1}{4}\sum_{j=0}^n \sum_{k=0}^{n-1}
282+
\frac{{n \choose j} {{n - 1} \choose k}}
283+
{{{2n - 1} \choose {j+k}}}
284+
[P_j^{(0)} (P_{k+1}^{(1)} - P_{k}^{(1)})
285+
- P_j^{(1)} (P_{k+1}^{(0)} - P_{k}^{(0)})]
286+
"""
287+
n=self.degree
288+
area=0
289+
P=self.control_points
290+
dP=np.diff(P,axis=0)
291+
forjinrange(n+1):
292+
forkinrange(n):
293+
area+=_comb(n,j)*_comb(n-1,k)/_comb(2*n-1,j+k) \
294+
* (P[j,0]*dP[k,1]-P[j,1]*dP[k,0])
295+
return (1/4)*area
296+
297+
@classmethod
298+
defdifferentiate(cls,B):
299+
"""Return the derivative of a BezierSegment, itself a BezierSegment"""
300+
dcontrol_points=B.degree*np.diff(B.control_points,axis=0)
301+
returncls(dcontrol_points)
203302

204303
@property
205304
defcontrol_points(self):
@@ -270,6 +369,7 @@ def axis_aligned_extrema(self):
270369
"""
271370
n=self.degree
272371
Cj=self.polynomial_coefficients
372+
# much faster than .differentiate(self).polynomial_coefficients
273373
dCj=np.atleast_2d(np.arange(1,n+1)).T*Cj[1:]
274374
iflen(dCj)==0:
275375
returnnp.array([]),np.array([])

‎lib/matplotlib/path.py

Lines changed: 35 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -642,6 +642,41 @@ def intersects_bbox(self, bbox, filled=True):
642642
return_path.path_intersects_rectangle(
643643
self,bbox.x0,bbox.y0,bbox.x1,bbox.y1,filled)
644644

645+
defsigned_area(self,**kwargs):
646+
"""
647+
Get signed area filled by path.
648+
649+
All sub paths are treated as if they had been closed. That is, if there
650+
is a MOVETO without a preceding CLOSEPOLY, one is added.
651+
652+
Signed area means that if a path is self-intersecting, the drawing rule
653+
"even-odd" is used and only the filled area is counted.
654+
655+
Returns
656+
-------
657+
area : float
658+
The (signed) enclosed area of the path.
659+
"""
660+
area=0
661+
prev_point=None
662+
prev_code=None
663+
start_point=None
664+
forB,codeinself.iter_bezier(**kwargs):
665+
ifcode==Path.MOVETO:
666+
ifprev_codeisnotNoneandprev_codeisnotPath.CLOSEPOLY:
667+
Bclose=BezierSegment(np.array([prev_point,start_point]))
668+
area+=Bclose.arc_area()
669+
start_point=B.control_points[0]
670+
area+=B.arc_area()
671+
prev_point=B.control_points[-1]
672+
prev_code=code
673+
# add final implied CLOSEPOLY, if necessary
674+
ifstart_pointisnotNone \
675+
andnotnp.all(np.isclose(start_point,prev_point)):
676+
B=BezierSegment(np.array([prev_point,start_point]))
677+
area+=B.arc_area()
678+
returnarea
679+
645680
definterpolated(self,steps):
646681
"""
647682
Return a new path resampled to length N x steps.

‎lib/matplotlib/tests/test_bezier.py

Lines changed: 35 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,35 @@
1+
frommatplotlibimportbezier
2+
3+
importnumpyasnp
4+
importscipy.integrate
5+
6+
7+
test_curves= [
8+
# cubic, degenerate derivative
9+
bezier.BezierSegment([[0,0], [1,0], [1,1], [0,1]]),
10+
# cubic, could be quadratic
11+
bezier.BezierSegment([[0,0], [1,1/2], [1,1/2], [0,1]]),
12+
# interior extrema determine extents
13+
bezier.BezierSegment([[0,0], [4,0], [-4,1], [0,1]]),
14+
# a quadratic curve
15+
bezier.BezierSegment([[0,0], [0,1], [1,1]]),
16+
# a linear curve
17+
bezier.BezierSegment([[0,1], [1,1]]),
18+
# a point
19+
bezier.BezierSegment([[1,2]])
20+
]
21+
22+
23+
def_integral_arc_area(B):
24+
"""(Signed) area swept out by ray from origin to curve."""
25+
dB=B.differentiate(B)
26+
defintegrand(t):
27+
dr=dB(t).T
28+
n=np.array([dr[1],-dr[0]])
29+
return (B(t).T @n)/2
30+
returnscipy.integrate.quad(integrand,0,1)[0]
31+
32+
33+
deftest_area_formula():
34+
forBintest_curves:
35+
assert(np.isclose(_integral_arc_area(B),B.arc_area()))

‎lib/matplotlib/tests/test_path.py

Lines changed: 18 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -50,10 +50,25 @@ def test_contains_points_negative_radius():
5050
np.testing.assert_equal(result, [True,False,False])
5151

5252

53+
# convenient curve (xmax=0.75, area=0.6, length=2.0)
54+
_hard_curve=Path([[0,0], [1,0], [1,1], [0,1]],
55+
[Path.MOVETO,Path.CURVE4,Path.CURVE4,Path.CURVE4])
5356
deftest_exact_extents_cubic():
54-
hard_curve=Path([[0,0], [1,0], [1,1], [0,1]],
55-
[Path.MOVETO,Path.CURVE4,Path.CURVE4,Path.CURVE4])
56-
assert(hard_curve.get_extents().bounds== (0.,0.,0.75,1.))
57+
assert(_hard_curve.get_extents().bounds== (0.,0.,0.75,1.))
58+
59+
60+
deftest_signed_area_curve():
61+
assert(np.isclose(_hard_curve.signed_area(),0.6))
62+
# now counter-clockwise
63+
rverts=_hard_curve.vertices[:0:-1]
64+
rverts=np.append(rverts,np.atleast_2d(_hard_curve.vertices[0]),axis=0)
65+
rcurve=Path(rverts,_hard_curve.codes)
66+
assert(np.isclose(rcurve.signed_area(),-0.6))
67+
68+
69+
deftest_signed_area_unit_rectangle():
70+
rect=Path.unit_rectangle()
71+
assert(np.isclose(rect.signed_area(),1))
5772

5873

5974
deftest_point_in_path_nan():

0 commit comments

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