Refactor Bézier intersection methods across modules:

- Rename `bezier_intersect_certified_full` to `bez_ccx` for brevity in `_bez_ccx3.py`.
- Update all associated references and docstrings in test cases and examples.
- Replace `bezier_curve_surface_intersect_certified` with `bez_csx` in `_ncsx2.py` for consistency.
- Integrate `bez_csx` utility in `_ssx4.py` for surface-surface intersection refinements.
- Clean up unused imports and deprecated utilities in `_ssx4.py`.

Author: sth-v

mmcore/geom/_nurbs_eval.py

@@ -393,6 +393,7 @@ def evaluate_nurbs_surface(surface:NURBSSurfaceTuple, u, v, d_order=2)->Evaluate
     nv = len(surface1.control_points[0])
     U = surface1.knot_u[:]  # assume these are already lists/numpy arrays
     V = surface1.knot_v[:]
+
     span_u = _find_span_linear(p, U, nu, u)
     span_v = _find_span_linear(q, V, nv, v)
     # print(p, U, span_u, u, d_order)
@@ -712,13 +713,13 @@ def _join_weights(pts,weights):
 def _nurbs_to_tuple(s1:nurbs.NURBSCurve | nurbs.NURBSSurface)->NURBSCurveTuple | NURBSSurfaceTuple:
     if isinstance(s1,nurbs.NURBSSurface):
 
