Refactor SSX intersection logic: simplify curve evaluation, improve parameter initialization, add overlap deduplication, and optimize branch interpolation with centripetal method.

Author: sth-v

mmcore/numeric/intersection/ssx/_ssx4.py

@@ -1792,8 +1792,8 @@ def compute_branch_curves_hermite(branch: SSXBranch, surface1:NURBSSurfaceTuple,
     def _eval_curve(t):
         stuv = evaluate_nurbs_curve(branch.curve, t, d_order=0)["C"]
         se1=evaluate_nurbs_surface(surface1,stuv[0],stuv[1],d_order=0)
-        se2=evaluate_nurbs_surface(surface2,stuv[2],stuv[3],d_order=0)
-        return (se2['S']-se1['S'])*0.5+se1['S']
+
+        return se1['S']
 
     eval_curve_ders=fdm(_eval_curve)
     _points = []
@@ -1802,16 +1802,15 @@ def _eval_curve(t):
     params=_greville_abscissae(branch.curve.knot, branch.curve.degree)
     # params=np.unique(branch.curve.knot)
     #p_cur=0
-
+    params=branch.curve.control_points
     for s in params:
         #d1=eval_curve_ders(s)
         #d1/=np.linalg.norm(d1)
 
         #_ders.append(d1)
-        _points.append(_eval_curve(s))
 
-        # p_cur+=np.dot(d1   ,d1)
 
+        _points.append(evaluate_nurbs_surface(surface1, s[0], s[1], d_order=0)['S'])
 
 
     #from more_itertools import pairwise
@@ -1822,7 +1821,7 @@ def _eval_curve(t):
     #
     #
     #    crvs.append(bern_to_nurbs_bezier(np.array([ps,ps+ds*h,pe-ds*h,pe]),interval=(ts,te),rational=False))
-    branch.curve_xyz=interpolate_nurbs_curve(np.array(_points),branch.curve.degree)
+    branch.curve_xyz=interpolate_nurbs_curve(np.array(_points),degree=branch.curve.degree,use_centripetal=True,method='lu',remove_duplicates=True,tol=atol)
     #branch.curve_xyz=remove_knots_after_merge(crv,interior,atol)
 
 def compute_branch_curves(branch: SSXBranch, surface1:NURBSSurfaceTuple, surface2:NURBSSurfaceTuple,**kwargs):
@@ -1840,6 +1839,7 @@ def compute_branch_curves(branch: SSXBranch, surface1:NURBSSurfaceTuple, surface
 
 def compute_point_xyz(point:SSXPoint, surface1:NURBSSurfaceTuple, surface2:NURBSSurfaceTuple,**kwargs):
     point.xyz=evaluate_nurbs_surface(surface1,point.stuv[0],point.stuv[1],d_order=0)['S']
+
 def bez_ssx(H1: NDArray[np.float64], H2: NDArray[np.float64], *, atol: float = 0.001, angle_tol: float = 0.052,
             tol: float = 1e-8, gjk_max_iter: int = 64, gm_eps: float = 1e-5, gm_tol: float = 1e-8, max_depth: int = 24,
             magic_start_depth: int = 6, flat_angle: float = 0.015, march_samples: int = 8,
@@ -1864,7 +1864,9 @@ def nurbs_ssx(surf1, surf2, *, atol: float = 0.001, angle_tol: float = 0.052, to
               deflate_hard_case: bool = True) -> tuple[list[SSXBranch], list[SSXPoint]]:
     s1 = surf1
     s2 =surf2
-
+    if s1.control_points.shape ==s2.control_points.shape:
+        if np.allclose(s1.control_points,s2.control_points):
+            return
     s1d = decompose_surface(s1)
     s2d = decompose_surface(s2)
 

mmcore/numeric/intersection/ssx/trace_inter_segm.py

@@ -689,28 +689,39 @@ def adaptive_refine_bruteforce(
 
         stuv_start.next=new_node
         stuv_end.prev=new_node
+        xyzs = evaluate_nurbs_surface(s1, stuv_start.s_eval['stuv'][0], stuv_start.s_eval['stuv'][1], d_order=0)['S']
+        xyze = evaluate_nurbs_surface(s1, stuv_end.s_eval['stuv'][0], stuv_end.s_eval['stuv'][1], d_order=0)['S']
 
         new_node.compute_normals_tangent_and_curvature()
-        T = (new_node.s_eval['s1']["S"] - stuv_start.s_eval['s1']["S"]) / 2
-        dT=np.linalg.norm(T)
 
         s = np.sqrt(8.0 * spt / np.linalg.norm(new_node.s_eval['K']))
 
+        delta = np.linalg.norm(p_eval['S'] - xyz1)
+        s=np.clip(s,spt,np.linalg.norm(xyzs-xyze))
 
+        if delta< spt:
+
+            continue
+        xyzs = evaluate_nurbs_surface(s1, stuv_start.s_eval['stuv'] [0], stuv_start.s_eval['stuv'] [1], d_order=0)['S']
+        xyze = evaluate_nurbs_surface(s1, stuv_end.s_eval['stuv'] [0], stuv_end.s_eval['stuv'] [1], d_order=0)['S']
+        dT1=np.linalg.norm(xyzs-p_eval['S'])
+        dT2=np.linalg.norm(xyze - p_eval['S'])
+        if dT1 <s or dT2<s:
+
+            continue
+
+        stack.append((stuv_start, new_node, cur_depth + 1, delta))
+
+        stack.append((new_node, stuv_end, cur_depth + 1, delta))
         #print(s,dT)
-        if dT<=s:
-                continue
-        delta=np.linalg.norm(p_eval['S']-xyz1)
 
 
 
-        if delta>= spt:
 
-            stack.append((stuv_start,new_node,cur_depth+1,delta))
-            stack.append((new_node,stuv_end,cur_depth+1,delta))
-        else:
-            continue
+
+
     ll=list(node_a)
+
     #print(len(ll))
     return interpolate_nurbs_curve(np.array(ll),degree=min(len(ll)-1,3))