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

Commit296baeb

Browse files
committed
add function to compute (signed) area of path
1 parent34d3763 commit296baeb

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
@@ -200,15 +200,114 @@ def __init__(self, control_points):
200200
coeff= [math.factorial(self._N-1)
201201
// (math.factorial(i)*math.factorial(self._N-1-i))
202202
foriinrange(self._N)]
203-
self._px=self._cpoints.T*coeff
203+
self._px=(self._cpoints.T*coeff).T
204204

205205
def__call__(self,t):
206-
returnself.point_at_t(t)
206+
t=np.array(t)
207+
orders_shape= (1,)*t.ndim+self._orders.shape
208+
t_shape=t.shape+ (1,)# self._orders.ndim == 1
209+
orders=np.reshape(self._orders,orders_shape)
210+
rev_orders=np.reshape(self._orders[::-1],orders_shape)
211+
t=np.reshape(t,t_shape)
212+
return ((1-t)**rev_orders*t**orders) @self._px
207213

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

213312
@property
214313
defcontrol_points(self):
@@ -279,6 +378,7 @@ def axis_aligned_extrema(self):
279378
"""
280379
n=self.degree
281380
Cj=self.polynomial_coefficients
381+
# much faster than .differentiate(self).polynomial_coefficients
282382
dCj=np.atleast_2d(np.arange(1,n+1)).T*Cj[1:]
283383
iflen(dCj)==0:
284384
returnnp.array([]),np.array([])

‎lib/matplotlib/path.py

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

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

‎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