-        surf1 = NURBSSurfaceTuple(order_u=s1.degree[0] + 1, order_v=s1.degree[1] + 1, knot_u=s1.knots_u.tolist(),
-                              knot_v=s1.knots_v.tolist(), control_points=np.array(s1.control_points),
+        surf1 = NURBSSurfaceTuple(order_u=s1.degree[0] + 1, order_v=s1.degree[1] + 1, knot_u=np.array(s1.knots_u),
+                              knot_v=np.array(s1.knots_v), control_points=np.array(s1.control_points),
                               weights=np.ascontiguousarray(s1.control_points_w[..., -1]))
         return surf1
     elif isinstance(s1,nurbs.NURBSCurve):
         cpts,weights=from_homogeneous_1d(np.array(s1.control_pointsw))
-        curve1 = NURBSCurveTuple(order=s1.degree + 1,knot=s1.knots.tolist(),
+        curve1 = NURBSCurveTuple(order=s1.degree + 1,knot=np.array(s1.knots),
                                   control_points=cpts,
                                   weights=weights)
         return curve1

mmcore/geom/_nurbs_interp.py

@@ -105,7 +105,7 @@ def _remove_adjacent_duplicates(points, tol=1e-5):
     return points[mask]
 
 from typing import Literal
-def interpolate_curve(points, degree,  use_centripetal=False, method:Literal['lstsq','lu']='lstsq',remove_duplicates:bool=False,tol=1e-6):
+def interpolate_curve(points, degree,  use_centripetal=False, method:Literal['lstsq','lu']='lstsq',remove_duplicates:bool=False,tol=1e-6,params=None):
     """ Curve interpolation through the data points.
 
     Please refer to Algorithm A9.1 on The NURBS Book (2nd Edition), pp.369-370 for details.
@@ -119,7 +119,12 @@ def interpolate_curve(points, degree,  use_centripetal=False, method:Literal['ls
     num_points = len(points)
 
     # Get uk
-    uk,total_chord_length = compute_params_curve(points, use_centripetal)
+    if params is None:
+
+        uk,total_chord_length = compute_params_curve(points, use_centripetal)
+    else:
+        uk = params
+
 
     # Compute knot vector
     kv = compute_knot_vector(degree, num_points, uk)

mmcore/numeric/_bern_homog.py

@@ -308,6 +308,7 @@ def eval_bezier_surface_homog_with_derivs(Pw: np.ndarray, u: float, v: float, wa
         Shuv : (d_h,) mixed partial
         Shvv : (d_h,) second partial wrt v
     """
+
     Pw = np.asarray(Pw, dtype=float, order="C")
     u = float(u)
     v = float(v)

mmcore/numeric/cbern.pyx

@@ -82,9 +82,8 @@ def bernstein_basis(int n, double t):
         # Our C logic handles all indices via the recurrence or memset, so empty is safe).
 
         cdef double[:] B = np.empty(n + 1, dtype=np.float64)
-        with nogil:
-            # Pass the raw pointer to the nogil C function
-            bernstein_basis_funs_c(n, t, &B[0])
+
+        bernstein_basis_funs_c(n, t, &B[0])
 
         return B
 

mmcore/numeric/intersection/ccx/_bez_ccx3.py

@@ -1387,7 +1387,7 @@ def _bez_get_tol_adapter(c, tol, rational=False, interval=None):
     return nurbs_curve_param_tolerance(bern_to_nurbs_bezier(c, rational=rational, interval=interval), tol)
 
 
-def bezier_intersect_certified_full(C1: NDArray, C2: NDArray,  sv_thresh: float = 1e-8, atol: float = 1e-3, rational: bool = False,**kwargs) -> IntersectionResult:
+def bez_ccx(C1: NDArray, C2: NDArray, sv_thresh: float = 1e-8, atol: float = 1e-3, rational: bool = False, **kwargs) -> IntersectionResult:
     """Certified intersection for (possibly rational) Bézier curve pairs.
 
     Parameters
@@ -1434,10 +1434,10 @@ def bezier_intersect_certified_full(C1: NDArray, C2: NDArray,  sv_thresh: float
     Planar polynomial case::
 
         >>> import numpy as np
-        >>> from mmcore.numeric.intersection.ccx._bez_ccx3 import bezier_intersect_certified_full
+        >>> from mmcore.numeric.intersection.ccx._bez_ccx3 import bez_ccx
         >>> c1 = np.array([[0., 0., 0.], [1., 1., 0.], [2., 0., 0.]])
         >>> c2 = np.array([[0., 0., 0.], [1., 0., 0.], [2., 0., 0.]])
-        >>> res = bezier_intersect_certified_full(c1, c2)
+        >>> res = bez_ccx(c1, c2)
         >>> len(res['isolated'])
         2
 
@@ -1446,7 +1446,7 @@ def bezier_intersect_certified_full(C1: NDArray, C2: NDArray,  sv_thresh: float
         >>> w = np.sqrt(0.5)
         >>> arc = np.array([[1., 0., 1.], [w, w, w], [0., 1., 1.]])
         >>> line = np.array([[0., 0., 1.], [0.5, 0.5, 1.], [1., 1., 1.]])
-        >>> res = bezier_intersect_certified_full(arc, line, rational=True)
+        >>> res = bez_ccx(arc, line, rational=True)
         >>> np.allclose((res['isolated'][0]['u'], res['isolated'][0]['v']), (0.5, np.sqrt(0.5)))
         True
     """
@@ -1627,7 +1627,8 @@ def cell_contains_known_isolated(u0, u1, v0, v1, margin=1e-9):
             min_d = float(np.min(dnet))
             if  cell_contains_known_isolated(u0, u1, v0, v1)  :
                 res = {"type": "none"}
-                print('bae',    (u0, u1), (v0, v1))
+                print('bae',    (u0, u1), (v0, v1), )
+                continue
             else:
                 res = contact_detect_and_extract(Pseg, Qseg, seed_uv=(0.5, 0.5), sv_thresh=sv_thresh, rational=rational)
 
@@ -1886,7 +1887,7 @@ def get_interv(n) -> tuple[float, float]:
 
     # case 1: One true overlap, no isolated points
     s = time.perf_counter()
-    inter12 = bezier_intersect_certified_full(curve1, curve2)
+    inter12 = bez_ccx(curve1, curve2)
     print(time.perf_counter() - s, "pruned:", inter12["stats"]["pruned"])  # 0.07984466700008852 (current perf)
     print("\n\n case1\n", "-" * 80, "\n")
     assert len(inter12["isolated"]) == 0
@@ -1924,7 +1925,7 @@ def get_interv(n) -> tuple[float, float]:
     # case 2: Two isolated points, no overlaps
     print("\n\n case2\n", "-" * 80, "\n")
     s = time.perf_counter()
-    inter34 = bezier_intersect_certified_full(curve3, curve4)
+    inter34 = bez_ccx(curve3, curve4)
     # 0.02212008300011803  (current perf)
     print(time.perf_counter() - s, "pruned:", inter34["stats"]["pruned"])
     #print(inter34)
@@ -1947,7 +1948,7 @@ def get_interv(n) -> tuple[float, float]:
     print("\n\n case3.1\n", "-" * 80, "\n")
     # case 3.1: Two isolated points, but the curves are close to each other and are almost parallel, which makes this example similar to overlap.
     s = time.perf_counter()
-    inter36 = bezier_intersect_certified_full(curve3, curve6)
+    inter36 = bez_ccx(curve3, curve6)
     print(time.perf_counter() - s, "pruned:", inter36["stats"]["pruned"])  # 1.4956505000000107 current perf
     #print(inter36)  # WRONG! Only one isolated point is returning.
 
@@ -1977,7 +1978,7 @@ def get_interv(n) -> tuple[float, float]:
     print("\n\n case3.2\n", "-" * 80, "\n")
     # case 3.2: Everything is the same, they just swapped the curves around.
     s = time.perf_counter()
-    inter63 = bezier_intersect_certified_full(curve6, curve3)
+    inter63 = bez_ccx(curve6, curve3)
     print(time.perf_counter() - s, "pruned:", inter63["stats"]["pruned"])
     #print(
     #    inter63,
@@ -2015,7 +2016,7 @@ def get_interv(n) -> tuple[float, float]:
     curve8 = np.array([[0, 0,0.], [1, 1,0.], [2, 0,0.]], float) # quadratic arc
 
     s = time.perf_counter()
-    inter78 = bezier_intersect_certified_full(curve7, curve8)
+    inter78 = bez_ccx(curve7, curve8)
     print(time.perf_counter() - s, "pruned:", inter78["stats"]["pruned"])
     #print(
     #    inter78,
@@ -2044,7 +2045,7 @@ def get_interv(n) -> tuple[float, float]:
     curve_line_h = np.array([[0.0, 0.0, 1.0], [1.0, 1.0, 1.0]])
 
     s = time.perf_counter()
-    inter_rational = bezier_intersect_certified_full(curve_arc_h, curve_line_h, rational=True)
+    inter_rational = bez_ccx(curve_arc_h, curve_line_h, rational=True)
     print(time.perf_counter() - s, "pruned:", inter_rational["stats"]["pruned"])
     assert len(inter_rational["isolated"]) == 1
     assert np.allclose(inter_rational["isolated"][0]["point"], [np.sqrt(0.5), np.sqrt(0.5)], atol=1e-9)
@@ -2080,7 +2081,7 @@ def get_interv(n) -> tuple[float, float]:
        [52.914001, 42.130837,  0.      ,  0.707   ],
        [49.76    , 45.703   ,  0.      ,  1.      ]])
     s = time.perf_counter()
