8000 add function to compute (signed) area of path · matplotlib/matplotlib@37d533f · GitHub
[go: up one dir, main page]
Skip to content

Commit 37d533f

Browse files
brunobeltranbrunobeltran
committed
add function to compute (signed) area of path
1 parent d7084de commit 37d533f

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
for i in range(self._N)]
203-
self._px = self._cpoints.T * coeff
203+
self._px = (self._cpoints.T * coeff).T
204204

205205
def __call__(self, t):
206-
return self.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
def point_at_t(self, t):
209215
"""Return the point on the Bezier curve for parameter *t*."""
210-
return tuple(
211-
self._px @ (((1 - t) ** self._orders)[::-1] * t ** self._orders))
216+
return tuple(self(t))
217+
218+
def arc_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 fr 10000 om 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+
for j in range(n + 1):
301+
for k in range(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+
def differentiate(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+
return cls(dcontrol_points)
212311

213312
@property
214313
def control_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
if len(dCj) == 0:
284384
return np.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+
def signed_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+
for cpoints, code in self.iter_curves(**kwargs):
681+
if code == Path.MOVETO:
682+
if prev_code is not None and prev_code is not Path.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+
if start_point is not None \
691+
and not np.all(np.isclose(start_point, prev_point)):
692+
B = BezierSegment(np.array([prev_point, start_point]))
693+
area += B.arc_area()
694+
return area
695+
661696
def interpolated(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+
def test_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+
def test_signed_area_unit_rectangle():
70+
rect = Path.unit_rectangle()
71+
assert(np.isclose(rect.signed_area(), 1))
72+
73+
5874
def test_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)
0