From bf53eafbad06571f189ba19eb79092c31c46d22a Mon Sep 17 00:00:00 2001 From: Bruno Beltran Date: Mon, 23 Mar 2020 11:45:53 -0700 Subject: [PATCH] code to compute bezier segment / path lengths --- lib/matplotlib/bezier.py | 88 +++++++++++++++++++++++++++++ lib/matplotlib/bezier.pyi | 3 + lib/matplotlib/path.py | 22 ++++++++ lib/matplotlib/path.pyi | 1 + lib/matplotlib/tests/test_bezier.py | 16 ++++++ lib/matplotlib/tests/test_path.py | 23 +++++--- 6 files changed, 146 insertions(+), 7 deletions(-) diff --git a/lib/matplotlib/bezier.py b/lib/matplotlib/bezier.py index 1c5061c65c2d..39a15fe07ab1 100644 --- a/lib/matplotlib/bezier.py +++ b/lib/matplotlib/bezier.py @@ -4,6 +4,7 @@ import math import warnings +from collections import deque import numpy as np @@ -220,6 +221,93 @@ def point_at_t(self, t): """ return tuple(self(t)) + def split_at_t(self, t): + """ + Split into two Bezier curves using de casteljau's algorithm. + + Parameters + ---------- + t : float + Point in [0,1] at which to split into two curves + + Returns + ------- + B1, B2 : BezierSegment + The two sub-curves. + """ + new_cpoints = split_de_casteljau(self._cpoints, t) + return BezierSegment(new_cpoints[0]), BezierSegment(new_cpoints[1]) + + def control_net_length(self): + """Sum of lengths between control points""" + L = 0 + N, d = self._cpoints.shape + for i in range(N - 1): + L += np.linalg.norm(self._cpoints[i+1] - self._cpoints[i]) + return L + + def arc_length(self, rtol=1e-4, atol=1e-6): + """ + Estimate the length using iterative refinement. + + Our estimate is just the average between the length of the chord and + the length of the control net. + + Since the chord length and control net give lower and upper bounds + (respectively) on the length, this maximum possible error is tested + against an absolute tolerance threshold at each subdivision. + + However, sometimes this estimator converges much faster than this error + estimate would suggest. Therefore, the relative change in the length + estimate between subdivisions is compared to a relative error tolerance + after each set of subdivisions. + + Parameters + ---------- + rtol : float, default 1e-4 + If :code:`abs(est[i+1] - est[i]) <= rtol * est[i+1]`, we return + :code:`est[i+1]`. + atol : float, default 1e-6 + If the distance between chord length and control length at any + point falls below this number, iteration is terminated. + """ + chord = np.linalg.norm(self._cpoints[-1] - self._cpoints[0]) + net = self.control_net_length() + max_err = (net - chord)/2 + curr_est = chord + max_err + # early exit so we don't try to "split" paths of zero length + if max_err < atol: + return curr_est + + prev_est = np.inf + curves = deque([self]) + errs = deque([max_err]) + lengths = deque([curr_est]) + while np.abs(curr_est - prev_est) > rtol * curr_est: + # subdivide the *whole* curve before checking relative convergence + # again + prev_est = curr_est + num_curves = len(curves) + for i in range(num_curves): + curve = curves.popleft() + new_curves = curve.split_at_t(0.5) + max_err -= errs.popleft() + curr_est -= lengths.popleft() + for ncurve in new_curves: + chord = np.linalg.norm( + ncurve._cpoints[-1] - ncurve._cpoints[0]) + net = ncurve.control_net_length() + nerr = (net - chord)/2 + nlength = chord + nerr + max_err += nerr + curr_est += nlength + curves.append(ncurve) + errs.append(nerr) + lengths.append(nlength) + if max_err < atol: + return curr_est + return curr_est + @property def arc_area(self): r""" diff --git a/lib/matplotlib/bezier.pyi b/lib/matplotlib/bezier.pyi index 5900287a318e..0447057ebda0 100644 --- a/lib/matplotlib/bezier.pyi +++ b/lib/matplotlib/bezier.pyi @@ -35,6 +35,9 @@ class BezierSegment: def __init__(self, control_points: ArrayLike) -> None: ... def __call__(self, t: ArrayLike) -> np.ndarray: ... def point_at_t(self, t: float) -> tuple[float, ...]: ... + def split_at_t(self, t: float) -> tuple[BezierSegment, BezierSegment]: ... + def control_net_length(self) -> float: ... + def arc_length(self, rtol: float, atol: float) -> float: ... @property def arc_area(self) -> float: ... @classmethod diff --git a/lib/matplotlib/path.py b/lib/matplotlib/path.py index 0ef152c7dd8d..658ce2f82730 100644 --- a/lib/matplotlib/path.py +++ b/lib/matplotlib/path.py @@ -666,6 +666,28 @@ def intersects_bbox(self, bbox, filled=True): return _path.path_intersects_rectangle( self, bbox.x0, bbox.y0, bbox.x1, bbox.y1, filled) + def length(self, rtol=None, atol=None, **kwargs): + r""" + Get length of Path. + + Equivalent to (but not computed as) + + .. math:: + + \sum_{j=1}^N \int_0^1 ||B'_j(t)|| dt + + where the sum is over the :math:`N` Bezier curves that comprise the + Path. Notice that this measure of length will assign zero weight to all + isolated points on the Path. + + Returns + ------- + length : float + The path length. + """ + return np.sum([B.arc_length(rtol, atol) + for B, code in self.iter_bezier(**kwargs)]) + def signed_area(self): """ Get signed area of the filled path. diff --git a/lib/matplotlib/path.pyi b/lib/matplotlib/path.pyi index 3bd0689e79a2..b520643c4f50 100644 --- a/lib/matplotlib/path.pyi +++ b/lib/matplotlib/path.pyi @@ -90,6 +90,7 @@ class Path: def get_extents(self, transform: Transform | None = ..., **kwargs) -> Bbox: ... def intersects_path(self, other: Path, filled: bool = ...) -> bool: ... def intersects_bbox(self, bbox: Bbox, filled: bool = ...) -> bool: ... + def length(self, rtol: float | None, atol: float | None, **kwargs) -> float: ... def signed_area(self) -> float: ... def interpolated(self, steps: int) -> Path: ... def to_polygons( diff --git a/lib/matplotlib/tests/test_bezier.py b/lib/matplotlib/tests/test_bezier.py index 1b255ca0542c..5a2b064970c4 100644 --- a/lib/matplotlib/tests/test_bezier.py +++ b/lib/matplotlib/tests/test_bezier.py @@ -26,6 +26,15 @@ def test_split_bezier_with_large_values(): _test_curves = [list(tc.path.iter_bezier())[-1][0] for tc in _test_curves] +def _integral_arc_length(B): + dB = B.differentiate(B) + def integrand(t): + return np.linalg.norm(dB(t), axis=1) + x = np.linspace(0, 1, 1000) + y = integrand(x) + return np.trapz(y, x) + + def _integral_arc_area(B): """(Signed) area swept out by ray from origin to curve.""" dB = B.differentiate(B) @@ -39,3 +48,10 @@ def integrand(t): @pytest.mark.parametrize("B", _test_curves) def test_area_formula(B): assert np.isclose(_integral_arc_area(B), B.arc_area) + + +@pytest.mark.parametrize("B", _test_curves) +def test_length_iteration(B): + assert np.isclose(_integral_arc_length(B), + B.arc_length(rtol=1e-5, atol=1e-8), + rtol=1e-5, atol=1e-8) diff --git a/lib/matplotlib/tests/test_path.py b/lib/matplotlib/tests/test_path.py index 78ddc2756e2c..789364a5a724 100644 --- a/lib/matplotlib/tests/test_path.py +++ b/lib/matplotlib/tests/test_path.py @@ -88,22 +88,24 @@ def test_contains_points_negative_radius(): np.testing.assert_equal(result, [True, False, False]) -_ExampleCurve = namedtuple('_ExampleCurve', ['path', 'extents', 'area']) +_ExampleCurve = namedtuple('_ExampleCurve', + ['path', 'extents', 'area', 'length']) _test_curves = [ # interior extrema determine extents and degenerate derivative _ExampleCurve(Path([[0, 0], [1, 0], [1, 1], [0, 1]], [Path.MOVETO, Path.CURVE4, Path.CURVE4, Path.CURVE4]), - extents=(0., 0., 0.75, 1.), area=0.6), + extents=(0., 0., 0.75, 1.), area=0.6, length=2.0), # a quadratic curve, clockwise - _ExampleCurve(Path([[0, 0], [0, 1], [1, 0]], - [Path.MOVETO, Path.CURVE3, Path.CURVE3]), - extents=(0., 0., 1., 0.5), area=-1/3), + _ExampleCurve(Path([[0, 0], [0, 1], [1, 0]], [Path.MOVETO, Path.CURVE3, + Path.CURVE3]), extents=(0., 0., 1., 0.5), area=-1/3, + length=(1/25)*(10 + 15*np.sqrt(2) + np.sqrt(5) + * (np.arcsinh(2) + np.arcsinh(3)))), # a linear curve, degenerate vertically _ExampleCurve(Path([[0, 1], [1, 1]], [Path.MOVETO, Path.LINETO]), - extents=(0., 1., 1., 1.), area=0.), + extents=(0., 1., 1., 1.), area=0., length=1.0), # a point _ExampleCurve(Path([[1, 2]], [Path.MOVETO]), extents=(1., 2., 1., 2.), - area=0.), + area=0., length=0.0), ] @@ -154,6 +156,13 @@ def test_signed_area_unit_circle(): assert np.isclose(circ.signed_area(), np.pi) +@pytest.mark.parametrize('precomputed_curve', _test_curves) +def test_length_curve(precomputed_curve): + path, length = precomputed_curve.path, precomputed_curve.length + assert np.isclose(path.length(rtol=1e-5, atol=1e-8), length, rtol=1e-5, + atol=1e-8) + + def test_point_in_path_nan(): box = np.array([[0, 0], [1, 0], [1, 1], [0, 1], [0, 0]]) p = Path(box)