-    inter_rational2 = bezier_intersect_certified_full(elliptical_arc_pts_h, arc_pts2_h, rational=True)
+    inter_rational2 = bez_ccx(elliptical_arc_pts_h, arc_pts2_h, rational=True)
     print(time.perf_counter() - s, "pruned:", inter_rational2["stats"]["pruned"])
     assert len(inter_rational2["isolated"]) == 2
 
@@ -2098,7 +2099,7 @@ def get_interv(n) -> tuple[float, float]:
 
     c3d1,c3d2=np.array([[-2.5, -5.0, 0.826], [-2.5, -3.333, 0.826], [-2.5, -1.667, 0.001], [-2.5, 0.0, -0.102]]),np.array([[-5.0, -2.5, 0.052], [-3.333, -2.5, 0.052], [-1.667, -2.5, 0.75], [0.0, -2.5, 0.75]])
     s = time.perf_counter()
-    inter_3d_1 = bezier_intersect_certified_full(c3d1,c3d2, rational=False)
+    inter_3d_1 = bez_ccx(c3d1, c3d2, rational=False)
     print(time.perf_counter() - s, "pruned:", inter_3d_1["stats"]["pruned"])
     print("pruned_by:", set(inter_3d_1["stats"]["pruned_by"]))
     make_rich_table(inter_3d_1, {"isolated": [{'u':0.499982,'v':0.499986,'point':[-2.500035, -2.500044, 0.400817]}], "overlaps": [],'stats':{}}, 'case 7 (3d)',oracle_name='Expected')
