From 4fc7e8592ea8139c3ca9a294c62bf7331263d270 Mon Sep 17 00:00:00 2001 From: mozman Date: Sat, 11 Jan 2025 19:07:44 +0100 Subject: [PATCH] add BSpline.point_inversion() method --- docs/source/math/core.rst | 2 + notes/pages/CHANGELOG.md | 1 + src/ezdxf/math/bspline.py | 85 ++++++++++++++++++++ tests/test_06_math/test_628_bspline_tools.py | 26 +++++- 4 files changed, 111 insertions(+), 3 deletions(-) diff --git a/docs/source/math/core.rst b/docs/source/math/core.rst index e32051e11..3bc772678 100644 --- a/docs/source/math/core.rst +++ b/docs/source/math/core.rst @@ -1143,6 +1143,8 @@ BSpline .. automethod:: degree_elevation + .. automethod:: point_inversion + Bezier ------ diff --git a/notes/pages/CHANGELOG.md b/notes/pages/CHANGELOG.md index 4e08c49b3..242e17da1 100644 --- a/notes/pages/CHANGELOG.md +++ b/notes/pages/CHANGELOG.md @@ -5,6 +5,7 @@ - NEW: `Polyline.points_in_wcs()` method - NEW: `EdgePath.close_gaps()` method - NEW: `BSpline.degree_elevation()` method + - NEW: `BSpline.point_inversion()` method - NEW: `BSpline.insert_knot()` and `BSpline.knot_refinement()` supports rational splines - BUGFIX: Exported `MESH` entities without vertices or faces create invalid DXF files - {{issue 1219}} diff --git a/src/ezdxf/math/bspline.py b/src/ezdxf/math/bspline.py index 459e3000e..0e8a9edbb 100644 --- a/src/ezdxf/math/bspline.py +++ b/src/ezdxf/math/bspline.py @@ -1363,6 +1363,28 @@ def degree_elevation(self, t: int) -> BSpline: """ return degree_elevation(self, t) + def point_inversion( + self, point: UVec, *, epsilon=1e-8, max_iterations=100, init=8 + ) -> float: + """Returns the parameter t for a point on the curve that is closest to the input + point. + + This is an iterative search using Newton's method, so there is no guarantee + of success, especially for splines with many turns. + + Args: + point(UVec): point on the curve or near the curve + epsilon(float): desired precision (distance input point to point on curve) + max_iterations(int): max iterations for Newton's method + init(int): number of points to calculate in the initialization phase + + .. versionadded:: 1.4 + + """ + return point_inversion( + self, Vec3(point), epsilon=epsilon, max_iterations=max_iterations, init=init + ) + def subdivide_params(p: list[float]) -> Iterable[float]: for i in range(len(p) - 1): @@ -1824,3 +1846,66 @@ def from_homogeneous_points( points.append(Vec3(point[:3]) / w) weights.append(w) return points, weights + + +def point_inversion( + spline: BSpline, point: Vec3, *, epsilon=1e-8, max_iterations=100, init=8 +) -> float: + """Returns the parameter t for a point on the curve that is closest to the input + point. + + This is an iterative search using Newton's method, so there is no guarantee + of success, especially for splines with many turns. + + Args: + spline(BSpline): curve + point(Vec3): point on the curve or near the curve + epsilon(float): desired precision (distance input point to point on curve) + max_iterations(int): max iterations for Newton's method + init(int): number of points to calculate in the initialization phase + + .. versionadded:: 1.4 + + """ + max_t = spline.max_t + prev_distance = float("inf") + u = max_t / 2 + + # Initialization phase + t = np.linspace(0, max_t, init, endpoint=True) + chk_points = list(spline.points(t)) + for u1, p in zip(t, chk_points): + distance = point.distance(p) + if distance < prev_distance: + u = float(u1) + prev_distance = distance + + no_progress_counter: int = 0 + for iteration in range(max_iterations): + # Evaluate the B-spline curve at the current parameter value + p, dpdu = spline.derivative(u, n=1) + + # Calculate the difference between the current point and the target point + diff = p - point + + # Check if the difference is within the desired epsilon + distance = diff.magnitude + if distance < epsilon: + break # goal reached + if math.isclose(prev_distance, distance): + no_progress_counter += 1 + if no_progress_counter > 2: + break + else: + no_progress_counter = 0 + prev_distance = distance + + # Update the parameter value using Newton's method + u -= diff.dot(dpdu) / dpdu.dot(dpdu) + + # Clamp the parameter value within the valid range + if u < 0.0: + u = 0.0 + elif u > max_t: + u = max_t + return u diff --git a/tests/test_06_math/test_628_bspline_tools.py b/tests/test_06_math/test_628_bspline_tools.py index b7710c077..2f16af6ef 100644 --- a/tests/test_06_math/test_628_bspline_tools.py +++ b/tests/test_06_math/test_628_bspline_tools.py @@ -5,7 +5,7 @@ import math from ezdxf.math import BSpline from ezdxf.math.bspline import to_homogeneous_points, from_homogeneous_points - +import numpy as np CONTROL_POINTS = [(0, 0), (1, -1), (2, 0), (3, 2), (4, 0), (5, -2)] WEIGHTS = [1, 2, 3, 3, 2, 1] @@ -54,7 +54,7 @@ EXPECTED_WEIGHTS = [1.0, 1.75, 2.25, 2.875, 3.0, 2.8749999999999996, 2.25, 1.75, 1.0] -class TestNonRationalBSpline: +class TestDegreeElevationNonRationalBSpline: @pytest.fixture(scope="class") def result(self): spline = BSpline(CONTROL_POINTS) @@ -85,7 +85,7 @@ def test_elevation_0_times(self): assert spline.degree_elevation(-1) is spline -class TestRationalBSpline: +class TestDegreeElevationRationalBSpline: @pytest.fixture(scope="class") def result(self): spline = BSpline(CONTROL_POINTS, weights=WEIGHTS) @@ -145,5 +145,25 @@ def test_revert_hg_points(self): assert all(math.isclose(w, e) for w, e in zip(weights, WEIGHTS)) is True +class TestPointInversion: + @pytest.fixture(scope="class") + def spline(self): + return BSpline(CONTROL_POINTS) + + def test_start_point(self, spline): + assert math.isclose(0, spline.point_inversion((0, 0))) + + def test_end_point(self, spline): + result = spline.point_inversion(CONTROL_POINTS[-1]) + assert math.isclose(spline.max_t, result) + + def test_many_points(self, spline): + t = np.linspace(0, spline.max_t, 13) + points = spline.points(t) + assert all( + math.isclose(spline.point_inversion(p), u) for p, u in zip(points, t) + ) + + if __name__ == "__main__": pytest.main([__file__])