Refactor NURBS intersection logic: enhance interpolation with interpolate_nurbs_curve, fix camera initialization in viewer, and improve SSX result handling with new example script.

sth-v · sth-v · commit 7b47c4d9bd59 · 2026-02-09T08:05:15.000+03:00
diff --git a/examples/ssx/common_helpers.py b/examples/ssx/common_helpers.py
@@ -125,7 +125,7 @@ def draw_ssx(
 
     bb = AABB.from_points(s1.control_points.reshape(-1, 3)).merge(AABB.from_points(s2.control_points.reshape(-1, 3)))
     if viewer is None:
-        viewer = Viewer(camera=OrbitCamera(target=bb.centroid(), distance=np.linalg.norm(bb.diag())*2, near=1.0))
+        viewer = Viewer(camera=OrbitCamera(target=bb.centroid(), distance=np.linalg.norm(bb.diag())**2, near=1.0))
 
     viewer.add_nurbs_surface(s1, color=surf1_material.wires_material.color, surface_color=surf1_material.color, u_count=surf1_material.wires_material.u_count, v_count=surf1_material.wires_material.v_count,show_edges=surf1_material.show_wires,show_isocurves=surf1_material.show_wires,)
     viewer.add_nurbs_surface(
diff --git a/examples/ssx/nurbs_nurbs_intersection_11.py b/examples/ssx/nurbs_nurbs_intersection_11.py
@@ -0,0 +1,108 @@
+"""
+ 
+ This example demonstrates the intersection of two NURBS spheres.
+"""
+import time
+
+from mmcore.construction import cylinder_surface_2pt
+from mmcore.geom._nurbs_eval import _tuple_to_nurbs, NURBSSurfaceTuple
+from mmcore.geom._nurbs_transform import transform_nurbs
+from mmcore.geom.bvh.lbvh import AABB
+from mmcore.geom.nurbs import NURBSSurface
+from mmcore.construction import nurbs_surface
+from mmcore.numeric.intersection.ssx import ssx
+
+# Creating intersection objects
+import numpy as np
+
+import numpy as np
+from mmcore.geom._nurbs_eval import NURBSSurfaceTuple
+
+
+s1 = NURBSSurfaceTuple(
+    order_u=4,
+    order_v=4,
+    knot_u=np.array([   0.00000000,    0.00000000,    0.00000000,    0.00000000,
+            584.56714755, 1169.13429511, 1169.13429511, 1169.13429511,
+           1169.13429511]),
+    knot_v=np.array([  0.00000000,   0.00000000,   0.00000000,   0.00000000,
+           675.00000000, 675.00000000, 675.00000000, 675.00000000]),
+    control_points=np.array([[[2488.56941006, 1499.08843254,    0.00000000],
+            [2394.60779838, 1603.70965174,    0.00000000],
+            [2394.60779838, 1675.09483487,    0.00000000],
+            [2488.56941006, 1779.71605406,    0.00000000]],
+
+           [[2582.53102174, 1394.46721335,    0.00000000],
+            [2550.34379877, 1515.98031620,  152.12435146],
+            [2550.34379877, 1762.82417041,  152.12435146],
+            [2582.53102174, 1884.33727326,    0.00000000]],
+
+           [[2754.10779838, 1284.40224330,    0.00000000],
+            [2754.10779838, 1466.26757810,  257.50000000],
+            [2754.10779838, 1812.53690851,  257.50000000],
+            [2754.10779838, 1994.40224330,    0.00000000]],
+
+           [[2925.68457502, 1394.46721335,    0.00000000],
+            [2957.87179798, 1515.98031620,  152.12435146],
+            [2957.87179798, 1762.82417041,  152.12435146],
+            [2925.68457502, 1884.33727326,    0.00000000]],
+
+           [[3019.64618670, 1499.08843254,    0.00000000],
+            [3113.60779838, 1603.70965174,    0.00000000],
+            [3113.60779838, 1675.09483487,    0.00000000],
+            [3019.64618670, 1779.71605406,    0.00000000]]]),
+    weights=np.array([[1.00000000, 1.00000000, 1.00000000, 1.00000000],
+           [1.00000000, 1.00000000, 1.00000000, 1.00000000],
+           [1.00000000, 1.00000000, 1.00000000, 1.00000000],
+           [1.00000000, 1.00000000, 1.00000000, 1.00000000],
+           [1.00000000, 1.00000000, 1.00000000, 1.00000000]])
+)
+import numpy as np
+from mmcore.geom._nurbs_eval import NURBSSurfaceTuple
+
+
+s2 = NURBSSurfaceTuple(
+    order_u=2,
+    order_v=2,
+    knot_u=np.array([  0.00000000,   0.00000000, 854.10934693, 854.10934693]),
+    knot_v=np.array([   0.00000000,    0.00000000, 1119.84031048, 1119.84031048]),
+    control_points=np.array([[[2372.10906381, 1470.47565841,  170.90189109],
+            [2614.15055984, 2054.83790061,  228.66096126]],
+
+           [[2953.02855072, 1248.46099847,  -17.28737813],
+            [3195.07004676, 1832.82324067,   40.47169204]]]),
+    weights=np.array([[1.00000000, 1.00000000],
+           [1.00000000, 1.00000000]])
+)
+
+
+s1,s2=s2,s1
+import logging
+from examples.ssx.common_helpers import parse_args, save_pkl, draw_ssx, VIEWER_INSTALLED, CurveMaterial, ControlNetMaterial, PointMaterial
+args = parse_args()
+logging.basicConfig(level=getattr(logging, args.loglevel, logging.INFO))
+from mmcore.numeric.intersection.ssx import nurbs_ssx
+
+s = time.time()
+result = nurbs_ssx(s1, s2, atol=0.001, angle_tol=0.052)
+
+print(f"intersection computed at: {time.time() - s} sec.")
+print(len(result[0]), "branch(s)")
+print(len(result[1]), "pts(s)")
+
+if args.save_pkl or args.pkl_path is not None:
+    path = save_pkl(s1, s2, result, fp=args.pkl_path)
+    print(path.absolute().as_posix())
+
+RENDER = args.viewer
+
+if RENDER:
+    inter_curves_mat = CurveMaterial(
+        (0.0, 1.0, 0.5, 1.0),
+        show_control_net=args.show_inter_cpts,
+        control_net_material=ControlNetMaterial((0.0, 1.0, 0.5, 0.7), control_point_material=PointMaterial((0.0, 1.0, 0.5, 0.4), size=8)),
+    )
+
+    viewer = draw_ssx(s1, s2, result, intersection_curves_material=inter_curves_mat)
+
+    viewer.run()
diff --git a/mmcore/geom/_nurbs_interp.py b/mmcore/geom/_nurbs_interp.py
@@ -310,9 +310,18 @@ def fair_interpolate_curve(points, degree, lambda_reg=1e-3)->tuple[NDArray[float
     control_points = np.vstack([d0, x, d_end])
 
     return control_points, knot_vector
+def interpolate_nurbs_curve(points, degree,  use_centripetal=False,rational=False,**kwargs):
 
+    cp,kv=interpolate_curve(points, degree, use_centripetal=use_centripetal,**kwargs)
+    cp=np.array(cp)
+    if rational:
+        cp,w=from_homogeneous_1d(cp)
+    else:
+        w=np.ones_like(cp[:,0])
+    return NURBSCurveTuple(order=degree+1, knot=np.array(kv), control_points=cp, weights=w)
 from mmcore.geom.nurbs import ders_basis_funs as _ders_basis_funs
 from mmcore.geom.nurbs import find_span as _find_span
+from mmcore.geom._nurbs_eval import from_homogeneous_1d,to_homogeneous_1d
 import numpy as np
 from scipy.special import comb, factorial
 from typing import NamedTuple, List, Optional
diff --git a/mmcore/numeric/intersection/ssx/_ssx4.py b/mmcore/numeric/intersection/ssx/_ssx4.py
@@ -10,7 +10,7 @@
 
 from mmcore.construction import nurbs_curve
 from mmcore.geom._nurbs_eval import evaluate_nurbs_curve, evaluate_nurbs_surface, NURBSSurfaceTuple
-from mmcore.geom._nurbs_interp import interpolate_curve
+from mmcore.geom._nurbs_interp import interpolate_curve, interpolate_nurbs_curve
 
 import functools
 import pickle
@@ -1822,7 +1822,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=nurbs_curve(np.array(_points),degree=1)
+    branch.curve_xyz=interpolate_nurbs_curve(np.array(_points),branch.curve.degree)
     #branch.curve_xyz=remove_knots_after_merge(crv,interior,atol)
 
 def compute_branch_curves(branch: SSXBranch, surface1:NURBSSurfaceTuple, surface2:NURBSSurfaceTuple,**kwargs):
diff --git a/mmcore/numeric/intersection/ssx/trace_inter_segm.py b/mmcore/numeric/intersection/ssx/trace_inter_segm.py
@@ -19,7 +19,9 @@
 import numpy as np
 from numpy.typing import NDArray
 
+from mmcore.geom._nurbs_interp import interpolate_curve
 from mmcore.geom._nurbs_knots import link_curves,remove_knot,remove_knot_curve_max
+from mmcore.geom._nurbs_interp import interpolate_nurbs_curve
 from mmcore.numeric.intersection.ssx.refine import refine_intersection_point
 from mmcore.numeric.sbern import bern_to_nurbs_bezier
 
@@ -518,10 +520,10 @@ def compute(self):
 
     def compute_normals_tangent_and_curvature(self):
         if self.s_eval['s1'].get('N') is None:
-            self.s_eval['s1']['N']=evaluate_normal2(   self.s_eval['s1']['Su'],self.s_eval['s1']['Sv'],self.s_eval['s1']['Suu'],self.s_eval['s1']['Suv'],self.s_eval['s1']['Svv'],self.s_eval['s1']['Svv'])
+            self.s_eval['s1']['N']=evaluate_normal2(   self.s_eval['s1']['Su'],self.s_eval['s1']['Sv'],self.s_eval['s1']['Suu'],self.s_eval['s1']['Suv'],self.s_eval['s1']['Svv'])
             self.s_eval['s1']['N']/=np.linalg.norm(self.s_eval['s1']['N'])
         if self.s_eval['s2'].get('N') is None:
-            self.s_eval['s2']['N']=evaluate_normal2(   self.s_eval['s2']['Su'],self.s_eval['s2']['Sv'],self.s_eval['s2']['Suu'],self.s_eval['s2']['Suv'],self.s_eval['s2']['Svv'],self.s_eval['s2']['Svv'])
+            self.s_eval['s2']['N']=evaluate_normal2(   self.s_eval['s2']['Su'],self.s_eval['s2']['Sv'],self.s_eval['s2']['Suu'],self.s_eval['s2']['Suv'],self.s_eval['s2']['Svv'])
             self.s_eval['s2']['N']/=np.linalg.norm(self.s_eval['s2']['N'])
 
         self.s_eval['T']=np.cross(self.s_eval['s1']['N'],
@@ -573,12 +575,16 @@ def __init__(self,start:InterpNode,end:InterpNode):
         end: InterpNode
         start.curve_stuv: NDArray = None
         end.curve_xyz: NDArray = None
+        self._curve=None
 
     def stuv_mid(self) -> NDArray[np.float64]: ...
 
     def xyz_mid(self)->NDArray[np.float64]:...
 
-    def build_curves(self):...
+    def build_curves(self):
+
+        self._curve=hinterp(np.array((self.start.s_eval['s1']['S'],self.end.s_eval['s1']['S'])),np.array((self.start.s_eval['T'],self.end.s_eval['T'])))
+
 
     @property
     def next_segm(self):
@@ -649,11 +655,13 @@ def adaptive_refine_bruteforce(
 
     node_a = InterpNode(  {'stuv':stuv_a,'s1':xyz_s,'s2':xyz_s2,"T":None,"pN":None,"K":None},prev=None, next=None)
     node_b = InterpNode({'stuv':stuv_b,'s1':xyz_e,'s2':xyz_e2,"T":None,"pN":None,"K":None}, prev=node_a, next=None)
+    #max_h=np.linalg.norm(stuv_a-stuv_b)/8
 
     node_a.next = node_b
     node_a.compute()
     node_b.compute()
     stack = [(node_a, node_b, depth, np.inf)]
+
     while stack:
         stuv_start, stuv_end, cur_depth,delta = stack.pop()
         if cur_depth >= fit_max_depth:
@@ -685,9 +693,13 @@ def adaptive_refine_bruteforce(
         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']))
+
+
+        print(s,dT)
         if dT<=s:
-            continue
+                continue
         delta=np.linalg.norm(p_eval['S']-xyz1)
 
 
@@ -698,8 +710,9 @@ def adaptive_refine_bruteforce(
             stack.append((new_node,stuv_end,cur_depth+1,delta))
         else:
             continue
-
-    return nurbs_curve(np.array(list(node_a)),1)
+    ll=list(node_a)
+    print(len(ll))
+    return interpolate_nurbs_curve(np.array(ll),degree=min(len(ll)-1,3))
 
 
 def _from_homogeneous(H: NDArray[np.float64]) -> tuple[NDArray[np.float64], NDArray[np.float64]]:

Original file line number	Diff line number	Diff line change
`@@ -10,7 +10,7 @@`
`10`	`10`
`11`	`11`	`from mmcore.construction import nurbs_curve`
`12`	`12`	`from mmcore.geom._nurbs_eval import evaluate_nurbs_curve, evaluate_nurbs_surface, NURBSSurfaceTuple`
`13`		`-from mmcore.geom._nurbs_interp import interpolate_curve`
	`13`	`+from mmcore.geom._nurbs_interp import interpolate_curve, interpolate_nurbs_curve`
`14`	`14`
`15`	`15`	`import functools`
`16`	`16`	`import pickle`
`@@ -1822,7 +1822,7 @@ def _eval_curve(t):`
`1822`	`1822`	`#`
`1823`	`1823`	`#`
`1824`	`1824`	`# crvs.append(bern_to_nurbs_bezier(np.array([ps,ps+dsh,pe-dsh,pe]),interval=(ts,te),rational=False))`
`1825`		`- branch.curve_xyz=nurbs_curve(np.array(_points),degree=1)`
	`1825`	`+ branch.curve_xyz=interpolate_nurbs_curve(np.array(_points),branch.curve.degree)`
`1826`	`1826`	`#branch.curve_xyz=remove_knots_after_merge(crv,interior,atol)`
`1827`	`1827`
`1828`	`1828`	`def compute_branch_curves(branch: SSXBranch, surface1:NURBSSurfaceTuple, surface2:NURBSSurfaceTuple,**kwargs):`