@@ -2114,7 +2115,7 @@ def get_interv(n) -> tuple[float, float]:
               [-0.18421653, 0.21499867, 0.],
               [-0.18858128, 0.2141638, 0.]])
     s = time.perf_counter()
-    inter910 = bezier_intersect_certified_full(crv9, crv10, rational=False)
+    inter910 = bez_ccx(crv9, crv10, rational=False)
     print(time.perf_counter() - s, "pruned:", inter910["stats"]["pruned"])
 
     for i in inter910["isolated"]:
@@ -2140,7 +2141,7 @@ def get_interv(n) -> tuple[float, float]:
     crv12 = np.array([[-0.4507911, 0.11900211, 0.],
                      [0.32897481, 0.35801962, 0.],
                      [0.25497926, 0.73647834, 0.]])
-    inter1112 = bezier_intersect_certified_full(crv11, crv12, rational=False)
+    inter1112 = bez_ccx(crv11, crv12, rational=False)
     print(time.perf_counter() - s, "pruned:", inter1112["stats"]["pruned"])
     for i in inter1112["isolated"]:
         print('uv:',[float(i['u']),float(i['v'])], 'point:',i['point'].tolist())

mmcore/numeric/intersection/ccx/_nccx.py

@@ -227,7 +227,7 @@ def _eq(x):
     return ints
 def _is_rational(curve:NURBSCurveTuple):
     return not np.allclose(curve.weights,1)
-from mmcore.numeric.intersection.ccx._bez_ccx3 import bezier_intersect_certified_full,map_local_to_global,IsolatedIntersection,OverlapIntersection
+from mmcore.numeric.intersection.ccx._bez_ccx3 import bez_ccx,map_local_to_global,IsolatedIntersection,OverlapIntersection
 
 
 def nurbs_ccx(curve1: NURBSCurve | NURBSCurveTuple, curve2: NURBSCurve | NURBSCurveTuple, tol: float = 1e-3,
@@ -261,7 +261,7 @@ def nurbs_ccx(curve1: NURBSCurve | NURBSCurveTuple, curve2: NURBSCurve | NURBSCu
             pts1 = _c1.control_points
             pts2 = _c2.control_points
         #print(tol)
-        result=bezier_intersect_certified_full(pts1,  pts2, atol=tol, rational=rational)
+        result=bez_ccx(pts1, pts2, atol=tol, rational=rational)
         if len(result['isolated'])==0 and len(result['overlaps'])==0:
             #print(set(result['stats']['pruned_by']))
             #print(pts1.tolist(),pts2.tolist())
@@ -359,7 +359,7 @@ def get_segm_curve_index(segment_index):
             pts1 = segm1.control_points
             pts2 = segm2.control_points
 
-        result = bezier_intersect_certified_full(pts1, pts2, atol=tol, rational=rational)
+        result = bez_ccx(pts1, pts2, atol=tol, rational=rational)
         if len(result['isolated'])==0 and len(result['overlaps'])==0:
             #print(segm1)
             #print('$')

mmcore/numeric/intersection/csx/__init__.py

@@ -1,11 +1,11 @@
 
 from mmcore.numeric.intersection.csx._ncsx import nurbs_csx
 from mmcore.numeric.intersection.csx._ncsx2 import nurbs_csx_v2
-from mmcore.numeric.intersection.csx._bez_csx3 import bezier_curve_surface_intersect_certified
+from mmcore.numeric.intersection.csx._bez_csx3 import bez_csx
 
 
 
 
 
 
-__all__ = ['nurbs_csx','nurbs_csx_v2','bezier_curve_surface_intersect_certified']
+__all__ = ['nurbs_csx','nurbs_csx_v2', 'bez_csx']

mmcore/numeric/intersection/csx/_bez_csx3.py

@@ -837,13 +837,14 @@ def _check_interval(t_a, t_b, uu, vv):
         iso_bnd = sorted(iso_bnd, key=lambda it: it["t"])
         iso_unique = [iso_bnd[0]]
         for it in iso_bnd[1:]:
+
             if abs(float(it["t"]) - float(iso_unique[-1]["t"])) > t_eps:
                 iso_unique.append(it)
         iso_bnd = iso_unique
     if ovl_bnd:
         if len(ovl_bnd) > 1:
             logger.error(
-                "Boundary overlap returned multiple segments (len=%d) for span confirm.", len(ovl_bnd)
+                "Boundary overlap returned multiple segments (len=%d) for span confirm.", len(ovl_bnd),ovl_bnd
             )
         ovl = dict(ovl_bnd[0])
         ovl["partial"] = True
@@ -885,10 +886,10 @@ def _check_interval(t_a, t_b, uu, vv):
                 "partial": True,
             }
         logger.error(
-            "Span confirm: one endpoint on surface, boundary intersects=%d (expected 1).",
-            len(iso_bnd),
+            "Span confirm: one endpoint on surface, boundary intersects=%d (expected 1).  %s",
+            len(iso_bnd),repr(iso_bnd)
         )
-        raise RuntimeError("curve_surface_boundary_intersect returned unexpected count for span confirm")
+        raise RuntimeError("curve_surface_boundary_intersect returned unexpected count for span confirm {}".format(len(iso_bnd)))
 
     if not on0 and not on1:
         if len(iso_bnd) == 2:
@@ -912,17 +913,28 @@ def _check_interval(t_a, t_b, uu, vv):
                 "end": "boundary",
                 "partial": True,
             }
