8000 code to compute bezier segment / path lengths · matplotlib/matplotlib@0dbef92 · GitHub
[go: up one dir, main page]
Skip to content

Commit 0dbef92

Browse files
brunobeltranbrunobeltran
committed
code to compute bezier segment / path lengths
1 parent a568413 commit 0dbef92

File tree

4 files changed

+141
-5
lines changed

4 files changed

+141
-5
lines changed

‎lib/matplotlib/bezier.py

Lines changed: 91 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -4,6 +4,7 @@
44

55
import math
66
import warnings
7+
from collections import deque
78

89
import numpy as np
910

@@ -219,6 +220,96 @@ def point_at_t(self, t):
219220
"""Evaluate curve at a single point *t*. Returns a Tuple[float*d]."""
220221
return tuple(self(t))
221222

223+
def split_at_t(self, t):
224+
"""Split into two Bezier curves using de casteljau's algorithm.
225+
226+
Parameters
227+
----------
228+
t : float
229+
Point in [0,1] at which to split into two curves
230+
231+
Returns
232+
-------
233+
B1, B2 : BezierSegment
234+
The two sub-curves.
235+
"""
236+
new_cpoints = split_de_casteljau(self._cpoints, t)
237+
return BezierSegment(new_cpoints[0]), BezierSegment(new_cpoints[1])
238+
239+
def control_net_length(self):
240+
"""Sum of lengths between control points"""
241+
L = 0
242+
N, d = self._cpoints.shape
243+
for i in range(N - 1):
244+
L += np.linalg.norm(self._cpoints[i+1] - self._cpoints[i])
245+
return L
246+
247+
def arc_length(self, rtol=None, atol=None):
248+
"""Estimate the length using iterative refinement.
249+
250+
Our estimate is just the average between the length of the chord and
251+
the length of the control net.
252+
253+
Since the chord length and control net give lower and upper bounds
254+
(respectively) on the length, this maximum possible error is tested
255+
against an absolute tolerance threshold at each subdivision.
256+
257+
However, sometimes this estimator converges much faster than this error
258+
esimate would suggest. Therefore, the relative change in the length
259+
estimate between subdivisions is compared to a relative error tolerance
260+
after each set of subdivisions.
261+
262+
Parameters
263+
----------
264+
rtol : float, default 1e-4
265+
If :code:`abs(est[i+1] - est[i]) <= rtol * est[i+1]`, we return
266+
:code:`est[i+1]`.
267+
atol : float, default 1e-6
268+
If the distance between chord length and control length at any
269+
point falls below this number, iteration is terminated.
270+
"""
271+
if rtol is None:
272+
rtol = 1e-4
273+
if atol is None:
274+
atol = 1e-6
275+
276+
chord = np.linalg.norm(self._cpoints[-1] - self._cpoints[0])
277+
net = self.control_net_length()
278+
max_err = (net - chord)/2
279+
curr_est = chord + max_err
280+
# early exit so we don't try to "split" paths of zero length
281+
if max_err < atol:
282+
return curr_est
283+
284+
prev_est = np.inf
285+
curves = deque([self])
286+
errs = deque([max_err])
287+
lengths = deque([curr_est])
288+
while np.abs(curr_est - prev_est) > rtol * curr_est:
289+
# subdivide the *whole* curve before checking relative convergence
290+
# again
291+
prev_est = curr_est
292+
num_curves = len(curves)
293+
for i in range(num_curves):
294+
curve = curves.popleft()
295+
new_curves = curve.split_at_t(0.5)
296+
max_err -= errs.popleft()
297+
curr_est -= lengths.popleft()
298+
for ncurve in new_curves:
299+
chord = np.linalg.norm(
300+
ncurve._cpoints[-1] - ncurve._cpoints[0])
301+
net = ncurve.control_net_length()
302+
nerr = (net - chord)/2
303+
nlength = chord + nerr
304+
max_err += nerr
305+
curr_est += nlength
306+
curves.append(ncurve)
307+
errs.append(nerr)
308+
lengths.append(nlength)
309+
if max_err < atol:
310+
return curr_est
311+
return curr_est
312+
222313
def arc_area(self):
223314
r"""
224315
(Signed) area swept out by ray from origin to curve.

‎lib/matplotlib/path.py

Lines changed: 21 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -642,6 +642,27 @@ 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+
def length(self, rtol=None, atol=None, **kwargs):
646+
r"""Get length of Path.
647+
648+
Equivalent to (but not computed as)
649+
650+
.. math::
651+
652+
\sum_{j=1}^N \int_0^1 ||B'_j(t)|| dt
653+
654+
where the sum is over the :math:`N` Bezier curves that comprise the
655+
Path. Notice that this measure of length will assign zero weight to all
656+
isolated points on the Path.
657+
658+
Returns
659+
-------
660+
length : float
661+
The path length.
662+
"""
663+
return np.sum([B.arc_length(rtol, atol)
664+
for B, code in self.iter_bezier(**kwargs)])
665+
645666
def signed_area(self, **kwargs):
646667
"""
647668
Get signed area of the filled path.

