add BSpline.point_inversion() method

docs/source/math/core.rst

@@ -1143,6 +1143,8 @@ BSpline
 
     .. automethod:: degree_elevation
 
+    .. automethod:: point_inversion
+
 
 Bezier
 ------

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}}

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

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__])