+        elif len(iso_bnd) == 1:
+            return {
+                "t": iso_bnd[0]["t"],
+                "u": iso_bnd[0]["u"],
+                "v": iso_bnd[0]["v"],
+                "point":  eval_bezier_curve(C,  iso_bnd[0]["t"], rational=rational),
+                "isolated": True,
+
+            }
         logger.error(
-            "Span confirm: no endpoints on surface, boundary intersects=%d (expected 2).", len(iso_bnd)
+            "Span confirm: no endpoints on surface, boundary intersects=%d (expected 2). %s", len(iso_bnd),repr(iso_bnd)
         )
-        raise RuntimeError("curve_surface_boundary_intersect returned unexpected count for span confirm")
+        raise RuntimeError("curve_surface_boundary_intersect returned unexpected count for span confirm {}".format(len(iso_bnd)))
 
     logger.error(
         "Span confirm: unexpected boundary state (on0=%s, on1=%s, iso=%d, ovl=%d).",
         on0,
         on1,
         len(iso_bnd),
         len(ovl_bnd),
+        iso_bnd,
+        ovl_bnd,
     )
     raise RuntimeError("Unexpected span-confirm boundary state")
 
@@ -1654,7 +1666,7 @@ def ch_separability(pts1,pts2, atol,rational=False):
     # print(res,(pts1.tolist(), pts2.tolist()))
     return int(res)
 
-from mmcore.numeric.intersection.ccx._bez_ccx3 import bezier_intersect_certified_full
+from mmcore.numeric.intersection.ccx._bez_ccx3 import bez_ccx
 
 
 # ---------------------------------------------------------------------------