‎lib/matplotlib/tests/test_bezier.py

Lines changed: 14 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -10,6 +10,13 @@
1010
_test_curves = [list(tc.path.iter_bezier())[-1][0] for tc in _test_curves]
1111

1212

13+
def _integral_arc_length(B):
14+
dB = B.differentiate(B)
15+
def integrand(t):
16+
return np.linalg.norm(dB(t))
17+
return scipy.integrate.quad(integrand, 0, 1)[0]
18+
19+
1320
def _integral_arc_area(B):
1421
"""(Signed) area swept out by ray from origin to curve."""
1522
dB = B.differentiate(B)
@@ -23,3 +30,10 @@ def integrand(t):
2330
@pytest.mark.parametrize("B", _test_curves)
2431
def test_area_formula(B):
2532
assert np.isclose(_integral_arc_area(B), B.arc_area())
33+
34+
35+
@pytest.mark.parametrize("B", _test_curves)
36+
def test_length_iteration(B):
37+
assert np.isclose(_integral_arc_length(B),
38+
B.arc_length(rtol=1e-5, atol=1e-8),
39+
rtol=1e-5, atol=1e-8)

‎lib/matplotlib/tests/test_path.py

Lines changed: 15 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -49,21 +49,24 @@ def test_contains_points_negative_radius():
4949
np.testing.assert_equal(result, [True, False, False])
5050

5151

52-
_ExampleCurve = namedtuple('ExampleCurve', ['path', 'extents', 'area'])
52+
_ExampleCurve = namedtuple('ExampleCurve',
53+
['path', 'extents', 'area', 'length'])
5354
_test_curves = [
5455
# interior extrema determine extents and degenerate derivative
5556
_ExampleCurve(Path([[0, 0], [1, 0], [1, 1], [0, 1]],
5657
[Path.MOVETO, Path.CURVE4, Path.CURVE4, Path.CURVE4]),
57-
extents=(0., 0., 0.75, 1.), area=0.6),
58+
extents=(0., 0., 0.75, 1.), area=0.6, length=2.0),
5859
# a quadratic curve, clockwise
5960
_ExampleCurve(Path([[0, 0], [0, 1], [1, 0]], [Path.MOVETO, Path.CURVE3,
60-
Path.CURVE3]), extents=(0., 0., 1., 0.5), area=-1/3),
61+
Path.CURVE3]), extents=(0., 0., 1., 0.5), area=-1/3,
62+
length=(1/25)*(10 + 15*np.sqrt(2) + np.sqrt(5) \
63+
* (np.arcsinh(2) + np.arcsinh(3)))),
6164
# a linear curve, degenerate vertically
6265
_ExampleCurve(Path([[0, 1], [1, 1]], [Path.MOVETO, Path.LINETO]),
63-
extents=(0., 1., 1., 1.), area=0.),
66+
extents=(0., 1., 1., 1.), area=0., length=1.0),
6467
# a point
6568
_ExampleCurve(Path([[1, 2]], [Path.MOVETO]), extents=(1., 2., 1., 2.),
66-
area=0.),
69+
area=0., length=0.0),
6770
]
6871

6972

@@ -105,6 +108,13 @@ def test_signed_area_unit_circle():
105108
assert np.isclose(circ.signed_area(), 3.1415935732517166)
106109

107110

111+
@pytest.mark.parametrize('precomputed_curve', _test_curves)
112+
def test_length_curve(precomputed_curve):
113+
path, length = precomputed_curve.path, precomputed_curve.length
114+
assert np.isclose(path.length(rtol=1e-5, atol=1e-8), length, rtol=1e-5,
115+
atol=1e-8)
116+
117+
108118
def test_point_in_path_nan():
109119
box = np.array([[0, 0], [1, 0], [1, 1], [0, 1], [0, 0]])
110120
p = Path(box)

0 commit comments

Comments
 (0)
0