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

Commit08e61bc

Browse files
committed
add function to compute (signed) area of path
1 parent6654144 commit08e61bc

File tree

3 files changed

+155
-4
lines changed

3 files changed

+155
-4
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
@@ -654,6 +654,41 @@ def intersects_bbox(self, bbox, filled=True):
654654
return_path.path_intersects_rectangle(
655655
self,bbox.x0,bbox.y0,bbox.x1,bbox.y1,filled)
656656

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

‎lib/matplotlib/tests/test_path.py

Lines changed: 16 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -55,6 +55,22 @@ def test_exact_extents_cubic():
5555
np.testing.assert_equal(hard_curve.get_exact_extents(), [0.,0.,0.75,1.])
5656

5757

58+
deftest_signed_area_unit_circle():
59+
circ=Path.unit_circle()
60+
# not quite pi...since it's not quite a circle!
61+
assert(np.isclose(circ.signed_area(),3.1415935732517166))
62+
# now counter-clockwise
63+
rverts=circ.vertices[-2::-1]
64+
rverts=np.append(rverts,np.atleast_2d(circ.vertices[0]),axis=0)
65+
rcirc=Path(rverts,circ.codes)
66+
assert(np.isclose(rcirc.signed_area(),-3.1415935732517166))
67+
68+
69+
deftest_signed_area_unit_rectangle():
70+
rect=Path.unit_rectangle()
71+
assert(np.isclose(rect.signed_area(),1))
72+
73+
5874
deftest_point_in_path_nan():
5975
box=np.array([[0,0], [1,0], [1,1], [0,1], [0,0]])
6076
p=Path(box)

0 commit comments

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