@@ -1737,7 +1749,7 @@ def curve_surface_boundary_intersect(
 
     for bnd_curve, fixed_axis, fixed_val in boundaries:
         # Intersect curve C with boundary curve
-        ccx_result = bezier_intersect_certified_full(
+        ccx_result = bez_ccx(
             C, bnd_curve,
             sv_thresh=sv_thresh,
             atol=atol,
@@ -1808,7 +1820,7 @@ def curve_surface_boundary_intersect(
 check_boundaries_inter = curve_surface_boundary_intersect
 
 
-def bezier_curve_surface_intersect_certified(
+def bez_csx(
     C: NDArray,
     S: NDArray,
     sv_thresh: float = 1e-8,
@@ -1876,6 +1888,7 @@ def _bbox_diag_len(Cc, Ss):
         rel_thresh=1e-6,
     )
     pre_overlaps: list["OverlapIntersection"] = []
+    pre_isolated: list["IsolatedIntersection"] = []
     stats = IntersectionStats(cells=0, pruned=0, overlap_traces=0, pruned_by=[])
     if span_overlap is not None:
         if span_overlap.get("partial"):
@@ -1884,13 +1897,19 @@ def _bbox_diag_len(Cc, Ss):
             stats = IntersectionStats(
                 cells=0, pruned=0, overlap_traces=1, pruned_by=["span_confirm_partial"]
             )
+        elif span_overlap.get('isolated'):
+            stats = IntersectionStats(cells=1, pruned=0, overlap_traces=0, pruned_by=["span_isolated_confirm"])
+            del span_overlap['isolated']
+            pre_isolated.append(span_overlap)
+            return IntersectionResult(isolated=pre_isolated, overlaps=[], stats=stats)
+
         else:
             stats = IntersectionStats(cells=1, pruned=0, overlap_traces=1, pruned_by=["span_confirm"])
             return IntersectionResult(isolated=[], overlaps=[span_overlap], stats=stats)
-    #isolated,overlaps=curve_surface_boundary_intersect(C,S,rational=rational,atol=atol,sv_thresh=sv_thresh)
-    #print(isolated)
-    #print(overlaps)
-    isolated: list["IsolatedIntersection"] = []
+    # isolated,overlaps=curve_surface_boundary_intersect(C,S,rational=rational,atol=atol,sv_thresh=sv_thresh)
+    # print(isolated)
+    # print(overlaps)
+    isolated: list["IsolatedIntersection"] = pre_isolated
     overlaps: list["OverlapIntersection"] = pre_overlaps
 
     atol_sq = atol * atol
@@ -2376,7 +2395,7 @@ def make_rich_table(inter1, inter2, *args, **kwargs):
            [ -6.52097953,  26.04573929,   0.        ],
            [  5.81667121,  -9.0084582 ,  17.47802526]])
     s=time.perf_counter()
-    res=bezier_curve_surface_intersect_certified(crv1, surf1,rational=False)
+    res=bez_csx(crv1, surf1, rational=False)
 
     gt={'isolated':[{
         "t":0.223246,
@@ -2391,7 +2410,7 @@ def make_rich_table(inter1, inter2, *args, **kwargs):
 
     crv2 = np.array([[-8.089928672303788, -9.35078823406258, 15.194922104381954, 1.0], [-5.720443423501492, -2.6029262297235642, 10.744432459609845, 0.7071067811865476], [-13.75962333043464, -3.6810935759317296, 15.194922104381954, 1.0]])
     s=time.perf_counter()
-    res=bezier_curve_surface_intersect_certified(crv2, surf2,rational=True)
+    res=bez_csx(crv2, surf2, rational=True)
     print("case2", time.perf_counter() - s)
     gt={
         "isolated": [],
@@ -2411,7 +2430,7 @@ def make_rich_table(inter1, inter2, *args, **kwargs):
     crv3 = np.array([[-13.704782222797013, 20.157097244338125, 4.980117410327855, 1.0], [-13.780432085162484, 14.253220150508163, 3.521474791948015, 0.7071067811865476], [-19.48847395019814, 14.373405516937003, 4.980117410327855, 1.0]]
                     )
     s=time.perf_counter()
-    res=bezier_curve_surface_intersect_certified(crv3, surf3,rational=True)
+    res=bez_csx(crv3, surf3, rational=True)
     print("case3", time.perf_counter() - s)
     gt={
       "isolated": [],
@@ -2448,7 +2467,7 @@ def make_rich_table(inter1, inter2, *args, **kwargs):
     )
 
     s=time.perf_counter()
-    res = bezier_curve_surface_intersect_certified(crv4, surf3, rational=True)
+    res = bez_csx(crv4, surf3, rational=True)
     print("case4", time.perf_counter() - s)
     gt={
       "isolated": [