From 47e039997100033050a267bbf97a489c5ff800c0 Mon Sep 17 00:00:00 2001
From: adam-urbanczyk <13981538+adam-urbanczyk@users.noreply.github.com>
Date: Mon, 10 Mar 2025 08:09:55 +0100
Subject: [PATCH 01/55] Adding fig

---
 cadquery/fig.py | 90 +++++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 90 insertions(+)
 create mode 100644 cadquery/fig.py

diff --git a/cadquery/fig.py b/cadquery/fig.py
new file mode 100644
index 000000000..0c7843a7f
--- /dev/null
+++ b/cadquery/fig.py
@@ -0,0 +1,90 @@
+from trame.app import get_server, Server
+from trame.widgets import html, vtk as vtk_widgets, client
+from trame.ui.html import DivLayout
+
+from cadquery import Shape
+from cadquery.vis import style
+
+from vtkmodules.vtkRenderingCore import (
+    vtkRenderer,
+    vtkRenderWindow,
+    vtkRenderWindowInteractor,
+    vtkActor,
+)
+
+FULL_SCREEN = "position:absolute; left:0; top:0; width:100vw; height:100vh;"
+
+
+class Figure:
+
+    server: Server
+    win: vtkRenderWindow
+    ren: vtkRenderer
+    shapes: dict[Shape, tuple[vtkActor, ...]]
+    actors: list[vtkActor]
+
+    def __init__(self, port: int):
+
+        # vtk boilerplate
+        renderer = vtkRenderer()
+        win = vtkRenderWindow()
+        win.AddRenderer(renderer)
+        win.OffScreenRenderingOn()
+
+        inter = vtkRenderWindowInteractor()
+        inter.SetRenderWindow(win)
+        inter.GetInteractorStyle().SetCurrentStyleToTrackballCamera()
+
+        self.win = win
+        self.ren = renderer
+
+        self.shapes = {}
+        self.actors = []
+
+        # server
+        server = get_server(123)
+        server.client_type = "vue3"
+
+        # layout
+        with DivLayout(server):
+            client.Style("body { margin: 0; }")
+
+            with html.Div(style=FULL_SCREEN):
+                vtk_widgets.VtkRemoteView(
+                    win, interactive_ratio=1, interactive_quality=100
+                )
+
+        server.state.flush()
+        server.start(thread=True, exec_mode="task", port=port, open_browser=True)
+
+        self.server = server
+
+    def show(self, s: Shape | vtkActor, *args, **kwargs):
+
+        if isinstance(s, Shape):
+            actors = style(s, *args, **kwargs)
+            self.shapes[s] = actors
+
+            for a in actors:
+                self.ren.AddActor(a)
+        else:
+            self.actors.append(s)
+            self.ren.AddActor(s)
+
+        self.ren.ResetCamera()
+
+    def hide(self, s: Shape | vtkActor):
+
+        if isinstance(s, Shape):
+            actors = self.shapes[s]
+
+            for a in actors:
+                self.ren.RemoveActor(a)
+
+            del self.shapes[s]
+
+        else:
+            self.actors.remove(s)
+            self.ren.RemoveActor(s)
+
+        self.ren.ResetCamera()

From a918517c3d99d4fec910c1f641df70f163a850ee Mon Sep 17 00:00:00 2001
From: adam-urbanczyk <13981538+adam-urbanczyk@users.noreply.github.com>
Date: Tue, 11 Mar 2025 08:10:08 +0100
Subject: [PATCH 02/55] Mypy fixes

---
 cadquery/fig.py               | 69 +++++++++++++++++++++++++++--------
 cadquery/occ_impl/assembly.py | 11 +++---
 cadquery/vis.py               | 64 +++++++++++++++++---------------
 mypy.ini                      |  3 ++
 4 files changed, 96 insertions(+), 51 deletions(-)

diff --git a/cadquery/fig.py b/cadquery/fig.py
index 0c7843a7f..759777e17 100644
--- a/cadquery/fig.py
+++ b/cadquery/fig.py
@@ -1,3 +1,7 @@
+import asyncio
+
+from typing import Any, cast
+
 from trame.app import get_server, Server
 from trame.widgets import html, vtk as vtk_widgets, client
 from trame.ui.html import DivLayout
@@ -10,8 +14,15 @@
     vtkRenderWindow,
     vtkRenderWindowInteractor,
     vtkActor,
+    vtkProp3D,
 )
 
+
+from vtkmodules.vtkInteractionWidgets import vtkOrientationMarkerWidget
+from vtkmodules.vtkRenderingAnnotation import vtkAxesActor
+
+from vtkmodules.vtkInteractionStyle import vtkInteractorStyleTrackballCamera
+
 FULL_SCREEN = "position:absolute; left:0; top:0; width:100vw; height:100vh;"
 
 
@@ -20,21 +31,44 @@ class Figure:
     server: Server
     win: vtkRenderWindow
     ren: vtkRenderer
-    shapes: dict[Shape, tuple[vtkActor, ...]]
-    actors: list[vtkActor]
+    view: vtk_widgets.VtkRemoteView
+    shapes: dict[Shape, vtkProp3D]
+    actors: list[vtkProp3D]
+    loop: Any
 
     def __init__(self, port: int):
 
+        self.loop = asyncio.new_event_loop()
+        asyncio.set_event_loop(self.loop)
+
         # vtk boilerplate
         renderer = vtkRenderer()
         win = vtkRenderWindow()
+        w, h = win.GetScreenSize()
+        win.SetSize(w, h)
         win.AddRenderer(renderer)
         win.OffScreenRenderingOn()
 
         inter = vtkRenderWindowInteractor()
+        inter.SetInteractorStyle(vtkInteractorStyleTrackballCamera())
         inter.SetRenderWindow(win)
-        inter.GetInteractorStyle().SetCurrentStyleToTrackballCamera()
 
+        # axes
+        axes = vtkAxesActor()
+        axes.SetDragable(0)
+
+        orient_widget = vtkOrientationMarkerWidget()
+
+        orient_widget.SetOrientationMarker(axes)
+        orient_widget.SetViewport(0.9, 0.0, 1.0, 0.2)
+        orient_widget.SetZoom(1.1)
+        orient_widget.SetInteractor(inter)
+        orient_widget.SetCurrentRenderer(renderer)
+        orient_widget.EnabledOn()
+        orient_widget.InteractiveOff()
+
+        self.axes = axes
+        self.orient_widget = orient_widget
         self.win = win
         self.ren = renderer
 
@@ -50,7 +84,7 @@ def __init__(self, port: int):
             client.Style("body { margin: 0; }")
 
             with html.Div(style=FULL_SCREEN):
-                vtk_widgets.VtkRemoteView(
+                self.view = vtk_widgets.VtkRemoteView(
                     win, interactive_ratio=1, interactive_quality=100
                 )
 
@@ -59,27 +93,29 @@ def __init__(self, port: int):
 
         self.server = server
 
-    def show(self, s: Shape | vtkActor, *args, **kwargs):
+    def show(self, s: Shape | vtkActor | list[vtkProp3D], *args, **kwargs):
 
         if isinstance(s, Shape):
-            actors = style(s, *args, **kwargs)
-            self.shapes[s] = actors
-
-            for a in actors:
-                self.ren.AddActor(a)
-        else:
+            actor = style(s, *args, **kwargs)[0]
+            self.shapes[s] = actor
+            self.ren.AddActor(actor)
+        elif isinstance(s, vtkActor):
             self.actors.append(s)
             self.ren.AddActor(s)
+        else:
+            self.actors.extend(s)
+
+            for el in s:
+                self.ren.AddActor(el)
 
         self.ren.ResetCamera()
+        self.view.update()
 
-    def hide(self, s: Shape | vtkActor):
+    def clear(self, s: Shape | vtkActor):
 
         if isinstance(s, Shape):
-            actors = self.shapes[s]
-
-            for a in actors:
-                self.ren.RemoveActor(a)
+            actor = self.shapes[s]
+            self.ren.RemoveActor(actor)
 
             del self.shapes[s]
 
@@ -88,3 +124,4 @@ def hide(self, s: Shape | vtkActor):
             self.ren.RemoveActor(s)
 
         self.ren.ResetCamera()
+        self.view.update()
diff --git a/cadquery/occ_impl/assembly.py b/cadquery/occ_impl/assembly.py
index 24a9518fa..2cdb0e3a7 100644
--- a/cadquery/occ_impl/assembly.py
+++ b/cadquery/occ_impl/assembly.py
@@ -31,7 +31,7 @@
     vtkActor,
     vtkPolyDataMapper as vtkMapper,
     vtkRenderer,
-    vtkAssembly,
+    vtkProp3D,
 )
 
 from vtkmodules.vtkFiltersExtraction import vtkExtractCellsByType
@@ -295,9 +295,9 @@ def toVTKAssy(
     linewidth: float = 2,
     tolerance: float = 1e-3,
     angularTolerance: float = 0.1,
-) -> vtkAssembly:
+) -> List[vtkProp3D]:
 
-    rv = vtkAssembly()
+    rv: List[vtkProp3D] = []
 
     for shape, _, loc, col_ in assy:
 
@@ -338,7 +338,7 @@ def toVTKAssy(
         actor.GetProperty().SetColor(*col[:3])
         actor.GetProperty().SetOpacity(col[3])
 
-        rv.AddPart(actor)
+        rv.append(actor)
 
         mapper = vtkMapper()
         mapper.AddInputDataObject(data_edges)
@@ -350,9 +350,8 @@ def toVTKAssy(
         actor.GetProperty().SetLineWidth(linewidth)
         actor.SetVisibility(edges)
         actor.GetProperty().SetColor(*edgecolor[:3])
-        actor.GetProperty().SetLineWidth(edgecolor[3])
 
-        rv.AddPart(actor)
+        rv.append(actor)
 
     return rv
 
diff --git a/cadquery/vis.py b/cadquery/vis.py
index 5fbedb12f..c6c1f42d5 100644
--- a/cadquery/vis.py
+++ b/cadquery/vis.py
@@ -145,12 +145,12 @@ def _to_vtk_pts(
     return rv
 
 
-def _to_vtk_axs(locs: List[Location], scale: float = 0.1) -> vtkAssembly:
+def _to_vtk_axs(locs: List[Location], scale: float = 0.1) -> List[vtkProp3D]:
     """
     Convert Locations to vtkActor.
     """
 
-    rv = vtkAssembly()
+    rv: List[vtkProp3D] = []
 
     for l in locs:
         trans, rot = _loc2vtk(l)
@@ -161,7 +161,7 @@ def _to_vtk_axs(locs: List[Location], scale: float = 0.1) -> vtkAssembly:
         ax.SetOrientation(*rot)
         ax.SetScale(scale)
 
-        rv.AddPart(ax)
+        rv.append(ax)
 
     return rv
 
@@ -174,7 +174,7 @@ def _to_vtk_shapes(
     linewidth: float = 2,
     alpha: float = 1,
     tolerance: float = 1e-3,
-) -> vtkAssembly:
+) -> List[vtkProp3D]:
     """
     Convert Shapes to vtkAssembly.
     """
@@ -279,21 +279,18 @@ def ctrlPts(
     return rv
 
 
-def _iterate_actors(obj: Union[vtkProp3D, vtkActor, vtkAssembly]) -> Iterable[vtkActor]:
+def _iterate_actors(
+    obj: Union[vtkProp3D, vtkActor, List[vtkProp3D]]
+) -> Iterable[vtkActor]:
     """
     Iterate over vtkActors, other props are ignored.
     """
     if isinstance(obj, vtkActor):
         yield obj
-    elif isinstance(obj, vtkAssembly):
-        coll = vtkPropCollection()
-        obj.GetActors(coll)
-
-        coll.InitTraversal()
-        for i in range(0, coll.GetNumberOfItems()):
-            prop = coll.GetNextProp()
-            if isinstance(prop, vtkActor):
-                yield prop
+    elif isinstance(obj, list):
+        for el in obj:
+            if isinstance(el, vtkActor):
+                yield el
 
 
 def style(
@@ -313,7 +310,7 @@ def style(
     meshcolor: str = "lightgrey",
     vertexcolor: str = "cyan",
     **kwargs,
-) -> Union[vtkActor, vtkAssembly]:
+) -> List[vtkProp3D]:
     """
     Apply styling to CQ objects. To be used in conjunction with show.
     """
@@ -343,7 +340,7 @@ def _apply_color(actor):
     shapes, vecs, locs, actors = _split_showables([obj,])
 
     # convert to a prop
-    rv: Union[vtkActor, vtkAssembly]
+    rv: Union[vtkActor, List[vtkProp3D]]
 
     if shapes:
         rv = _to_vtk_shapes(
@@ -361,19 +358,22 @@ def _apply_color(actor):
             _apply_style(a)
 
     elif vecs:
-        rv = _to_vtk_pts(vecs)
-        _apply_style(rv)
-        _apply_color(rv)
+        tmp = _to_vtk_pts(vecs)
+        _apply_style(tmp)
+        _apply_color(tmp)
+        rv = [tmp]
+
     elif locs:
         rv = _to_vtk_axs(locs, scale=scale)
+
     else:
-        rv = vtkAssembly()
+        rv = []
 
         for p in actors:
             for a in _iterate_actors(p):
                 _apply_style(a)
                 _apply_color(a)
-                rv.AddPart(a)
+                rv.append(a)
 
     return rv
 
@@ -400,7 +400,7 @@ def show(
     gradient: bool = True,
     xpos: Union[int, float] = 0,
     ypos: Union[int, float] = 0,
-):
+) -> vtkRenderWindow:
     """
     Show CQ objects using VTK. This functions optionally allows to make screenshots.
     """
@@ -417,7 +417,9 @@ def show(
 
     # assy+renderer
     renderer = vtkRenderer()
-    renderer.AddActor(toVTKAssy(assy, tolerance=tolerance))
+
+    for act in toVTKAssy(assy, tolerance=tolerance):
+        renderer.AddActor(act)
 
     # VTK window boilerplate
     win = vtkRenderWindow()
@@ -455,12 +457,12 @@ def show(
     # construct an axes indicator
     axes = vtkAxesActor()
     axes.SetDragable(0)
+    axes.SetAxisLabels(0)
+    # tp = axes.GetXAxisCaptionActor2D().GetCaptionTextProperty()
+    # tp.SetColor(0, 0, 0)
 
-    tp = axes.GetXAxisCaptionActor2D().GetCaptionTextProperty()
-    tp.SetColor(0, 0, 0)
-
-    axes.GetYAxisCaptionActor2D().GetCaptionTextProperty().ShallowCopy(tp)
-    axes.GetZAxisCaptionActor2D().GetCaptionTextProperty().ShallowCopy(tp)
+    # axes.GetYAxisCaptionActor2D().GetCaptionTextProperty().ShallowCopy(tp)
+    # axes.GetZAxisCaptionActor2D().GetCaptionTextProperty().ShallowCopy(tp)
 
     # add to an orientation widget
     if trihedron:
@@ -483,7 +485,9 @@ def show(
 
     # add pts and locs
     renderer.AddActor(pts)
-    renderer.AddActor(axs)
+
+    for ax in axs:
+        renderer.AddActor(ax)
 
     # add other vtk actors
     for p in props:
@@ -538,6 +542,8 @@ def show(
     if interact:
         inter.Start()
 
+    return win
+
 
 # alias
 show_object = show
diff --git a/mypy.ini b/mypy.ini
index 97bbf2b5d..7bc958faf 100644
--- a/mypy.ini
+++ b/mypy.ini
@@ -37,3 +37,6 @@ ignore_missing_imports = True
 [mypy-casadi.*]
 ignore_missing_imports = True
 
+[mypy-trame.*]
+ignore_missing_imports = True
+

From 3b0a8709021a1eddb116b6ecd0cbd28bae5f3c90 Mon Sep 17 00:00:00 2001
From: adam-urbanczyk <13981538+adam-urbanczyk@users.noreply.github.com>
Date: Tue, 11 Mar 2025 19:28:03 +0100
Subject: [PATCH 03/55] Fix tests

---
 tests/test_vis.py | 22 +++++++++++++---------
 1 file changed, 13 insertions(+), 9 deletions(-)

diff --git a/tests/test_vis.py b/tests/test_vis.py
index 513ec2671..e22ff2eba 100644
--- a/tests/test_vis.py
+++ b/tests/test_vis.py
@@ -10,6 +10,7 @@
     vtkWindowToImageFilter,
     vtkActor,
     vtkAssembly,
+    vtkProp3D,
 )
 from vtkmodules.vtkRenderingAnnotation import vtkAnnotatedCubeActor
 from vtkmodules.vtkIOImage import vtkPNGWriter
@@ -17,6 +18,9 @@
 from pytest import fixture, raises
 from path import Path
 
+from typish import instance_of
+from typing import List
+
 
 @fixture(scope="module")
 def tmpdir(tmp_path_factory):
@@ -178,39 +182,39 @@ def test_style(wp, assy):
 
     # Shape
     act = style(t, color="red", alpha=0.5, tubes=True, spheres=True)
-    assert isinstance(act, (vtkActor, vtkAssembly))
+    assert instance_of(act, List[vtkProp3D])
 
     # Assy
     act = style(assy, color="red", alpha=0.5, tubes=True, spheres=True)
-    assert isinstance(act, (vtkActor, vtkAssembly))
+    assert instance_of(act, List[vtkProp3D])
 
     # Workplane
     act = style(wp, color="red", alpha=0.5, tubes=True, spheres=True)
-    assert isinstance(act, (vtkActor, vtkAssembly))
+    assert instance_of(act, List[vtkProp3D])
 
     # Shape
     act = style(e)
-    assert isinstance(act, (vtkActor, vtkAssembly))
+    assert instance_of(act, List[vtkProp3D])
 
     # Sketch
     act = style(Sketch().circle(1))
-    assert isinstance(act, (vtkActor, vtkAssembly))
+    assert instance_of(act, List[vtkProp3D])
 
     # list[Vector]
     act = style(pts)
-    assert isinstance(act, (vtkActor, vtkAssembly))
+    assert instance_of(act, List[vtkProp3D])
 
     # list[Location]
     act = style(locs)
-    assert isinstance(act, (vtkActor, vtkAssembly))
+    assert instance_of(act, List[vtkProp3D])
 
     # vtkAssembly
     act = style(style(t))
-    assert isinstance(act, (vtkActor, vtkAssembly))
+    assert instance_of(act, List[vtkProp3D])
 
     # vtkActor
     act = style(ctrlPts(e.toNURBS()))
-    assert isinstance(act, (vtkActor, vtkAssembly))
+    assert instance_of(act, List[vtkProp3D])
 
 
 def test_camera_position(wp, patch_vtk):

From fef77aa023a426e2cdb84511e3e8318c9e05ae7d Mon Sep 17 00:00:00 2001
From: adam-urbanczyk <13981538+adam-urbanczyk@users.noreply.github.com>
Date: Tue, 11 Mar 2025 19:29:43 +0100
Subject: [PATCH 04/55] Run the asyncio loop in a thread

---
 cadquery/fig.py | 106 +++++++++++++++++++++++++++++++++---------------
 1 file changed, 74 insertions(+), 32 deletions(-)

diff --git a/cadquery/fig.py b/cadquery/fig.py
index 759777e17..46423a947 100644
--- a/cadquery/fig.py
+++ b/cadquery/fig.py
@@ -1,6 +1,10 @@
-import asyncio
-
-from typing import Any, cast
+from asyncio import (
+    new_event_loop,
+    set_event_loop,
+    run_coroutine_threadsafe,
+    AbstractEventLoop,
+)
+from threading import Thread
 
 from trame.app import get_server, Server
 from trame.widgets import html, vtk as vtk_widgets, client
@@ -32,14 +36,15 @@ class Figure:
     win: vtkRenderWindow
     ren: vtkRenderer
     view: vtk_widgets.VtkRemoteView
-    shapes: dict[Shape, vtkProp3D]
+    shapes: dict[Shape, list[vtkProp3D]]
     actors: list[vtkProp3D]
-    loop: Any
+    loop: AbstractEventLoop
+    thread: Thread
 
     def __init__(self, port: int):
 
-        self.loop = asyncio.new_event_loop()
-        asyncio.set_event_loop(self.loop)
+        self.loop = new_event_loop()
+        set_event_loop(self.loop)
 
         # vtk boilerplate
         renderer = vtkRenderer()
@@ -76,7 +81,7 @@ def __init__(self, port: int):
         self.actors = []
 
         # server
-        server = get_server(123)
+        server = get_server("CQ")
         server.client_type = "vue3"
 
         # layout
@@ -89,39 +94,76 @@ def __init__(self, port: int):
                 )
 
         server.state.flush()
-        server.start(thread=True, exec_mode="task", port=port, open_browser=True)
+        coro = server.start(
+            thread=True, exec_mode="coroutine", port=port, open_browser=True
+        )
 
         self.server = server
+        self.loop = new_event_loop()
+
+        def _run():
+            set_event_loop(self.loop)
+            self.loop.run_forever()
+
+        self.thread = Thread(target=_run, daemon=True)
+        self.thread.start()
+
+        run_coroutine_threadsafe(coro, self.loop)
+
+    def _run(self, coro):
+
+        run_coroutine_threadsafe(coro, self.loop)
 
     def show(self, s: Shape | vtkActor | list[vtkProp3D], *args, **kwargs):
+        async def _show():
+
+            if isinstance(s, Shape):
+                # do not show markers by default
+                if "markersize" not in kwargs:
+                    kwargs["markersize"] = 0
+
+                actors = style(s, *args, **kwargs)
+                self.shapes[s] = actors
+
+                for actor in actors:
+                    self.ren.AddActor(actor)
+
+            elif isinstance(s, vtkActor):
+                self.actors.append(s)
+                self.ren.AddActor(s)
+            else:
+                self.actors.extend(s)
+
+                for el in s:
+                    self.ren.AddActor(el)
+
+            self.ren.ResetCamera()
+            self.view.update()
 
-        if isinstance(s, Shape):
-            actor = style(s, *args, **kwargs)[0]
-            self.shapes[s] = actor
-            self.ren.AddActor(actor)
-        elif isinstance(s, vtkActor):
-            self.actors.append(s)
-            self.ren.AddActor(s)
-        else:
-            self.actors.extend(s)
+        self._run(_show())
 
-            for el in s:
-                self.ren.AddActor(el)
+    def clear(self, *shapes: Shape | vtkActor):
+        async def _clear():
 
-        self.ren.ResetCamera()
-        self.view.update()
+            if len(shapes) == 0:
+                for a in self.actors:
+                    self.ren.RemoveActor(a)
 
-    def clear(self, s: Shape | vtkActor):
+                for actors in self.shapes.values():
+                    for a in actors:
+                        self.ren.RemoveActor(a)
 
-        if isinstance(s, Shape):
-            actor = self.shapes[s]
-            self.ren.RemoveActor(actor)
+            for s in shapes:
+                if isinstance(s, Shape):
+                    for a in self.shapes[s]:
+                        self.ren.RemoveActor(a)
 
-            del self.shapes[s]
+                    del self.shapes[s]
+                else:
+                    self.actors.remove(s)
+                    self.ren.RemoveActor(s)
 
-        else:
-            self.actors.remove(s)
-            self.ren.RemoveActor(s)
+            self.ren.ResetCamera()
+            self.view.update()
 
-        self.ren.ResetCamera()
-        self.view.update()
+        self._run(_clear())

From 524884e7674e7d2a0e6fab60ca759913f6c69794 Mon Sep 17 00:00:00 2001
From: adam-urbanczyk <13981538+adam-urbanczyk@users.noreply.github.com>
Date: Wed, 12 Mar 2025 19:01:26 +0100
Subject: [PATCH 05/55] Revert some changes in vis

---
 cadquery/vis.py | 18 ++++++++----------
 1 file changed, 8 insertions(+), 10 deletions(-)

diff --git a/cadquery/vis.py b/cadquery/vis.py
index c6c1f42d5..47b522370 100644
--- a/cadquery/vis.py
+++ b/cadquery/vis.py
@@ -150,7 +150,7 @@ def _to_vtk_axs(locs: List[Location], scale: float = 0.1) -> List[vtkProp3D]:
     Convert Locations to vtkActor.
     """
 
-    rv: List[vtkProp3D] = []
+    rv = vtkAssembly()
 
     for l in locs:
         trans, rot = _loc2vtk(l)
@@ -161,9 +161,9 @@ def _to_vtk_axs(locs: List[Location], scale: float = 0.1) -> List[vtkProp3D]:
         ax.SetOrientation(*rot)
         ax.SetScale(scale)
 
-        rv.append(ax)
+        rv.AddPart(ax)
 
-    return rv
+    return [rv]
 
 
 def _to_vtk_shapes(
@@ -400,7 +400,7 @@ def show(
     gradient: bool = True,
     xpos: Union[int, float] = 0,
     ypos: Union[int, float] = 0,
-) -> vtkRenderWindow:
+):
     """
     Show CQ objects using VTK. This functions optionally allows to make screenshots.
     """
@@ -458,11 +458,11 @@ def show(
     axes = vtkAxesActor()
     axes.SetDragable(0)
     axes.SetAxisLabels(0)
-    # tp = axes.GetXAxisCaptionActor2D().GetCaptionTextProperty()
-    # tp.SetColor(0, 0, 0)
+    tp = axes.GetXAxisCaptionActor2D().GetCaptionTextProperty()
+    tp.SetColor(0, 0, 0)
 
-    # axes.GetYAxisCaptionActor2D().GetCaptionTextProperty().ShallowCopy(tp)
-    # axes.GetZAxisCaptionActor2D().GetCaptionTextProperty().ShallowCopy(tp)
+    axes.GetYAxisCaptionActor2D().GetCaptionTextProperty().ShallowCopy(tp)
+    axes.GetZAxisCaptionActor2D().GetCaptionTextProperty().ShallowCopy(tp)
 
     # add to an orientation widget
     if trihedron:
@@ -542,8 +542,6 @@ def show(
     if interact:
         inter.Start()
 
-    return win
-
 
 # alias
 show_object = show

From 3bb10364f0d0e6e0add28f302a33f0ed574063ff Mon Sep 17 00:00:00 2001
From: adam-urbanczyk <13981538+adam-urbanczyk@users.noreply.github.com>
Date: Wed, 12 Mar 2025 19:02:09 +0100
Subject: [PATCH 06/55] Smaller coros and bg color

---
 cadquery/fig.py | 54 ++++++++++++++++++++++++++-----------------------
 1 file changed, 29 insertions(+), 25 deletions(-)

diff --git a/cadquery/fig.py b/cadquery/fig.py
index 46423a947..5b121b4de 100644
--- a/cadquery/fig.py
+++ b/cadquery/fig.py
@@ -17,7 +17,6 @@
     vtkRenderer,
     vtkRenderWindow,
     vtkRenderWindowInteractor,
-    vtkActor,
     vtkProp3D,
 )
 
@@ -58,6 +57,10 @@ def __init__(self, port: int):
         inter.SetInteractorStyle(vtkInteractorStyleTrackballCamera())
         inter.SetRenderWindow(win)
 
+        # background
+        renderer.SetBackground(1, 1, 1)
+        renderer.GradientBackgroundOn()
+
         # axes
         axes = vtkAxesActor()
         axes.SetDragable(0)
@@ -94,55 +97,56 @@ def __init__(self, port: int):
                 )
 
         server.state.flush()
-        coro = server.start(
-            thread=True, exec_mode="coroutine", port=port, open_browser=True
-        )
 
         self.server = server
         self.loop = new_event_loop()
 
-        def _run():
+        def _run_loop():
             set_event_loop(self.loop)
             self.loop.run_forever()
 
-        self.thread = Thread(target=_run, daemon=True)
+        self.thread = Thread(target=_run_loop, daemon=True)
         self.thread.start()
 
-        run_coroutine_threadsafe(coro, self.loop)
+        coro = server.start(
+            thread=True, exec_mode="coroutine", port=port, open_browser=True
+        )
+
+        self._run(coro)
 
     def _run(self, coro):
 
         run_coroutine_threadsafe(coro, self.loop)
 
-    def show(self, s: Shape | vtkActor | list[vtkProp3D], *args, **kwargs):
-        async def _show():
+    def show(self, s: Shape | vtkProp3D | list[vtkProp3D], *args, **kwargs):
 
-            if isinstance(s, Shape):
-                # do not show markers by default
-                if "markersize" not in kwargs:
-                    kwargs["markersize"] = 0
+        if isinstance(s, Shape):
+            # do not show markers by default
+            if "markersize" not in kwargs:
+                kwargs["markersize"] = 0
 
-                actors = style(s, *args, **kwargs)
-                self.shapes[s] = actors
+            actors = style(s, *args, **kwargs)
+            self.shapes[s] = actors
 
-                for actor in actors:
-                    self.ren.AddActor(actor)
+            for actor in actors:
+                self.ren.AddActor(actor)
 
-            elif isinstance(s, vtkActor):
-                self.actors.append(s)
-                self.ren.AddActor(s)
-            else:
-                self.actors.extend(s)
+        elif isinstance(s, vtkProp3D):
+            self.actors.append(s)
+            self.ren.AddActor(s)
+        else:
+            self.actors.extend(s)
 
-                for el in s:
-                    self.ren.AddActor(el)
+            for el in s:
+                self.ren.AddActor(el)
 
+        async def _show():
             self.ren.ResetCamera()
             self.view.update()
 
         self._run(_show())
 
-    def clear(self, *shapes: Shape | vtkActor):
+    def clear(self, *shapes: Shape | vtkProp3D):
         async def _clear():
 
             if len(shapes) == 0:

From ca12abaa5b34d7b6e05fe85421488fec644fb8e7 Mon Sep 17 00:00:00 2001
From: adam-urbanczyk <13981538+adam-urbanczyk@users.noreply.github.com>
Date: Fri, 14 Mar 2025 22:11:50 +0100
Subject: [PATCH 07/55] Refactor into a singleton

---
 cadquery/fig.py | 98 ++++++++++++++++++++++++++++++++++++++-----------
 1 file changed, 77 insertions(+), 21 deletions(-)

diff --git a/cadquery/fig.py b/cadquery/fig.py
index 5b121b4de..d7528feff 100644
--- a/cadquery/fig.py
+++ b/cadquery/fig.py
@@ -5,13 +5,15 @@
     AbstractEventLoop,
 )
 from threading import Thread
+from itertools import chain
+from uuid import uuid1 as uuid
 
 from trame.app import get_server, Server
 from trame.widgets import html, vtk as vtk_widgets, client
 from trame.ui.html import DivLayout
 
-from cadquery import Shape
-from cadquery.vis import style
+from . import Shape
+from .vis import style, Showable, ShapeLike, _split_showables, _to_vtk_pts, _to_vtk_axs
 
 from vtkmodules.vtkRenderingCore import (
     vtkRenderer,
@@ -35,12 +37,26 @@ class Figure:
     win: vtkRenderWindow
     ren: vtkRenderer
     view: vtk_widgets.VtkRemoteView
-    shapes: dict[Shape, list[vtkProp3D]]
+    shapes: dict[ShapeLike, list[vtkProp3D]]
     actors: list[vtkProp3D]
     loop: AbstractEventLoop
     thread: Thread
+    empty: bool
 
-    def __init__(self, port: int):
+    _instance = None
+    _initialized: bool = False
+
+    def __new__(cls, *args, **kwargs):
+
+        if not cls._instance:
+            cls._instance = object.__new__(cls)
+
+        return cls._instance
+
+    def __init__(self, port: int = 18081):
+
+        if self._initialized:
+            return
 
         self.loop = new_event_loop()
         set_event_loop(self.loop)
@@ -84,7 +100,7 @@ def __init__(self, port: int):
         self.actors = []
 
         # server
-        server = get_server("CQ")
+        server = get_server("CQ-server")
         server.client_type = "vue3"
 
         # layout
@@ -112,33 +128,65 @@ def _run_loop():
             thread=True, exec_mode="coroutine", port=port, open_browser=True
         )
 
-        self._run(coro)
+        if coro:
+            self._run(coro)
+
+        # prevent reinitialization
+        self._initialized = True
+
+        # view is initialized as empty
+        self.empty = True
 
     def _run(self, coro):
 
         run_coroutine_threadsafe(coro, self.loop)
 
-    def show(self, s: Shape | vtkProp3D | list[vtkProp3D], *args, **kwargs):
+    def show(self, *showables: Showable | vtkProp3D | list[vtkProp3D], **kwargs):
+        """
+        Show objects.
+        """
+
+        # split objects
+        shapes, vecs, locs, props = _split_showables(showables)
+
+        pts = _to_vtk_pts(vecs)
+        axs = _to_vtk_axs(locs)
 
-        if isinstance(s, Shape):
+        for s in shapes:
             # do not show markers by default
             if "markersize" not in kwargs:
                 kwargs["markersize"] = 0
 
-            actors = style(s, *args, **kwargs)
+            actors = style(s, **kwargs)
             self.shapes[s] = actors
 
             for actor in actors:
                 self.ren.AddActor(actor)
 
-        elif isinstance(s, vtkProp3D):
-            self.actors.append(s)
-            self.ren.AddActor(s)
-        else:
-            self.actors.extend(s)
+        for prop in chain(props, axs):
+            self.actors.append(prop)
+            self.ren.AddActor(prop)
 
-            for el in s:
-                self.ren.AddActor(el)
+        if vecs:
+            self.actors.append(pts)
+            self.ren.AddActor(pts)
+
+        async def _show():
+            self.view.update()
+
+        self._run(_show())
+
+        # zoom to fit on 1st object added
+        if self.empty:
+            self.fit()
+            self.empty = False
+
+        return self
+
+    def fit(self):
+        """
+        Update view to fit all objects.
+        """
 
         async def _show():
             self.ren.ResetCamera()
@@ -146,16 +194,23 @@ async def _show():
 
         self._run(_show())
 
+        return self
+
     def clear(self, *shapes: Shape | vtkProp3D):
+        """
+        Clear specified objects. If no arguments are passed, clears all objects.
+        """
+
         async def _clear():
 
             if len(shapes) == 0:
-                for a in self.actors:
+                for a in self.ren.GetActors():
                     self.ren.RemoveActor(a)
 
-                for actors in self.shapes.values():
-                    for a in actors:
-                        self.ren.RemoveActor(a)
+                self.actors.clear()
+                self.shapes.clear()
+
+                self.empty = True
 
             for s in shapes:
                 if isinstance(s, Shape):
@@ -167,7 +222,8 @@ async def _clear():
                     self.actors.remove(s)
                     self.ren.RemoveActor(s)
 
-            self.ren.ResetCamera()
             self.view.update()
 
         self._run(_clear())
+
+        return self

From 31f8c2854abd1447a4b7ad8aeb5237a63b564fca Mon Sep 17 00:00:00 2001
From: adam-urbanczyk <13981538+adam-urbanczyk@users.noreply.github.com>
Date: Fri, 14 Mar 2025 22:42:51 +0100
Subject: [PATCH 08/55] Update deps

---
 conda/meta.yaml | 2 ++
 environment.yml | 2 ++
 setup.py        | 2 ++
 3 files changed, 6 insertions(+)

diff --git a/conda/meta.yaml b/conda/meta.yaml
index e3f6d3a28..2f976f844 100644
--- a/conda/meta.yaml
+++ b/conda/meta.yaml
@@ -26,6 +26,8 @@ requirements:
     - multimethod >=1.11,<2.0
     - casadi
     - typish
+    - trame
+    - trame-vtk
 
 test:
   requires:
diff --git a/environment.yml b/environment.yml
index 6d03a2359..3d5f4ef69 100644
--- a/environment.yml
+++ b/environment.yml
@@ -25,6 +25,8 @@ dependencies:
   - pathspec
   - click
   - appdirs
+  - trame
+  - trame-vtk
   - pip
   - pip:
     - --editable=.
diff --git a/setup.py b/setup.py
index 45442d3e1..c8f162df8 100644
--- a/setup.py
+++ b/setup.py
@@ -33,6 +33,8 @@
         "typish",
         "casadi",
         "path",
+        "trame",
+        "trame-vtk"
     ]
 
 

From 44fbb64b06b861d976e01047132126ecdaa437bc Mon Sep 17 00:00:00 2001
From: adam-urbanczyk <13981538+adam-urbanczyk@users.noreply.github.com>
Date: Sat, 15 Mar 2025 15:52:35 +0100
Subject: [PATCH 09/55] Implemented pop and fixed clear

---
 cadquery/fig.py | 46 ++++++++++++++++++++++++++++++++++++++++++----
 1 file changed, 42 insertions(+), 4 deletions(-)

diff --git a/cadquery/fig.py b/cadquery/fig.py
index d7528feff..f013dcb20 100644
--- a/cadquery/fig.py
+++ b/cadquery/fig.py
@@ -4,9 +4,9 @@
     run_coroutine_threadsafe,
     AbstractEventLoop,
 )
+from typing import Optional
 from threading import Thread
 from itertools import chain
-from uuid import uuid1 as uuid
 
 from trame.app import get_server, Server
 from trame.widgets import html, vtk as vtk_widgets, client
@@ -20,6 +20,7 @@
     vtkRenderWindow,
     vtkRenderWindowInteractor,
     vtkProp3D,
+    vtkActor,
 )
 
 
@@ -42,6 +43,9 @@ class Figure:
     loop: AbstractEventLoop
     thread: Thread
     empty: bool
+    last: Optional[
+        tuple[list[ShapeLike], list[vtkProp3D], Optional[vtkActor], list[vtkProp3D]]
+    ]
 
     _instance = None
     _initialized: bool = False
@@ -136,6 +140,7 @@ def _run_loop():
 
         # view is initialized as empty
         self.empty = True
+        self.last = None
 
     def _run(self, coro):
 
@@ -150,7 +155,9 @@ def show(self, *showables: Showable | vtkProp3D | list[vtkProp3D], **kwargs):
         shapes, vecs, locs, props = _split_showables(showables)
 
         pts = _to_vtk_pts(vecs)
-        axs = _to_vtk_axs(locs)
+        axs = _to_vtk_axs(
+            locs, **({"scale": kwargs["scale"]} if "scale" in kwargs else {})
+        )
 
         for s in shapes:
             # do not show markers by default
@@ -171,6 +178,9 @@ def show(self, *showables: Showable | vtkProp3D | list[vtkProp3D], **kwargs):
             self.actors.append(pts)
             self.ren.AddActor(pts)
 
+        # store to enable pop
+        self.last = (shapes, axs, pts if vecs else None, props)
+
         async def _show():
             self.view.update()
 
@@ -204,8 +214,7 @@ def clear(self, *shapes: Shape | vtkProp3D):
         async def _clear():
 
             if len(shapes) == 0:
-                for a in self.ren.GetActors():
-                    self.ren.RemoveActor(a)
+                self.ren.RemoveAllViewProps()
 
                 self.actors.clear()
                 self.shapes.clear()
@@ -224,6 +233,35 @@ async def _clear():
 
             self.view.update()
 
+        # reset last, bc we don't want to keep track of what was removed
+        self.last = None
         self._run(_clear())
 
         return self
+
+    def pop(self):
+        """
+        Clear the last showable.
+        """
+
+        async def _pop():
+
+            (shapes, axs, pts, props) = self.last
+
+            for s in shapes:
+                for act in self.shapes.pop(s):
+                    self.ren.RemoveActor(act)
+
+            for act in chain(axs, props):
+                self.ren.RemoveActor(act)
+                self.actors.remove(act)
+
+            if pts:
+                self.ren.RemoveActor(pts)
+                self.actors.remove(pts)
+
+            self.view.update()
+
+        self._run(_pop())
+
+        return self

From 125de0921f133d1dae122d0cb307d77f3bffe197 Mon Sep 17 00:00:00 2001
From: adam-urbanczyk <13981538+adam-urbanczyk@users.noreply.github.com>
Date: Sat, 15 Mar 2025 16:10:55 +0100
Subject: [PATCH 10/55] Styling fix

---
 cadquery/fig.py | 14 +++++++-------
 cadquery/vis.py |  9 ++++++++-
 2 files changed, 15 insertions(+), 8 deletions(-)

diff --git a/cadquery/fig.py b/cadquery/fig.py
index f013dcb20..9e5e625eb 100644
--- a/cadquery/fig.py
+++ b/cadquery/fig.py
@@ -44,7 +44,9 @@ class Figure:
     thread: Thread
     empty: bool
     last: Optional[
-        tuple[list[ShapeLike], list[vtkProp3D], Optional[vtkActor], list[vtkProp3D]]
+        tuple[
+            list[ShapeLike], list[vtkProp3D], Optional[list[vtkProp3D]], list[vtkProp3D]
+        ]
     ]
 
     _instance = None
@@ -154,10 +156,8 @@ def show(self, *showables: Showable | vtkProp3D | list[vtkProp3D], **kwargs):
         # split objects
         shapes, vecs, locs, props = _split_showables(showables)
 
-        pts = _to_vtk_pts(vecs)
-        axs = _to_vtk_axs(
-            locs, **({"scale": kwargs["scale"]} if "scale" in kwargs else {})
-        )
+        pts = style(vecs, **kwargs)
+        axs = style(locs, **kwargs)
 
         for s in shapes:
             # do not show markers by default
@@ -175,8 +175,8 @@ def show(self, *showables: Showable | vtkProp3D | list[vtkProp3D], **kwargs):
             self.ren.AddActor(prop)
 
         if vecs:
-            self.actors.append(pts)
-            self.ren.AddActor(pts)
+            self.actors.append(*pts)
+            self.ren.AddActor(*pts)
 
         # store to enable pop
         self.last = (shapes, axs, pts if vecs else None, props)
diff --git a/cadquery/vis.py b/cadquery/vis.py
index 47b522370..93e9194be 100644
--- a/cadquery/vis.py
+++ b/cadquery/vis.py
@@ -53,7 +53,14 @@
 
 ShapeLike = Union[Shape, Workplane, Assembly, Sketch, TopoDS_Shape]
 Showable = Union[
-    ShapeLike, List[ShapeLike], Vector, List[Vector], vtkProp3D, List[vtkProp3D]
+    ShapeLike,
+    List[ShapeLike],
+    Vector,
+    List[Vector],
+    vtkProp3D,
+    List[vtkProp3D],
+    Location,
+    List[Location],
 ]
 
 

From c3a4c634572981a9cbea5fe8fff0188906d59420 Mon Sep 17 00:00:00 2001
From: adam-urbanczyk <13981538+adam-urbanczyk@users.noreply.github.com>
Date: Sat, 15 Mar 2025 16:16:50 +0100
Subject: [PATCH 11/55] Fix pop

---
 cadquery/fig.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/cadquery/fig.py b/cadquery/fig.py
index 9e5e625eb..ea6991aa0 100644
--- a/cadquery/fig.py
+++ b/cadquery/fig.py
@@ -257,8 +257,8 @@ async def _pop():
                 self.actors.remove(act)
 
             if pts:
-                self.ren.RemoveActor(pts)
-                self.actors.remove(pts)
+                self.ren.RemoveActor(*pts)
+                self.actors.remove(*pts)
 
             self.view.update()
 

From 6f78783191eac0a8f9030916ab03b59e47cecf35 Mon Sep 17 00:00:00 2001
From: adam-urbanczyk <13981538+adam-urbanczyk@users.noreply.github.com>
Date: Sun, 16 Mar 2025 13:51:11 +0100
Subject: [PATCH 12/55] Add simple data

---
 setup.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/setup.py b/setup.py
index c8f162df8..80b5c3812 100644
--- a/setup.py
+++ b/setup.py
@@ -34,7 +34,7 @@
         "casadi",
         "path",
         "trame",
-        "trame-vtk"
+        "trame-vtk",
     ]
 
 

From 4a98f6abe979f3b5aeec20b9edc38fb23404d453 Mon Sep 17 00:00:00 2001
From: adam-urbanczyk <13981538+adam-urbanczyk@users.noreply.github.com>
Date: Mon, 17 Mar 2025 07:47:29 +0100
Subject: [PATCH 13/55] Get rid of vtk msgs

---
 cadquery/fig.py | 10 +++++++++-
 1 file changed, 9 insertions(+), 1 deletion(-)

diff --git a/cadquery/fig.py b/cadquery/fig.py
index ea6991aa0..3784608eb 100644
--- a/cadquery/fig.py
+++ b/cadquery/fig.py
@@ -7,6 +7,7 @@
 from typing import Optional
 from threading import Thread
 from itertools import chain
+from webbrowser import open_new_tab
 
 from trame.app import get_server, Server
 from trame.widgets import html, vtk as vtk_widgets, client
@@ -131,7 +132,11 @@ def _run_loop():
         self.thread.start()
 
         coro = server.start(
-            thread=True, exec_mode="coroutine", port=port, open_browser=True
+            thread=True,
+            exec_mode="coroutine",
+            port=port,
+            open_browser=False,
+            show_connection_info=False,
         )
 
         if coro:
@@ -144,6 +149,9 @@ def _run_loop():
         self.empty = True
         self.last = None
 
+        # open webbrowser
+        open_new_tab(f"http://localhost:{port}")
+
     def _run(self, coro):
 
         run_coroutine_threadsafe(coro, self.loop)

From b45fa6928b6e26d8a55b5d421048af689ea36104 Mon Sep 17 00:00:00 2001
From: adam-urbanczyk <13981538+adam-urbanczyk@users.noreply.github.com>
Date: Mon, 31 Mar 2025 07:40:09 +0200
Subject: [PATCH 14/55] Wait for clear

---
 cadquery/fig.py | 8 +++++---
 1 file changed, 5 insertions(+), 3 deletions(-)

diff --git a/cadquery/fig.py b/cadquery/fig.py
index 3784608eb..bbfa5502d 100644
--- a/cadquery/fig.py
+++ b/cadquery/fig.py
@@ -3,6 +3,7 @@
     set_event_loop,
     run_coroutine_threadsafe,
     AbstractEventLoop,
+    Future,
 )
 from typing import Optional
 from threading import Thread
@@ -152,9 +153,9 @@ def _run_loop():
         # open webbrowser
         open_new_tab(f"http://localhost:{port}")
 
-    def _run(self, coro):
+    def _run(self, coro) -> Future:
 
-        run_coroutine_threadsafe(coro, self.loop)
+        return run_coroutine_threadsafe(coro, self.loop)
 
     def show(self, *showables: Showable | vtkProp3D | list[vtkProp3D], **kwargs):
         """
@@ -243,7 +244,8 @@ async def _clear():
 
         # reset last, bc we don't want to keep track of what was removed
         self.last = None
-        self._run(_clear())
+        future = self._run(_clear())
+        future.result()
 
         return self
 

From 56c6cfc79e8f874a63248a55cfe660c241132e81 Mon Sep 17 00:00:00 2001
From: adam-urbanczyk <13981538+adam-urbanczyk@users.noreply.github.com>
Date: Mon, 31 Mar 2025 07:42:47 +0200
Subject: [PATCH 15/55] Display axis labels

---
 cadquery/vis.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/cadquery/vis.py b/cadquery/vis.py
index 93e9194be..fb644386d 100644
--- a/cadquery/vis.py
+++ b/cadquery/vis.py
@@ -464,7 +464,7 @@ def show(
     # construct an axes indicator
     axes = vtkAxesActor()
     axes.SetDragable(0)
-    axes.SetAxisLabels(0)
+
     tp = axes.GetXAxisCaptionActor2D().GetCaptionTextProperty()
     tp.SetColor(0, 0, 0)
 

From dd7803e479d2f84b5e97ea3743048ecdd1be57cc Mon Sep 17 00:00:00 2001
From: adam-urbanczyk <13981538+adam-urbanczyk@users.noreply.github.com>
Date: Thu, 12 Jun 2025 17:50:31 +0200
Subject: [PATCH 16/55] Adding smoke test for fig

---
 tests/conftest.py | 22 ++++++++++++++++++++++
 tests/test_fig.py | 44 ++++++++++++++++++++++++++++++++++++++++++++
 2 files changed, 66 insertions(+)
 create mode 100644 tests/conftest.py
 create mode 100644 tests/test_fig.py

diff --git a/tests/conftest.py b/tests/conftest.py
new file mode 100644
index 000000000..d30b369de
--- /dev/null
+++ b/tests/conftest.py
@@ -0,0 +1,22 @@
+import pytest
+
+
+def pytest_addoption(parser):
+    parser.addoption("--gui", action="store_true", default=False, help="run gui tests")
+
+
+def pytest_configure(config):
+    config.addinivalue_line("markers", "gui: mark gui test")
+
+
+def pytest_collection_modifyitems(config, items):
+
+    # run gui tests --gui option is proveded
+    if config.getoption("--gui"):
+        return
+
+    # skip gui tests otherwise
+    skip_gui = pytest.mark.skip(reason="need --gui option to run")
+    for item in items:
+        if "gui" in item.keywords:
+            item.add_marker(skip_gui)
diff --git a/tests/test_fig.py b/tests/test_fig.py
new file mode 100644
index 000000000..c9cc5bf37
--- /dev/null
+++ b/tests/test_fig.py
@@ -0,0 +1,44 @@
+from cadquery import Workplane, Assembly, Sketch, Vector, Location
+from cadquery.func import box
+from cadquery.vis import vtkAxesActor, ctrlPts
+from cadquery.fig import Figure
+
+from pytest import fixture, mark
+
+
+@fixture(scope="module")
+def fig():
+    return Figure()
+
+
+@mark.gui
+def test_fig(fig):
+
+    # showables
+    s = box(1, 1, 1)
+    wp = Workplane().box(1, 1, 1)
+    assy = Assembly().add(box(1, 1, 1))
+    sk = Sketch().rect(1, 1)
+    ctrl_pts = ctrlPts(sk.val().toNURBS())
+    v = Vector()
+    loc = Location()
+    act = vtkAxesActor()
+
+    # individual showables
+    fig.show(s, wp, assy, sk, ctrl_pts, v, loc, act)
+
+    # fit
+    fig.fit()
+
+    # clear
+    fig.clear()
+
+    # lists of showables
+    fig.show(s.Edges()).show([Vector(), Vector(0, 1)])
+
+    # displaying nonsense does not throw
+    fig.show("a").show(["a", 1234])
+
+    # pop
+    fig.show(s, color="red")
+    fig.pop()

From 34a499afa4044f66793b99cc24cc418cfca26cb8 Mon Sep 17 00:00:00 2001
From: adam-urbanczyk <13981538+adam-urbanczyk@users.noreply.github.com>
Date: Thu, 12 Jun 2025 17:52:15 +0200
Subject: [PATCH 17/55] Run gui tests in appveyor

---
 appveyor.yml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/appveyor.yml b/appveyor.yml
index 5364fb68a..6c1632e00 100644
--- a/appveyor.yml
+++ b/appveyor.yml
@@ -37,7 +37,7 @@ build: false
 test_script:
     - mamba run -n cadquery black . --diff --check
     - mamba run -n cadquery mypy cadquery
-    - mamba run -n cadquery pytest -v --cov
+    - mamba run -n cadquery pytest -v --gui --cov
 
 on_success:
     - mamba run -n cadquery codecov

From 7d3df36cb7b992301da4f197ba6ee4a3864aea1c Mon Sep 17 00:00:00 2001
From: adam-urbanczyk <13981538+adam-urbanczyk@users.noreply.github.com>
Date: Thu, 12 Jun 2025 18:15:20 +0200
Subject: [PATCH 18/55] mypy fix

---
 cadquery/fig.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/cadquery/fig.py b/cadquery/fig.py
index bbfa5502d..e9447ddf4 100644
--- a/cadquery/fig.py
+++ b/cadquery/fig.py
@@ -3,8 +3,8 @@
     set_event_loop,
     run_coroutine_threadsafe,
     AbstractEventLoop,
-    Future,
 )
+from concurrent.futures import Future
 from typing import Optional
 from threading import Thread
 from itertools import chain

From 30edd4a2ac5a51c2598a88e1573a3b3dc8712d61 Mon Sep 17 00:00:00 2001
From: adam-urbanczyk <13981538+adam-urbanczyk@users.noreply.github.com>
Date: Thu, 12 Jun 2025 18:32:05 +0200
Subject: [PATCH 19/55] Misc fixes

---
 cadquery/fig.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/cadquery/fig.py b/cadquery/fig.py
index e9447ddf4..5143794e5 100644
--- a/cadquery/fig.py
+++ b/cadquery/fig.py
@@ -10,19 +10,19 @@
 from itertools import chain
 from webbrowser import open_new_tab
 
-from trame.app import get_server, Server
+from trame.app import get_server
+from trame.app.core import Server
 from trame.widgets import html, vtk as vtk_widgets, client
 from trame.ui.html import DivLayout
 
 from . import Shape
-from .vis import style, Showable, ShapeLike, _split_showables, _to_vtk_pts, _to_vtk_axs
+from .vis import style, Showable, ShapeLike, _split_showables
 
 from vtkmodules.vtkRenderingCore import (
     vtkRenderer,
     vtkRenderWindow,
     vtkRenderWindowInteractor,
     vtkProp3D,
-    vtkActor,
 )
 
 

From 2ad988424266d9c86b7a50d3765b8e69256077d7 Mon Sep 17 00:00:00 2001
From: adam-urbanczyk <13981538+adam-urbanczyk@users.noreply.github.com>
Date: Thu, 12 Jun 2025 19:15:12 +0200
Subject: [PATCH 20/55] Test GUI only on win

---
 tests/test_fig.py | 19 ++++++++++++++++---
 1 file changed, 16 insertions(+), 3 deletions(-)

diff --git a/tests/test_fig.py b/tests/test_fig.py
index c9cc5bf37..6789a57c7 100644
--- a/tests/test_fig.py
+++ b/tests/test_fig.py
@@ -5,6 +5,8 @@
 
 from pytest import fixture, mark
 
+from sys import platform
+
 
 @fixture(scope="module")
 def fig():
@@ -12,6 +14,7 @@ def fig():
 
 
 @mark.gui
+@mark.skipif(platform != "win32", reason="CI with UI only works on win for now")
 def test_fig(fig):
 
     # showables
@@ -24,8 +27,10 @@ def test_fig(fig):
     loc = Location()
     act = vtkAxesActor()
 
+    showables = (s, wp, assy, sk, ctrl_pts, v, loc, act)
+
     # individual showables
-    fig.show(s, wp, assy, sk, ctrl_pts, v, loc, act)
+    fig.show(*showables)
 
     # fit
     fig.fit()
@@ -33,6 +38,9 @@ def test_fig(fig):
     # clear
     fig.clear()
 
+    # clear with an arg
+    fig.show(s).clear(s)
+
     # lists of showables
     fig.show(s.Edges()).show([Vector(), Vector(0, 1)])
 
@@ -40,5 +48,10 @@ def test_fig(fig):
     fig.show("a").show(["a", 1234])
 
     # pop
-    fig.show(s, color="red")
-    fig.pop()
+    for el in showables:
+        fig.show(el, color="red")
+        fig.pop()
+
+    # test singleton behavior of fig
+    fig2 = Figure()
+    assert fig is fig2

From d1c726a87ac6c13a4a5336827c15ab21bfd3a693 Mon Sep 17 00:00:00 2001
From: adam-urbanczyk <13981538+adam-urbanczyk@users.noreply.github.com>
Date: Thu, 12 Jun 2025 19:42:17 +0200
Subject: [PATCH 21/55] Coverage tweak

---
 tests/test_fig.py | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/tests/test_fig.py b/tests/test_fig.py
index 6789a57c7..3d810f840 100644
--- a/tests/test_fig.py
+++ b/tests/test_fig.py
@@ -39,7 +39,8 @@ def test_fig(fig):
     fig.clear()
 
     # clear with an arg
-    fig.show(s).clear(s)
+    for el in showables:
+        fig.show(el).clear(el)
 
     # lists of showables
     fig.show(s.Edges()).show([Vector(), Vector(0, 1)])

From 99af75fdd71e625ab9b0d8ef234ca6421083e486 Mon Sep 17 00:00:00 2001
From: adam-urbanczyk <13981538+adam-urbanczyk@users.noreply.github.com>
Date: Thu, 12 Jun 2025 21:26:44 +0200
Subject: [PATCH 22/55] Fix test

---
 cadquery/fig.py   | 2 +-
 tests/test_fig.py | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/cadquery/fig.py b/cadquery/fig.py
index 5143794e5..5a3298946 100644
--- a/cadquery/fig.py
+++ b/cadquery/fig.py
@@ -231,7 +231,7 @@ async def _clear():
                 self.empty = True
 
             for s in shapes:
-                if isinstance(s, Shape):
+                if isinstance(s, ShapeLike):
                     for a in self.shapes[s]:
                         self.ren.RemoveActor(a)
 
diff --git a/tests/test_fig.py b/tests/test_fig.py
index 3d810f840..ffa09264c 100644
--- a/tests/test_fig.py
+++ b/tests/test_fig.py
@@ -39,7 +39,7 @@ def test_fig(fig):
     fig.clear()
 
     # clear with an arg
-    for el in showables:
+    for el in (s, wp, assy, sk, ctrl_pts):
         fig.show(el).clear(el)
 
     # lists of showables

From f7bf1bf070ce421b11b5ff145b8eb7df3d0f954d Mon Sep 17 00:00:00 2001
From: adam-urbanczyk <13981538+adam-urbanczyk@users.noreply.github.com>
Date: Fri, 13 Jun 2025 08:05:07 +0200
Subject: [PATCH 23/55] Mypy fix

---
 cadquery/fig.py | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/cadquery/fig.py b/cadquery/fig.py
index 5a3298946..4d08201c9 100644
--- a/cadquery/fig.py
+++ b/cadquery/fig.py
@@ -10,6 +10,8 @@
 from itertools import chain
 from webbrowser import open_new_tab
 
+from typish import instance_of
+
 from trame.app import get_server
 from trame.app.core import Server
 from trame.widgets import html, vtk as vtk_widgets, client
@@ -231,7 +233,7 @@ async def _clear():
                 self.empty = True
 
             for s in shapes:
-                if isinstance(s, ShapeLike):
+                if instance_of(s, ShapeLike):
                     for a in self.shapes[s]:
                         self.ren.RemoveActor(a)
 

From 204df308ea090c553e01d589c050f58eb56c66e9 Mon Sep 17 00:00:00 2001
From: adam-urbanczyk <13981538+adam-urbanczyk@users.noreply.github.com>
Date: Mon, 16 Jun 2025 08:04:41 +0200
Subject: [PATCH 24/55] Change zoom reset behavior

---
 cadquery/fig.py | 2 --
 1 file changed, 2 deletions(-)

diff --git a/cadquery/fig.py b/cadquery/fig.py
index 4d08201c9..a68a43e59 100644
--- a/cadquery/fig.py
+++ b/cadquery/fig.py
@@ -230,8 +230,6 @@ async def _clear():
                 self.actors.clear()
                 self.shapes.clear()
 
-                self.empty = True
-
             for s in shapes:
                 if instance_of(s, ShapeLike):
                     for a in self.shapes[s]:

From 4d025d0db1c57394fc95d8af563ad91e26ed0e7e Mon Sep 17 00:00:00 2001
From: adam-urbanczyk <13981538+adam-urbanczyk@users.noreply.github.com>
Date: Fri, 27 Jun 2025 17:46:36 +0200
Subject: [PATCH 25/55] Initial commit of b-spline/nurbs tools

---
 cadquery/occ_impl/nurbs.py | 345 +++++++++++++++++++++++++++++++++++++
 conda/meta.yaml            |   1 +
 environment.yml            |   1 +
 mypy.ini                   |   3 +
 tests/test_nurbs.py        |  55 ++++++
 5 files changed, 405 insertions(+)
 create mode 100644 cadquery/occ_impl/nurbs.py
 create mode 100644 tests/test_nurbs.py

diff --git a/cadquery/occ_impl/nurbs.py b/cadquery/occ_impl/nurbs.py
new file mode 100644
index 000000000..8c97f5b5c
--- /dev/null
+++ b/cadquery/occ_impl/nurbs.py
@@ -0,0 +1,345 @@
+#%% imports
+from numba import njit, prange
+from numpy import linspace, array, empty_like, atleast_1d
+import numpy as np
+import math
+
+from numba import njit as _njit, prange
+from typing import NamedTuple, Optional
+from numpy.typing import NDArray
+
+njit = _njit(cache=False, error_model="numpy", fastmath=True, parallel=False)
+
+njiti = _njit(
+    cache=True, inline="always", error_model="numpy", fastmath=True, parallel=False
+)
+
+
+#%% vocabulary types
+
+Array = NDArray[np.float64]
+ArrayI = NDArray[np.int_]
+
+
+class COO(NamedTuple):
+    """
+    COO sparse matrix container.
+    """
+
+    i: ArrayI
+    j: ArrayI
+    v: Array
+
+
+#%% basis functions
+
+
+@njiti
+def nbFindSpan(
+    u: float,
+    order: int,
+    knots: Array,
+    low: Optional[int] = None,
+    high: Optional[int] = None,
+) -> int:
+    """
+    NURBS book A2.1 with modifications to handle periodic usecases.
+
+    Parameters
+    ----------
+    u : float
+        Parameter value.
+    order : int
+        Spline order.
+    knots : ndarray
+        Knot vectr.
+
+    Returns
+    -------
+    Span index.
+
+    """
+
+    if low is None:
+        low = order
+
+    if high is None:
+        high = knots.shape[0] - order - 1
+
+    mid = (low + high) // 2
+
+    if u >= knots[-1]:
+        return high - 1  # handle last span
+    elif u < knots[0]:
+        return low
+
+    while u < knots[mid] or u >= knots[mid + 1]:
+        if u < knots[mid]:
+            high = mid
+        else:
+            low = mid
+
+        mid = (low + high) // 2
+
+    return mid
+
+
+@njiti
+def nbBasis(i: int, u: float, order: int, knots: Array, out: Array):
+    """
+    NURBS book A2.2
+
+    Parameters
+    ----------
+    i : int
+        Span index.
+    u : float
+        Parameter value.
+    order : int
+        B-spline order.
+    knots : ndarray
+        Knot vectr.
+    out : ndarray
+        B-spline basis function values.
+
+    Returns
+    -------
+    None.
+
+    """
+
+    out[0] = 1.0
+
+    left = np.zeros_like(out)
+    right = np.zeros_like(out)
+
+    for j in range(1, order + 1):
+        left[j] = u - knots[i + 1 - j]
+        right[j] = knots[i + j] - u
+
+        saved = 0.0
+
+        for r in range(j):
+            temp = out[r] / (right[r + 1] + left[j - r])
+            out[r] = saved + right[r + 1] * temp
+            saved = left[j - r] * temp
+
+        out[j] = saved
+
+
+def nbBasisDer(i: int, u: float, order: int, dorder: int, knots: Array, out: Array):
+    """
+    NURBS book A2.3
+
+    Parameters
+    ----------
+    i : int
+        Span index.
+    u : float
+        Parameter value.
+    order : int
+        B-spline order.
+    dorder : int
+        Derivative order.
+    knots : ndarray
+        Knot vectr.
+    out : ndarray
+        B-spline basis function and derivative values.
+
+    Returns
+    -------
+    None.
+
+    """
+
+    ndu = np.zeros((order + 1, order + 1))
+
+    left = np.zeros(order + 1)
+    right = np.zeros(order + 1)
+
+    a = np.zeros((2, order + 1))
+
+    ndu[0, 0] = 1
+
+    for j in range(1, order + 1):
+        left[j] = u - knots[i + 1 - j]
+        right[j] = knots[i + j] - u
+
+        saved = 0.0
+
+        for r in range(j):
+            ndu[j, r] = right[r + 1] + left[j - r]
+            temp = ndu[r, j - 1] / ndu[j, r]
+
+            ndu[r, j] = saved + right[r + 1] * temp
+            saved = left[j - r] * temp
+
+        ndu[j, j] = saved
+
+    # store the basis functions
+    out[0, :] = ndu[:, order]
+
+    # calculate and store derivatives
+
+    # loop over basis functions
+    for r in range(order + 1):
+        s1 = 0
+        s2 = 1
+
+        a[0, 0] = 1
+
+        # loop over derivative orders
+        for k in range(1, dorder + 1):
+            d = 0.0
+            rk = r - k
+            pk = order - k
+
+            if r >= k:
+                a[s2, 0] = a[s1, 0] / ndu[pk + 1, rk]
+                d = a[s2, 0] * ndu[rk, pk]
+
+            if rk >= -1:
+                j1 = 1
+            else:
+                j1 = -rk
+
+            if r - 1 <= pk:
+                j2 = k - 1
+            else:
+                j2 = order - r
+
+            for j in range(j1, j2 + 1):
+                a[s2, j] = (a[s1, j] - a[s1, j - 1]) / ndu[pk + 1, rk + j]
+                d += a[s2, j] * ndu[rk + j, pk]
+
+            if r <= pk:
+                a[s2, k] = -a[s1, k - 1] / ndu[pk + 1, r]
+                d += a[s2, k] * ndu[r, pk]
+
+            # store the kth derivative of rth basis
+            out[k, r] = d
+
+            # switch
+            s1, s2 = s2, s1
+
+    # multiply recursively by the order
+    r = order
+
+    for k in range(1, dorder + 1):
+        out[k, :] *= r
+        r *= order - k
+
+
+@njit
+def designMatrix(u: Array, order: int, knots: Array) -> COO:
+    """
+    Create a sparse design matrix.
+    """
+
+    # number of param values
+    nu = np.size(u)
+
+    # chunck size
+    n = order + 1
+
+    # temp chunck storage
+    temp = np.zeros(n)
+
+    # initialize the empty matrix
+    rv = COO(
+        i=np.empty(n * nu, dtype=np.int64),
+        j=np.empty(n * nu, dtype=np.int64),
+        v=np.empty(n * nu),
+    )
+
+    # loop over param values
+    for i in range(nu):
+        ui = u[i]
+
+        # find the supporting span
+        span = nbFindSpan(ui, order, knots)
+
+        # evaluate non-zero functions
+        nbBasis(span, ui, order, knots, temp)
+
+        # update the matrix
+        rv.i[i * n : (i + 1) * n] = i
+        rv.j[i * n : (i + 1) * n] = span - order + np.arange(n)
+        rv.v[i * n : (i + 1) * n] = temp
+
+    return rv
+
+
+@njit
+def periodicDesignMatrix(u: Array, order: int, knots: Array) -> COO:
+    """
+    Create a sparse periodic design matrix.
+    """
+
+    # extend the knots
+    knots_ext = np.concat(
+        (knots[-order:-1] - knots[-1], knots, knots[-1] + knots[1:order])
+    )
+
+    # number of param values
+    nu = len(u)
+
+    # number of basis functions
+    nb = len(knots) - 1
+
+    # chunck size
+    n = order + 1
+
+    # temp chunck storage
+    temp = np.zeros(n)
+
+    # initialize the empty matrix
+    rv = COO(
+        i=np.empty(n * nu, dtype=np.int64),
+        j=np.empty(n * nu, dtype=np.int64),
+        v=np.empty(n * nu),
+    )
+
+    # loop over param values
+    for i in range(nu):
+        ui = u[i]
+
+        # find the supporting span
+        # span = np.clip(findSpan(ui, knots), None, nb - 1)  + order - 1
+        span = nbFindSpan(ui, order, knots, 0, nb) + order - 1
+
+        # evaluate non-zero functions
+        nbBasis(span, ui, order, knots_ext, temp)
+
+        # update the matrix
+        rv.i[i * n : (i + 1) * n] = i
+        rv.j[i * n : (i + 1) * n] = (
+            span - order + np.arange(n)
+        ) % nb  # NB: this is due to peridicity
+        rv.v[i * n : (i + 1) * n] = temp
+
+    return rv
+
+
+@njit
+def findSpan(v, knots):
+
+    return np.searchsorted(knots, v, "right") - 1
+
+
+@njit
+def findSpanLinear(v, knots):
+
+    for rv in range(len(knots)):
+        if knots[rv] <= v and knots[rv + 1] > v:
+            return rv
+
+    return -1
+
+
+@njit
+def periodicKnots(degree: int, n_pts: int):
+    rv = np.arange(0.0, n_pts + degree + 1, 1.0)
+    rv /= rv[-1]
+
+    return rv
diff --git a/conda/meta.yaml b/conda/meta.yaml
index e3f6d3a28..7fa50dc50 100644
--- a/conda/meta.yaml
+++ b/conda/meta.yaml
@@ -26,6 +26,7 @@ requirements:
     - multimethod >=1.11,<2.0
     - casadi
     - typish
+    - numba
 
 test:
   requires:
diff --git a/environment.yml b/environment.yml
index 6d03a2359..3e6dfefff 100644
--- a/environment.yml
+++ b/environment.yml
@@ -25,6 +25,7 @@ dependencies:
   - pathspec
   - click
   - appdirs
+  - numba
   - pip
   - pip:
     - --editable=.
diff --git a/mypy.ini b/mypy.ini
index 97bbf2b5d..a49c4d4c0 100644
--- a/mypy.ini
+++ b/mypy.ini
@@ -37,3 +37,6 @@ ignore_missing_imports = True
 [mypy-casadi.*]
 ignore_missing_imports = True
 
+[mypy-numba.*]
+ignore_missing_imports = True
+
diff --git a/tests/test_nurbs.py b/tests/test_nurbs.py
new file mode 100644
index 000000000..490d42d77
--- /dev/null
+++ b/tests/test_nurbs.py
@@ -0,0 +1,55 @@
+from cadquery.occ_impl.nurbs import (
+    designMatrix,
+    periodicDesignMatrix,
+    nbFindSpan,
+    nbBasis,
+    nbBasisDer,
+)
+
+import numpy as np
+import scipy.sparse as sp
+
+
+def test_periodic_dm():
+
+    knots = np.linspace(0, 1, 5)
+    params = np.linspace(0, 1, 100)
+    order = 3
+
+    res = periodicDesignMatrix(params, order, knots)
+
+    C = sp.coo_array((res.v, (res.i, res.j)))
+
+    assert C.shape[0] == len(params)
+    assert C.shape[1] == len(knots) - 1
+
+
+def test_dm():
+
+    knots = np.array([0, 0, 0, 0, 0.25, 0.5, 0.75, 1, 1, 1, 1])
+    params = np.linspace(0, 1, 100)
+    order = 3
+
+    res = designMatrix(params, order, knots)
+
+    C = sp.coo_array((res.v, (res.i, res.j)))
+
+    assert C.shape[0] == len(params)
+    assert C.shape[1] == len(knots) - order - 1
+
+
+def test_der():
+
+    knots = np.array([0, 0, 0, 0, 0.25, 0.5, 0.75, 1, 1, 1, 1])
+    params = np.linspace(0, 1, 100)
+    order = 3
+
+    out_der = np.zeros((order + 1, order + 1))
+    out = np.zeros(order + 1)
+
+    for p in params:
+        nbBasisDer(nbFindSpan(p, order, knots), p, order, order - 1, knots, out_der)
+        nbBasis(nbFindSpan(p, order, knots), p, order, knots, out)
+
+        # sanity check
+        assert np.allclose(out_der[0, :], out)

From cd52759083a05c5abb97e1ccc364d9885f048e37 Mon Sep 17 00:00:00 2001
From: adam-urbanczyk <13981538+adam-urbanczyk@users.noreply.github.com>
Date: Fri, 27 Jun 2025 17:58:07 +0200
Subject: [PATCH 26/55] Add scipy

---
 conda/meta.yaml | 1 +
 environment.yml | 1 +
 2 files changed, 2 insertions(+)

diff --git a/conda/meta.yaml b/conda/meta.yaml
index 7fa50dc50..70ffe6d94 100644
--- a/conda/meta.yaml
+++ b/conda/meta.yaml
@@ -27,6 +27,7 @@ requirements:
     - casadi
     - typish
     - numba
+    - scipy
 
 test:
   requires:
diff --git a/environment.yml b/environment.yml
index 3e6dfefff..01b92f9ec 100644
--- a/environment.yml
+++ b/environment.yml
@@ -26,6 +26,7 @@ dependencies:
   - click
   - appdirs
   - numba
+  - scipy
   - pip
   - pip:
     - --editable=.

From 048cef5c1d7becf6c11a07321a8338e11f04e4c0 Mon Sep 17 00:00:00 2001
From: adam-urbanczyk <13981538+adam-urbanczyk@users.noreply.github.com>
Date: Fri, 27 Jun 2025 18:36:26 +0200
Subject: [PATCH 27/55] Add derivative matrices

---
 cadquery/occ_impl/nurbs.py | 104 ++++++++++++++++++++++++++++++++++++-
 tests/test_nurbs.py        |   3 +-
 2 files changed, 105 insertions(+), 2 deletions(-)

diff --git a/cadquery/occ_impl/nurbs.py b/cadquery/occ_impl/nurbs.py
index 8c97f5b5c..4d488742c 100644
--- a/cadquery/occ_impl/nurbs.py
+++ b/cadquery/occ_impl/nurbs.py
@@ -127,6 +127,7 @@ def nbBasis(i: int, u: float, order: int, knots: Array, out: Array):
         out[j] = saved
 
 
+@njiti
 def nbBasisDer(i: int, u: float, order: int, dorder: int, knots: Array, out: Array):
     """
     NURBS book A2.3
@@ -305,7 +306,6 @@ def periodicDesignMatrix(u: Array, order: int, knots: Array) -> COO:
         ui = u[i]
 
         # find the supporting span
-        # span = np.clip(findSpan(ui, knots), None, nb - 1)  + order - 1
         span = nbFindSpan(ui, order, knots, 0, nb) + order - 1
 
         # evaluate non-zero functions
@@ -321,6 +321,108 @@ def periodicDesignMatrix(u: Array, order: int, knots: Array) -> COO:
     return rv
 
 
+@njit
+def derMatrix(u: Array, order: int, dorder: int, knots: Array) -> list[COO]:
+    """
+    Create a sparse design matrix and corresponding derivative matrices.
+    """
+
+    # number of param values
+    nu = np.size(u)
+
+    # chunck size
+    n = order + 1
+
+    # temp chunck storage
+    temp = np.zeros((n, n))
+
+    # initialize the empty matrix
+    rv = []
+
+    for _ in range(dorder + 1):
+        rv.append(
+            COO(
+                i=np.empty(n * nu, dtype=np.int64),
+                j=np.empty(n * nu, dtype=np.int64),
+                v=np.empty(n * nu),
+            )
+        )
+
+    # loop over param values
+    for i in range(nu):
+        ui = u[i]
+
+        # find the supporting span
+        span = nbFindSpan(ui, order, knots)
+
+        # evaluate non-zero functions
+        nbBasisDer(span, ui, order, dorder, knots, temp)
+
+        # update the matrices
+        for di in range(dorder + 1):
+            rv[di].i[i * n : (i + 1) * n] = i
+            rv[di].j[i * n : (i + 1) * n] = span - order + np.arange(n)
+            rv[di].v[i * n : (i + 1) * n] = temp[di, :]
+
+    return rv
+
+
+@njit
+def periodicDerMatrix(u: Array, order: int, dorder: int, knots: Array) -> list[COO]:
+    """
+    Create a sparse periodic design matrix and corresponding derivative matrices.
+    """
+
+    # extend the knots
+    knots_ext = np.concat(
+        (knots[-order:-1] - knots[-1], knots, knots[-1] + knots[1:order])
+    )
+
+    # number of param values
+    nu = len(u)
+
+    # number of basis functions
+    nb = len(knots) - 1
+
+    # chunck size
+    n = order + 1
+
+    # temp chunck storage
+    temp = np.zeros((n, n))
+
+    # initialize the empty matrix
+    rv = []
+
+    for _ in range(dorder + 1):
+        rv.append(
+            COO(
+                i=np.empty(n * nu, dtype=np.int64),
+                j=np.empty(n * nu, dtype=np.int64),
+                v=np.empty(n * nu),
+            )
+        )
+
+    # loop over param values
+    for i in range(nu):
+        ui = u[i]
+
+        # find the supporting span
+        span = nbFindSpan(ui, order, knots, 0, nb) + order - 1
+
+        # evaluate non-zero functions
+        nbBasisDer(span, ui, order, dorder, knots_ext, temp)
+
+        # update the matrices
+        for di in range(dorder + 1):
+            rv[di].i[i * n : (i + 1) * n] = i
+            rv[di].j[i * n : (i + 1) * n] = (
+                span - order + np.arange(n)
+            ) % nb  # NB: this is due to peridicity
+            rv[di].v[i * n : (i + 1) * n] = temp[di, :]
+
+    return rv
+
+
 @njit
 def findSpan(v, knots):
 
diff --git a/tests/test_nurbs.py b/tests/test_nurbs.py
index 490d42d77..6854995c1 100644
--- a/tests/test_nurbs.py
+++ b/tests/test_nurbs.py
@@ -48,7 +48,8 @@ def test_der():
     out = np.zeros(order + 1)
 
     for p in params:
-        nbBasisDer(nbFindSpan(p, order, knots), p, order, order - 1, knots, out_der)
+        nbBasisDer(nbFindSpan(p, order, knots), p,
+                   order, order - 1, knots, out_der)
         nbBasis(nbFindSpan(p, order, knots), p, order, knots, out)
 
         # sanity check

From 97fef3f72a83d93e525bb37c97452c07bdca2aee Mon Sep 17 00:00:00 2001
From: adam-urbanczyk <13981538+adam-urbanczyk@users.noreply.github.com>
Date: Fri, 27 Jun 2025 21:05:12 +0200
Subject: [PATCH 28/55] Black fix

---
 tests/test_nurbs.py | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/tests/test_nurbs.py b/tests/test_nurbs.py
index 6854995c1..490d42d77 100644
--- a/tests/test_nurbs.py
+++ b/tests/test_nurbs.py
@@ -48,8 +48,7 @@ def test_der():
     out = np.zeros(order + 1)
 
     for p in params:
-        nbBasisDer(nbFindSpan(p, order, knots), p,
-                   order, order - 1, knots, out_der)
+        nbBasisDer(nbFindSpan(p, order, knots), p, order, order - 1, knots, out_der)
         nbBasis(nbFindSpan(p, order, knots), p, order, knots, out)
 
         # sanity check

From be7d8db823f625220b8d14767159f605807ae628 Mon Sep 17 00:00:00 2001
From: adam-urbanczyk <13981538+adam-urbanczyk@users.noreply.github.com>
Date: Fri, 4 Jul 2025 18:02:02 +0200
Subject: [PATCH 29/55] Curve, Surface, lofting

---
 cadquery/occ_impl/nurbs.py | 479 ++++++++++++++++++++++++++++++++++++-
 tests/test_nurbs.py        | 158 ++++++++++++
 2 files changed, 635 insertions(+), 2 deletions(-)

diff --git a/cadquery/occ_impl/nurbs.py b/cadquery/occ_impl/nurbs.py
index 4d488742c..d9f07a92f 100644
--- a/cadquery/occ_impl/nurbs.py
+++ b/cadquery/occ_impl/nurbs.py
@@ -1,12 +1,28 @@
 #%% imports
-from numba import njit, prange
-from numpy import linspace, array, empty_like, atleast_1d
 import numpy as np
+import scipy.sparse as sp
 import math
 
 from numba import njit as _njit, prange
+
 from typing import NamedTuple, Optional
+
 from numpy.typing import NDArray
+from numpy import linspace, array, empty_like, atleast_1d
+
+from casadi import ldl, ldl_solve
+
+from OCP.Geom import Geom_BSplineCurve, Geom_BSplineSurface
+from OCP.TColgp import TColgp_Array1OfPnt, TColgp_Array2OfPnt
+from OCP.TColStd import (
+    TColStd_Array1OfInteger,
+    TColStd_Array1OfReal,
+    TColStd_Array2OfReal,
+)
+from OCP.gp import gp_Pnt
+from OCP.BRepBuilderAPI import BRepBuilderAPI_MakeEdge, BRepBuilderAPI_MakeFace
+
+from .shapes import Face, Edge
 
 njit = _njit(cache=False, error_model="numpy", fastmath=True, parallel=False)
 
@@ -15,6 +31,54 @@
 )
 
 
+#%% internal helpers
+
+
+def _colPtsArray(pts: NDArray):
+
+    rv = TColgp_Array1OfPnt(1, pts.shape[0])
+
+    for i, p in enumerate(pts):
+        rv.SetValue(i + 1, gp_Pnt(*p))
+
+    return rv
+
+
+def _colPtsArray2(pts: NDArray) -> TColStd_Array2OfReal:
+
+    assert pts.ndim == 3
+
+    nu, nv, _ = pts.shape
+
+    rv = TColgp_Array2OfPnt(1, len(pts), 1, len(pts[0]))
+
+    for i, row in enumerate(pts):
+        for j, pt in enumerate(row):
+            rv.SetValue(i + 1, j + 1, gp_Pnt(*pt))
+
+    return rv
+
+
+def _colRealArray(knots: NDArray):
+
+    rv = TColStd_Array1OfReal(1, len(knots))
+
+    for i, el in enumerate(knots):
+        rv.SetValue(i + 1, el)
+
+    return rv
+
+
+def _colIntArray(knots: NDArray):
+
+    rv = TColStd_Array1OfInteger(1, len(knots))
+
+    for i, el in enumerate(knots):
+        rv.SetValue(i + 1, el)
+
+    return rv
+
+
 #%% vocabulary types
 
 Array = NDArray[np.float64]
@@ -30,6 +94,119 @@ class COO(NamedTuple):
     j: ArrayI
     v: Array
 
+    def coo(self):
+
+        return sp.coo_matrix((self.v, (self.i, self.j)))
+
+    def csc(self):
+
+        return self.coo().tocsc()
+
+    def csr(self):
+
+        return self.coo().tocsr()
+
+
+class Curve(NamedTuple):
+    """
+    B-spline curve container.
+    """
+
+    pts: Array
+    knots: Array
+    order: int
+    periodic: bool
+
+    def curve(self) -> Geom_BSplineCurve:
+
+        if self.periodic:
+            mults = _colIntArray(np.ones_like(self.knots, dtype=int))
+            knots = _colRealArray(self.knots)
+        else:
+            unique_knots, mults_arr = np.unique(self.knots, return_counts=True)
+            knots = _colRealArray(unique_knots)
+            mults = _colIntArray(mults_arr)
+
+        return Geom_BSplineCurve(
+            _colPtsArray(self.pts), knots, mults, self.order, self.periodic,
+        )
+
+    def edge(self) -> Edge:
+
+        return Edge(BRepBuilderAPI_MakeEdge(self.curve()).Shape())
+
+    @classmethod
+    def fromEdge(cls, e: Edge):
+
+        assert (
+            e.geomType() == "BSPLINE"
+        ), "B-spline geometry required, try converting first."
+
+        g = e._geomAdaptor().BSpline()
+
+        knots = np.array(list(e._geomAdaptor().BSpline().KnotSequence()))
+        pts = np.array([(p.X(), p.Y(), p.Z()) for p in g.Poles()])
+        order = g.Degree()
+        periodic = g.IsPeriodic()
+
+        return cls(pts, knots, order, periodic)
+
+    def __call__(self, us: NDArray) -> NDArray:
+
+        pass
+
+    def der(self, us: NDArray) -> NDArray:
+
+        pass
+
+
+class Surface(NamedTuple):
+    """
+    B-spline surface container.
+    """
+
+    pts: Array
+    uknots: Array
+    vknots: Array
+    uorder: int
+    vorder: int
+    uperiodic: bool
+    vperiodic: bool
+
+    def surface(self) -> Geom_BSplineSurface:
+
+        if self.uperiodic:
+            umults = _colIntArray(np.ones_like(self.uknots, dtype=int))
+            uknots = _colRealArray(self.uknots)
+        else:
+            unique_knots, mults_arr = np.unique(self.uknots, return_counts=True)
+            uknots = _colRealArray(unique_knots)
+            umults = _colIntArray(mults_arr)
+
+        if self.vperiodic:
+            vmults = _colIntArray(np.ones_like(self.vknots, dtype=int))
+            vknots = _colRealArray(self.vknots)
+        else:
+            unique_knots, mults_arr = np.unique(self.vknots, return_counts=True)
+            vknots = _colRealArray(unique_knots)
+            vmults = _colIntArray(mults_arr)
+
+        return Geom_BSplineSurface(
+            _colPtsArray2(self.pts),
+            uknots,
+            vknots,
+            umults,
+            vmults,
+            self.uorder,
+            self.vorder,
+            self.uperiodic,
+            self.vperiodic,
+        )
+
+    def face(self, tol: float = 1e-3) -> Face:
+
+        return Face(BRepBuilderAPI_MakeFace(self.surface(), tol).Shape())
+
 
 #%% basis functions
 
@@ -423,6 +600,304 @@ def periodicDerMatrix(u: Array, order: int, dorder: int, knots: Array) -> list[C
     return rv
 
 
+@njit
+def periodicDiscretePenalty(us: Array, order: int) -> COO:
+
+    if order not in (1, 2):
+        raise ValueError(
+            f"Only 1st and 2nd order penalty is supported, requested order {order}"
+        )
+
+    # number of rows
+    nb = len(us)
+
+    # number of elements per row
+    ne = order + 1
+
+    # initialize the penlaty matrix
+    rv = COO(
+        i=np.empty(nb * ne, dtype=np.int64),
+        j=np.empty(nb * ne, dtype=np.int64),
+        v=np.empty(nb * ne),
+    )
+
+    if order == 1:
+        for ix in range(nb):
+            rv.i[ne * ix] = ix
+            rv.j[ne * ix] = (ix - 1) % nb
+            rv.v[ne * ix] = -0.5
+
+            rv.i[ne * ix + 1] = ix
+            rv.j[ne * ix + 1] = (ix + 1) % nb
+            rv.v[ne * ix + 1] = 0.5
+
+    elif order == 2:
+        for ix in range(nb):
+            rv.i[ne * ix] = ix
+            rv.j[ne * ix] = (ix - 1) % nb
+            rv.v[ne * ix] = 1
+
+            rv.i[ne * ix + 1] = ix
+            rv.j[ne * ix + 1] = ix
+            rv.v[ne * ix + 1] = -2
+
+            rv.i[ne * ix + 2] = ix
+            rv.j[ne * ix + 2] = (ix + 1) % nb
+            rv.v[ne * ix + 2] = 1
+
+    return rv
+
+
+@njit
+def discretePenalty(us: Array, order: int, splineorder: int = 3) -> COO:
+
+    if order not in (1, 2):
+        raise ValueError(
+            f"Only 1st and 2nd order penalty is supported, requested order {order}"
+        )
+
+    # number of rows
+    nb = len(us)
+
+    # number of elements per row
+    ne = order + 1
+
+    # initialize the penlaty matrix
+    rv = COO(
+        i=np.empty(nb * ne, dtype=np.int64),
+        j=np.empty(nb * ne, dtype=np.int64),
+        v=np.empty(nb * ne),
+    )
+
+    if order == 1:
+        for ix in range(nb):
+            if ix == 0:
+                rv.i[ne * ix] = ix
+                rv.j[ne * ix] = ix
+                rv.v[ne * ix] = -1
+
+                rv.i[ne * ix + 1] = ix
+                rv.j[ne * ix + 1] = ix + 1
+                rv.v[ne * ix + 1] = 1
+            elif ix < nb - 1:
+                rv.i[ne * ix] = ix
+                rv.j[ne * ix] = ix - 1
+                rv.v[ne * ix] = -0.5
+
+                rv.i[ne * ix + 1] = ix
+                rv.j[ne * ix + 1] = ix + 1
+                rv.v[ne * ix + 1] = 0.5
+            else:
+                rv.i[ne * ix] = ix
+                rv.j[ne * ix] = ix - 1
+                rv.v[ne * ix] = -1
+
+                rv.i[ne * ix + 1] = ix
+                rv.j[ne * ix + 1] = ix
+                rv.v[ne * ix + 1] = 1
+
+    elif order == 2:
+        for ix in range(nb):
+            if ix == 0:
+                rv.i[ne * ix] = ix
+                rv.j[ne * ix] = ix
+                rv.v[ne * ix] = 1
+
+                rv.i[ne * ix + 1] = ix
+                rv.j[ne * ix + 1] = ix + 1
+                rv.v[ne * ix + 1] = -2
+
+                rv.i[ne * ix + 2] = ix
+                rv.j[ne * ix + 2] = ix + 2
+                rv.v[ne * ix + 2] = 1
+            elif ix < nb - 1:
+                rv.i[ne * ix] = ix
+                rv.j[ne * ix] = ix - 1
+                rv.v[ne * ix] = 1
+
+                rv.i[ne * ix + 1] = ix
+                rv.j[ne * ix + 1] = ix
+                rv.v[ne * ix + 1] = -2
+
+                rv.i[ne * ix + 2] = ix
+                rv.j[ne * ix + 2] = ix + 1
+                rv.v[ne * ix + 2] = 1
+            else:
+                rv.i[ne * ix] = ix
+                rv.j[ne * ix] = ix - 2
+                rv.v[ne * ix] = 1
+
+                rv.i[ne * ix + 1] = ix
+                rv.j[ne * ix + 1] = ix - 1
+                rv.v[ne * ix + 1] = -2
+
+                rv.i[ne * ix + 2] = ix
+                rv.j[ne * ix + 2] = ix
+                rv.v[ne * ix + 2] = 1
+
+    return rv
+
+
+def periodicApproximate(
+    data: NDArray,
+    us: Optional[NDArray] = None,
+    knots: int | NDArray = 50,
+    order: int = 3,
+) -> Curve:
+
+    npts = data.shape[0]
+
+    # parametrize the points
+    us = linspace(0, 1, npts, endpoint=False)
+
+    # construct the knot vector
+    if isinstance(knots, int):
+        knots_ = linspace(0, 1, knots)
+    else:
+        knots_ = np.array(knots)
+
+    # construct the design matrix
+    C = periodicDesignMatrix(us, order, knots_).csc()
+    CtC = C.T @ C
+
+    # factorize
+    D, L, P = ldl(CtC, True)
+
+    # invert
+    pts = ldl_solve(C.T @ data, D, L, P).toarray()
+
+    # convert to an edge
+    rv = Curve(pts, knots_, order, periodic=True)
+
+    return rv
+
+
+def approximate(
+    data: NDArray,
+    us: Optional[NDArray] = None,
+    knots: int | NDArray = 50,
+    order: int = 3,
+    penalty: int = 4,
+    lam: float = 0,
+) -> Curve:
+
+    npts = data.shape[0]
+
+    # parametrize the points
+    us = linspace(0, 1, npts)
+
+    # construct the knot vector
+    if isinstance(knots, int):
+        knots_ = np.concatenate(
+            (np.repeat(0, order), linspace(0, 1, knots), np.repeat(1, order))
+        )
+    else:
+        knots_ = np.array(knots)
+
+    # construct the design matrix
+    C = designMatrix(us, order, knots_).csc()
+    CtC = C.T @ C
+
+    # add a penalty term if requested
+    if lam:
+        up = linspace(0, 1, order * npts)
+
+        assert penalty <= order + 2
+
+        # discrete + exact derivatives
+        if penalty > order:
+            Pexact = derMatrix(up, order, order - 1, knots_)[-1].csc()
+            Pdiscrete = discretePenalty(up, penalty - order, order).csc()
+
+            P = Pdiscrete @ Pexact
+
+        # only exact derivatives
+        else:
+            P = derMatrix(up, order, penalty, knots_)[-1].csc()
+
+        CtC += lam * P.T @ P
+
+    # clamp first and last point
+    Cc = C[[0, -1], :]
+    bc = data[[0, -1], :]
+
+    # final matrix and vector
+    Aug = sp.bmat([[CtC, Cc.T], [Cc, None]])
+    data_aug = np.vstack((C.T @ data, bc))
+
+    # factorize
+    D, L, P = ldl(Aug, False)
+
+    # invert
+    pts = ldl_solve(data_aug, D, L, P).toarray()[:-2, :]
+
+    # convert to an edge
+    rv = Curve(pts, knots_, order, periodic=False)
+
+    return rv
+
+
+def periodicLoft(*curves: Curve, order: int = 3) -> Curve:
+
+    nknots: int = len(curves) + 1
+
+    # collect control pts
+    pts = np.stack([c.pts for c in curves])
+
+    # approximate
+    pts_new = []
+
+    for j in range(pts.shape[1]):
+        pts_new.append(periodicApproximate(pts[:, j, :], knots=nknots, order=order).pts)
+
+    # construct the final surface
+    rv = Surface(
+        np.stack(pts_new).swapaxes(0, 1),
+        linspace(0, 1, nknots),
+        curves[0].knots,
+        order,
+        curves[0].order,
+        True,
+        curves[0].periodic,
+    )
+
+    return rv
+
+
+def loft(*curves: Curve, order: int = 3, lam: float = 1e-9, penalty: int = 4):
+
+    nknots: int = len(curves)
+
+    # collect control pts
+    pts = np.stack([c.pts for c in curves])
+
+    # approximate
+    pts_new = []
+
+    for j in range(pts.shape[1]):
+        pts_new.append(
+            approximate(
+                pts[:, j, :], knots=nknots, order=order, lam=lam, penalty=penalty
+            ).pts
+        )
+
+    # construct the final surface
+    rv = Surface(
+        np.stack(pts_new).swapaxes(0, 1),
+        np.concatenate(
+            (np.repeat(0, order), linspace(0, 1, nknots), np.repeat(1, order))
+        ),
+        curves[0].knots,
+        order,
+        curves[0].order,
+        False,
+        curves[0].periodic,
+    )
+
+    return rv
+
+
+#%% for removal?
 @njit
 def findSpan(v, knots):
 
diff --git a/tests/test_nurbs.py b/tests/test_nurbs.py
index 490d42d77..20308831c 100644
--- a/tests/test_nurbs.py
+++ b/tests/test_nurbs.py
@@ -4,11 +4,50 @@
     nbFindSpan,
     nbBasis,
     nbBasisDer,
+    Curve,
+    Surface,
+    approximate,
+    periodicApproximate,
+    periodicLoft,
+    loft,
 )
 
+from cadquery.func import circle
+
 import numpy as np
 import scipy.sparse as sp
 
+from pytest import approx, fixture
+
+
+@fixture
+def circles() -> list[Curve]:
+
+    # u,v periodic
+    c1 = circle(1).toSplines()
+    c2 = circle(5)
+
+    cs = [
+        Curve.fromEdge(c1.moved(loc))
+        for loc in c2.locations(np.linspace(0, 1, 10, False))
+    ]
+
+    return cs
+
+
+@fixture
+def trimmed_circles() -> list[Curve]:
+
+    c1 = circle(1).trim(0, 1).toSplines()
+    c2 = circle(5)
+
+    cs = [
+        Curve.fromEdge(c1.moved(loc))
+        for loc in c2.locations(np.linspace(0, 1, 10, False))
+    ]
+
+    return cs
+
 
 def test_periodic_dm():
 
@@ -53,3 +92,122 @@ def test_der():
 
         # sanity check
         assert np.allclose(out_der[0, :], out)
+
+
+def test_periodic_curve():
+
+    knots = np.linspace(0, 1, 5)
+    pts = np.array([[0, 0, 0], [0, 0, 1], [0, 1, 2], [0, 2, 0]])
+
+    crv = Curve(pts, knots, 3, True)
+
+    # is it indeed periodic?
+    assert crv.curve().IsPeriodic()
+
+    # convert to an edge
+    e = crv.edge()
+
+    assert e.isValid()
+    assert e.ShapeType() == "Edge"
+
+
+def test_curve():
+
+    knots = np.array([0, 0, 0, 0, 1, 1, 1, 1])
+    pts = np.array([[0, 0, 0], [0, 0, 1], [0, 1, 2], [0, 2, 0]])
+
+    crv = Curve(pts, knots, 3, False)
+
+    # sanity check
+    assert not crv.curve().IsPeriodic()
+
+    # convert to an edge
+    e = crv.edge()
+
+    assert e.isValid()
+    assert e.ShapeType() == "Edge"
+
+    # edge to curve
+    crv2 = Curve.fromEdge(e)
+    e2 = crv2.edge()
+
+    assert e2.isValid()
+
+
+def test_surface():
+
+    uknots = vknots = np.array([0, 0, 1, 1])
+    pts = np.array([[[0, 0, 0], [0, 1, 0]], [[1, 0, 0], [1, 1, 0]]])
+
+    srf = Surface(pts, uknots, vknots, 1, 1, False, False)
+
+    # convert to a face
+    f = srf.face()
+
+    assert f.isValid()
+    assert f.Area() == approx(1)
+
+
+def test_approximate():
+
+    pts_ = circle(1).trim(0, 1).sample(100)[0]
+    pts = np.array([list(p) for p in pts_])
+
+    # regular approximate
+    crv = approximate(pts)
+    e = crv.edge()
+
+    assert e.isValid()
+    assert e.Length() == approx(1)
+
+    # approximate with a  double penalty
+    crv = approximate(pts, penalty=4, lam=1e-9)
+    e = crv.edge()
+
+    assert e.isValid()
+    assert e.Length() == approx(1)
+
+    # approximate with a single penalty
+    crv = approximate(pts, penalty=2, lam=1e-9)
+    e = crv.edge()
+
+    assert e.isValid()
+    assert e.Length() == approx(1)
+
+
+def test_periodic_approximate():
+
+    pts_ = circle(1).sample(100)[0]
+    pts = np.array([list(p) for p in pts_])
+
+    crv = periodicApproximate(pts)
+    e = crv.edge()
+
+    assert e.isValid()
+    assert e.Length() == approx(2 * np.pi)
+
+
+def test_periodic_loft(circles, trimmed_circles):
+
+    # u,v periodic
+    surf1 = periodicLoft(*circles)
+
+    assert surf1.face().isValid()
+
+    # u periodic
+    surf2 = periodicLoft(*trimmed_circles)
+
+    assert surf2.face().isValid()
+
+
+def test_loft(circles, trimmed_circles):
+
+    # v periodic
+    surf1 = loft(*circles)
+
+    assert surf1.face().isValid()
+
+    # non-periodic
+    surf2 = loft(*trimmed_circles)
+
+    assert surf2.face().isValid()

From 4ffe9ff78e589dfc3fc6663442de01d0f9e04420 Mon Sep 17 00:00:00 2001
From: adam-urbanczyk <13981538+adam-urbanczyk@users.noreply.github.com>
Date: Fri, 4 Jul 2025 18:23:09 +0200
Subject: [PATCH 30/55] Mypy fixes

---
 cadquery/occ_impl/nurbs.py | 33 +++++++++++++++++++++++----------
 1 file changed, 23 insertions(+), 10 deletions(-)

diff --git a/cadquery/occ_impl/nurbs.py b/cadquery/occ_impl/nurbs.py
index d9f07a92f..5f37e2b49 100644
--- a/cadquery/occ_impl/nurbs.py
+++ b/cadquery/occ_impl/nurbs.py
@@ -5,7 +5,7 @@
 
 from numba import njit as _njit, prange
 
-from typing import NamedTuple, Optional
+from typing import NamedTuple, Optional, Tuple
 
 from numpy.typing import NDArray
 from numpy import linspace, array, empty_like, atleast_1d
@@ -34,7 +34,7 @@
 #%% internal helpers
 
 
-def _colPtsArray(pts: NDArray):
+def _colPtsArray(pts: NDArray) -> TColgp_Array1OfPnt:
 
     rv = TColgp_Array1OfPnt(1, pts.shape[0])
 
@@ -44,7 +44,7 @@ def _colPtsArray(pts: NDArray):
     return rv
 
 
-def _colPtsArray2(pts: NDArray) -> TColStd_Array2OfReal:
+def _colPtsArray2(pts: NDArray) -> TColgp_Array2OfPnt:
 
     assert pts.ndim == 3
 
@@ -59,7 +59,7 @@ def _colPtsArray2(pts: NDArray) -> TColStd_Array2OfReal:
     return rv
 
 
-def _colRealArray(knots: NDArray):
+def _colRealArray(knots: NDArray) -> TColStd_Array1OfReal:
 
     rv = TColStd_Array1OfReal(1, len(knots))
 
@@ -69,7 +69,7 @@ def _colRealArray(knots: NDArray):
     return rv
 
 
-def _colIntArray(knots: NDArray):
+def _colIntArray(knots: NDArray) -> TColStd_Array1OfInteger:
 
     rv = TColStd_Array1OfInteger(1, len(knots))
 
@@ -153,11 +153,11 @@ def fromEdge(cls, e: Edge):
 
     def __call__(self, us: NDArray) -> NDArray:
 
-        pass
+        raise NotImplementedError()
 
     def der(self, us: NDArray) -> NDArray:
 
-        pass
+        raise NotImplementedError()
 
 
 class Surface(NamedTuple):
@@ -779,6 +779,7 @@ def approximate(
     order: int = 3,
     penalty: int = 4,
     lam: float = 0,
+    tangents: Optional[Tuple[NDArray, NDArray]] = None,
 ) -> Curve:
 
     npts = data.shape[0]
@@ -820,6 +821,16 @@ def approximate(
     # clamp first and last point
     Cc = C[[0, -1], :]
     bc = data[[0, -1], :]
+    nc = 2  # number of constraints
+
+    # handle tangent constraints if needed
+    if tangents:
+        nc += 2
+
+        Cc2 = derMatrix(us[[0, -1]], order, 1, knots_)[-1].csc()
+
+        Cc = sp.vstack((Cc, Cc2))
+        bc = np.vstack((bc, *tangents))
 
     # final matrix and vector
     Aug = sp.bmat([[CtC, Cc.T], [Cc, None]])
@@ -829,7 +840,7 @@ def approximate(
     D, L, P = ldl(Aug, False)
 
     # invert
-    pts = ldl_solve(data_aug, D, L, P).toarray()[:-2, :]
+    pts = ldl_solve(data_aug, D, L, P).toarray()[:-nc, :]
 
     # convert to an edge
     rv = Curve(pts, knots_, order, periodic=False)
@@ -837,7 +848,7 @@ def approximate(
     return rv
 
 
-def periodicLoft(*curves: Curve, order: int = 3) -> Curve:
+def periodicLoft(*curves: Curve, order: int = 3) -> Surface:
 
     nknots: int = len(curves) + 1
 
@@ -864,7 +875,9 @@ def periodicLoft(*curves: Curve, order: int = 3) -> Curve:
     return rv
 
 
-def loft(*curves: Curve, order: int = 3, lam: float = 1e-9, penalty: int = 4):
+def loft(
+    *curves: Curve, order: int = 3, lam: float = 1e-9, penalty: int = 4
+) -> Surface:
 
     nknots: int = len(curves)
 

From 0ec8ee5d6618cd91f437af264b84179742065710 Mon Sep 17 00:00:00 2001
From: adam-urbanczyk <13981538+adam-urbanczyk@users.noreply.github.com>
Date: Fri, 4 Jul 2025 18:33:30 +0200
Subject: [PATCH 31/55] Add tangents to loft

---
 cadquery/occ_impl/nurbs.py | 15 ++++++++++++---
 1 file changed, 12 insertions(+), 3 deletions(-)

diff --git a/cadquery/occ_impl/nurbs.py b/cadquery/occ_impl/nurbs.py
index 5f37e2b49..434fb1371 100644
--- a/cadquery/occ_impl/nurbs.py
+++ b/cadquery/occ_impl/nurbs.py
@@ -5,7 +5,7 @@
 
 from numba import njit as _njit, prange
 
-from typing import NamedTuple, Optional, Tuple
+from typing import NamedTuple, Optional, Tuple, List
 
 from numpy.typing import NDArray
 from numpy import linspace, array, empty_like, atleast_1d
@@ -876,7 +876,11 @@ def periodicLoft(*curves: Curve, order: int = 3) -> Surface:
 
 
 def loft(
-    *curves: Curve, order: int = 3, lam: float = 1e-9, penalty: int = 4
+    *curves: Curve,
+    order: int = 3,
+    lam: float = 1e-9,
+    penalty: int = 4,
+    tangents: Optional[List[Tuple[NDArray, NDArray]]] = None,
 ) -> Surface:
 
     nknots: int = len(curves)
@@ -890,7 +894,12 @@ def loft(
     for j in range(pts.shape[1]):
         pts_new.append(
             approximate(
-                pts[:, j, :], knots=nknots, order=order, lam=lam, penalty=penalty
+                pts[:, j, :],
+                knots=nknots,
+                order=order,
+                lam=lam,
+                penalty=penalty,
+                tangents=tangents[j] if tangents else None,
             ).pts
         )
 

From 0a81707e425dcd19437808b6695350df52ce7bbc Mon Sep 17 00:00:00 2001
From: adam-urbanczyk <13981538+adam-urbanczyk@users.noreply.github.com>
Date: Sat, 5 Jul 2025 18:39:33 +0200
Subject: [PATCH 32/55] Mypy fix

---
 cadquery/occ_impl/nurbs.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/cadquery/occ_impl/nurbs.py b/cadquery/occ_impl/nurbs.py
index 434fb1371..e2ef464e2 100644
--- a/cadquery/occ_impl/nurbs.py
+++ b/cadquery/occ_impl/nurbs.py
@@ -81,7 +81,7 @@ def _colIntArray(knots: NDArray) -> TColStd_Array1OfInteger:
 
 #%% vocabulary types
 
-Array = NDArray[np.float64]
+Array = NDArray[np.floating]
 ArrayI = NDArray[np.int_]
 
 

From 1b7f1f71e5bf4c4184d557d62df682c5e76f7778 Mon Sep 17 00:00:00 2001
From: adam-urbanczyk <13981538+adam-urbanczyk@users.noreply.github.com>
Date: Thu, 10 Jul 2025 15:21:02 +0200
Subject: [PATCH 33/55] Add overloads

---
 cadquery/occ_impl/nurbs.py | 200 ++++++++++++++++++++++++++++++++++---
 1 file changed, 185 insertions(+), 15 deletions(-)

diff --git a/cadquery/occ_impl/nurbs.py b/cadquery/occ_impl/nurbs.py
index e2ef464e2..f9de5f506 100644
--- a/cadquery/occ_impl/nurbs.py
+++ b/cadquery/occ_impl/nurbs.py
@@ -1,14 +1,13 @@
 #%% imports
 import numpy as np
 import scipy.sparse as sp
-import math
 
-from numba import njit as _njit, prange
+from numba import njit as _njit
 
-from typing import NamedTuple, Optional, Tuple, List
+from typing import NamedTuple, Optional, Tuple, List, Union, cast, Type
 
 from numpy.typing import NDArray
-from numpy import linspace, array, empty_like, atleast_1d
+from numpy import linspace, ndarray, dtype
 
 from casadi import ldl, ldl_solve
 
@@ -17,13 +16,14 @@
 from OCP.TColStd import (
     TColStd_Array1OfInteger,
     TColStd_Array1OfReal,
-    TColStd_Array2OfReal,
 )
 from OCP.gp import gp_Pnt
 from OCP.BRepBuilderAPI import BRepBuilderAPI_MakeEdge, BRepBuilderAPI_MakeFace
 
 from .shapes import Face, Edge
 
+from multimethod import multidispatch, parametric
+
 njit = _njit(cache=False, error_model="numpy", fastmath=True, parallel=False)
 
 njiti = _njit(
@@ -81,8 +81,8 @@ def _colIntArray(knots: NDArray) -> TColStd_Array1OfInteger:
 
 #%% vocabulary types
 
-Array = NDArray[np.floating]
-ArrayI = NDArray[np.int_]
+Array = ndarray  # NDArray[np.floating]
+ArrayI = ndarray  # NDArray[np.int_]
 
 
 class COO(NamedTuple):
@@ -738,11 +738,14 @@ def discretePenalty(us: Array, order: int, splineorder: int = 3) -> COO:
     return rv
 
 
+@multidispatch
 def periodicApproximate(
-    data: NDArray,
-    us: Optional[NDArray] = None,
-    knots: int | NDArray = 50,
+    data: Array,
+    us: Optional[Array] = None,
+    knots: int | Array = 50,
     order: int = 3,
+    penalty: int = 4,
+    lam: float = 0,
 ) -> Curve:
 
     npts = data.shape[0]
@@ -760,6 +763,25 @@ def periodicApproximate(
     C = periodicDesignMatrix(us, order, knots_).csc()
     CtC = C.T @ C
 
+    # add the penalty if requested
+    if lam:
+        up = linspace(0, 1, order * npts, endpoint=False)
+
+        assert penalty <= order + 2
+
+        # discrete + exact derivatives
+        if penalty > order:
+            Pexact = periodicDerMatrix(up, order, order - 1, knots_)[-1].csc()
+            Pdiscrete = periodicDiscretePenalty(up, penalty - order).csc()
+
+            P = Pdiscrete @ Pexact
+
+        # only exact derivatives
+        else:
+            P = periodicDerMatrix(up, order, penalty, knots_)[-1].csc()
+
+        CtC += lam * P.T @ P
+
     # factorize
     D, L, P = ldl(CtC, True)
 
@@ -772,14 +794,74 @@ def periodicApproximate(
     return rv
 
 
+@periodicApproximate.register
+def _(
+    data: List[Array],
+    us: Optional[Array] = None,
+    knots: int | Array = 50,
+    order: int = 3,
+    penalty: int = 4,
+    lam: float = 0,
+) -> List[Curve]:
+
+    rv = []
+
+    npts = data[0].shape[0]
+
+    # parametrize the points
+    us = linspace(0, 1, npts, endpoint=False)
+
+    # construct the knot vector
+    if isinstance(knots, int):
+        knots_ = linspace(0, 1, knots)
+    else:
+        knots_ = np.array(knots)
+
+    # construct the design matrix
+    C = periodicDesignMatrix(us, order, knots_).csc()
+    CtC = C.T @ C
+
+    # add the penalty if requested
+    if lam:
+        up = linspace(0, 1, order * npts, endpoint=False)
+
+        assert penalty <= order + 2
+
+        # discrete + exact derivatives
+        if penalty > order:
+            Pexact = periodicDerMatrix(up, order, order - 1, knots_)[-1].csc()
+            Pdiscrete = periodicDiscretePenalty(up, penalty - order).csc()
+
+            P = Pdiscrete @ Pexact
+
+        # only exact derivatives
+        else:
+            P = periodicDerMatrix(up, order, penalty, knots_)[-1].csc()
+
+        CtC += lam * P.T @ P
+
+    # factorize
+    D, L, P = ldl(CtC, True)
+
+    # invert every dataset
+    for dataset in data:
+        pts = ldl_solve(C.T @ dataset, D, L, P).toarray()
+
+        # convert to an edge and store
+        rv.append(Curve(pts, knots_, order, periodic=True))
+
+    return rv
+
+
+@multidispatch
 def approximate(
-    data: NDArray,
-    us: Optional[NDArray] = None,
-    knots: int | NDArray = 50,
+    data: Array,
+    us: Optional[Array] = None,
+    knots: int | Array = 50,
     order: int = 3,
     penalty: int = 4,
     lam: float = 0,
-    tangents: Optional[Tuple[NDArray, NDArray]] = None,
+    tangents: Optional[Tuple[Array, Array]] = None,
 ) -> Curve:
 
     npts = data.shape[0]
@@ -848,6 +930,94 @@ def approximate(
     return rv
 
 
+@approximate.register
+def _(
+    data: List[Array],
+    us: Optional[Array] = None,
+    knots: int | Array = 50,
+    order: int = 3,
+    penalty: int = 4,
+    lam: float = 0,
+    tangents: Optional[Union[Tuple[Array, Array], List[Tuple[Array, Array]]]] = None,
+) -> List[Curve]:
+
+    rv = []
+
+    npts = data[0].shape[0]
+
+    # parametrize the points
+    us = linspace(0, 1, npts)
+
+    # construct the knot vector
+    if isinstance(knots, int):
+        knots_ = np.concatenate(
+            (np.repeat(0, order), linspace(0, 1, knots), np.repeat(1, order))
+        )
+    else:
+        knots_ = np.array(knots)
+
+    # construct the design matrix
+    C = designMatrix(us, order, knots_).csc()
+    CtC = C.T @ C
+
+    # add a penalty term if requested
+    if lam:
+        up = linspace(0, 1, order * npts)
+
+        assert penalty <= order + 2
+
+        # discrete + exact derivatives
+        if penalty > order:
+            Pexact = derMatrix(up, order, order - 1, knots_)[-1].csc()
+            Pdiscrete = discretePenalty(up, penalty - order, order).csc()
+
+            P = Pdiscrete @ Pexact
+
+        # only exact derivatives
+        else:
+            P = derMatrix(up, order, penalty, knots_)[-1].csc()
+
+        CtC += lam * P.T @ P
+
+    # clamp first and last point
+    Cc = C[[0, -1], :]
+
+    nc = 2  # number of constraints
+
+    # handle tangent constraints if needed
+    if tangents:
+        nc += 2
+        Cc2 = derMatrix(us[[0, -1]], order, 1, knots_)[-1].csc()
+        Cc = sp.vstack((Cc, Cc2))
+
+    # final matrix and vector
+    Aug = sp.bmat([[CtC, Cc.T], [Cc, None]])
+
+    # factorize
+    D, L, P = ldl(Aug, False)
+
+    # invert all datasets
+    for ix, dataset in enumerate(data):
+        bc = dataset[[0, -1], :]  # first and last point for clamping
+
+        if tangents:
+            if len(tangents) == len(data):
+                bc = np.vstack((bc, *tangents[ix]))
+            else:
+                bc = np.vstack((bc, *tangents))
+
+        # construct the LHS of the linear system
+        dataset_aug = np.vstack((C.T @ dataset, bc))
+
+        # actual solver
+        pts = ldl_solve(dataset_aug, D, L, P).toarray()[:-nc, :]
+
+        # convert to an edge
+        rv.append(Curve(pts, knots_, order, periodic=False))
+
+    return rv
+
+
 def periodicLoft(*curves: Curve, order: int = 3) -> Surface:
 
     nknots: int = len(curves) + 1
@@ -880,7 +1050,7 @@ def loft(
     order: int = 3,
     lam: float = 1e-9,
     penalty: int = 4,
-    tangents: Optional[List[Tuple[NDArray, NDArray]]] = None,
+    tangents: Optional[List[Tuple[Array, Array]]] = None,
 ) -> Surface:
 
     nknots: int = len(curves)

From 17e18c57a504c8eff042843e2d09a00ef4b933c8 Mon Sep 17 00:00:00 2001
From: adam-urbanczyk <13981538+adam-urbanczyk@users.noreply.github.com>
Date: Thu, 10 Jul 2025 16:24:50 +0200
Subject: [PATCH 34/55] Faster periodic loft

---
 cadquery/occ_impl/nurbs.py | 7 ++-----
 1 file changed, 2 insertions(+), 5 deletions(-)

diff --git a/cadquery/occ_impl/nurbs.py b/cadquery/occ_impl/nurbs.py
index f9de5f506..e8c7fb47b 100644
--- a/cadquery/occ_impl/nurbs.py
+++ b/cadquery/occ_impl/nurbs.py
@@ -1023,13 +1023,10 @@ def periodicLoft(*curves: Curve, order: int = 3) -> Surface:
     nknots: int = len(curves) + 1
 
     # collect control pts
-    pts = np.stack([c.pts for c in curves])
+    pts = [el for el in np.stack([c.pts for c in curves]).swapaxes(0, 1)]
 
     # approximate
-    pts_new = []
-
-    for j in range(pts.shape[1]):
-        pts_new.append(periodicApproximate(pts[:, j, :], knots=nknots, order=order).pts)
+    pts_new = [el.pts for el in periodicApproximate(pts, knots=nknots, order=order)]
 
     # construct the final surface
     rv = Surface(

From 6b82f65a5dc9bad2016b14b2430ab31fb73a9b51 Mon Sep 17 00:00:00 2001
From: adam-urbanczyk <13981538+adam-urbanczyk@users.noreply.github.com>
Date: Tue, 15 Jul 2025 12:46:45 +0200
Subject: [PATCH 35/55] Mypy fix

---
 cadquery/occ_impl/nurbs.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/cadquery/occ_impl/nurbs.py b/cadquery/occ_impl/nurbs.py
index e8c7fb47b..9fd22334b 100644
--- a/cadquery/occ_impl/nurbs.py
+++ b/cadquery/occ_impl/nurbs.py
@@ -795,7 +795,7 @@ def periodicApproximate(
 
 
 @periodicApproximate.register
-def _(
+def periodicApproximate(
     data: List[Array],
     us: Optional[Array] = None,
     knots: int | Array = 50,
@@ -931,7 +931,7 @@ def approximate(
 
 
 @approximate.register
-def _(
+def approximate(
     data: List[Array],
     us: Optional[Array] = None,
     knots: int | Array = 50,

From 553078762112c6a6650670983549d4323a9d5f24 Mon Sep 17 00:00:00 2001
From: adam-urbanczyk <13981538+adam-urbanczyk@users.noreply.github.com>
Date: Tue, 15 Jul 2025 16:12:29 +0200
Subject: [PATCH 36/55] Start with evaluation

---
 cadquery/occ_impl/nurbs.py | 199 ++++++++++++++++++++++++++++++++++++-
 1 file changed, 194 insertions(+), 5 deletions(-)

diff --git a/cadquery/occ_impl/nurbs.py b/cadquery/occ_impl/nurbs.py
index 9fd22334b..25f04cd77 100644
--- a/cadquery/occ_impl/nurbs.py
+++ b/cadquery/occ_impl/nurbs.py
@@ -211,6 +211,28 @@ def face(self, tol: float = 1e-3) -> Face:
 #%% basis functions
 
 
+@njiti
+def extendKnots(order: int, knots: Array) -> Array:
+    """
+    Knot vector extension for periodic b-splines.
+
+    Parameters
+    ----------
+    order : int
+        B-spline order.
+    knots : Array
+        Knot vector.
+
+    Returns
+    -------
+    knots_ext : Array
+        Extended knots vector.
+
+    """
+
+    return np.concat((knots[-order:-1] - knots[-1], knots, knots[-1] + knots[1:order]))
+
+
 @njiti
 def nbFindSpan(
     u: float,
@@ -229,7 +251,7 @@ def nbFindSpan(
     order : int
         Spline order.
     knots : ndarray
-        Knot vectr.
+        Knot vector.
 
     Returns
     -------
@@ -275,7 +297,7 @@ def nbBasis(i: int, u: float, order: int, knots: Array, out: Array):
     order : int
         B-spline order.
     knots : ndarray
-        Knot vectr.
+        Knot vector.
     out : ndarray
         B-spline basis function values.
 
@@ -320,7 +342,7 @@ def nbBasisDer(i: int, u: float, order: int, dorder: int, knots: Array, out: Arr
     dorder : int
         Derivative order.
     knots : ndarray
-        Knot vectr.
+        Knot vector.
     out : ndarray
         B-spline basis function and derivative values.
 
@@ -408,6 +430,170 @@ def nbBasisDer(i: int, u: float, order: int, dorder: int, knots: Array, out: Arr
         r *= order - k
 
 
+#%% evaluation
+
+
+@njit
+def nbCurve(
+    u: Array, order: int, knots: Array, pts: Array, periodic: bool = False
+) -> Array:
+    """
+    NURBS book A3.1 with modifications to handle periodicity.
+
+    Parameters
+    ----------
+    u : Array
+        Parameter values.
+    order : int
+        B-spline order.
+    knots : Array
+        Knot vector.
+    pts : Array
+        Control points.
+    periodic : bool, optional
+        Peridocity flag. The default is False.
+
+    Returns
+    -------
+    Array
+        Curve values.
+
+    """
+
+    # number of control points
+    nb = pts.shape[0]
+
+    # handle periodicity
+    if periodic:
+        period = knots[-1] - knots[0]
+        u_ = u % period
+        knots_ext = extendKnots(order, knots)
+        minspan = 0
+        maxspan = len(knots) - 1
+        deltaspan = order - 1
+    else:
+        u_ = u
+        knots_ext = knots
+        minspan = None
+        maxspan = None
+        deltaspan = 0
+
+    # number of param values
+    nu = np.size(u)
+
+    # chunck size
+    n = order + 1
+
+    # temp chunck storage
+    temp = np.zeros(n)
+
+    # initialize
+    out = np.zeros((nu, 3))
+
+    for i in range(nu):
+        ui = u_[i]
+
+        # find span
+        span = nbFindSpan(ui, order, knots, minspan, maxspan) + deltaspan
+
+        # evaluate chunk
+        nbBasis(span, ui, order, knots_ext, temp)
+
+        # multiply by ctrl points
+        for j in range(order + 1):
+            out[i, :] += temp[j] * pts[(span - order + j) % nb, :]
+
+    return out
+
+
+def nbCurveDer(
+    u: Array, order: int, dorder: int, knots: Array, pts: Array, periodic: bool = False
+) -> Array:
+    """
+    NURBS book A3.2 with modifications to handle periodicity.
+
+    Parameters
+    ----------
+    u : Array
+        Parameter values.
+    order : int
+        B-spline order.
+    dorder : int
+        Derivative order.
+    knots : Array
+        Knot vector.
+    pts : Array
+        Control points.
+    periodic : bool, optional
+        Peridocity flag. The default is False.
+
+
+    Returns
+    -------
+    Array
+        Curve values and derivatives.
+
+    """
+    # number of control points
+    nb = pts.shape[0]
+
+    # handle periodicity
+    if periodic:
+        period = knots[-1] - knots[0]
+        u_ = u % period
+        knots_ext = extendKnots(order, knots)
+        minspan = 0
+        maxspan = len(knots) - 1
+        deltaspan = order - 1
+    else:
+        u_ = u
+        knots_ext = knots
+        minspan = None
+        maxspan = None
+        deltaspan = 0
+
+    # number of param values
+    nu = np.size(u)
+
+    # chunck size
+    n = order + 1
+
+    # temp chunck storage
+    temp = np.zeros((dorder + 1, n))
+
+    # initialize
+    out = np.zeros((nu, dorder + 1, 3))
+
+    for i in range(nu):
+        ui = u_[i]
+
+        # find span
+        span = nbFindSpan(ui, order, knots, minspan, maxspan) + deltaspan
+
+        # evaluate chunk
+        nbBasisDer(span, ui, order, dorder, knots_ext, temp)
+
+        # multiply by ctrl points
+        for j in range(order + 1):
+            for k in range(dorder + 1):
+                out[i, k, :] += temp[k, j] * pts[(span - order + j) % nb, :]
+
+    return out
+
+
+def nbSurface():
+
+    pass
+
+
+def nbSurfaceDer():
+
+    pass
+
+
+#%% matrices
+
+
 @njit
 def designMatrix(u: Array, order: int, knots: Array) -> COO:
     """
@@ -511,7 +697,7 @@ def derMatrix(u: Array, order: int, dorder: int, knots: Array) -> list[COO]:
     n = order + 1
 
     # temp chunck storage
-    temp = np.zeros((n, n))
+    temp = np.zeros((dorder + 1, n))
 
     # initialize the empty matrix
     rv = []
@@ -565,7 +751,7 @@ def periodicDerMatrix(u: Array, order: int, dorder: int, knots: Array) -> list[C
     n = order + 1
 
     # temp chunck storage
-    temp = np.zeros((n, n))
+    temp = np.zeros((dorder + 1, n))
 
     # initialize the empty matrix
     rv = []
@@ -738,6 +924,9 @@ def discretePenalty(us: Array, order: int, splineorder: int = 3) -> COO:
     return rv
 
 
+#%% construction
+
+
 @multidispatch
 def periodicApproximate(
     data: Array,

From e96a10fe8473e16ffb6eade8e7d9ea0d30cd111a Mon Sep 17 00:00:00 2001
From: adam-urbanczyk <13981538+adam-urbanczyk@users.noreply.github.com>
Date: Tue, 15 Jul 2025 16:23:00 +0200
Subject: [PATCH 37/55] Add decorator

---
 cadquery/occ_impl/nurbs.py | 1 +
 1 file changed, 1 insertion(+)

diff --git a/cadquery/occ_impl/nurbs.py b/cadquery/occ_impl/nurbs.py
index 25f04cd77..756f78611 100644
--- a/cadquery/occ_impl/nurbs.py
+++ b/cadquery/occ_impl/nurbs.py
@@ -506,6 +506,7 @@ def nbCurve(
     return out
 
 
+@njit
 def nbCurveDer(
     u: Array, order: int, dorder: int, knots: Array, pts: Array, periodic: bool = False
 ) -> Array:

From 4f96f6be349b5b831fb914f9f6f5e9fd8b74af59 Mon Sep 17 00:00:00 2001
From: adam-urbanczyk <13981538+adam-urbanczyk@users.noreply.github.com>
Date: Wed, 16 Jul 2025 12:31:47 +0200
Subject: [PATCH 38/55] Add params to nbSurface

---
 cadquery/occ_impl/nurbs.py | 12 +++++++++++-
 1 file changed, 11 insertions(+), 1 deletion(-)

diff --git a/cadquery/occ_impl/nurbs.py b/cadquery/occ_impl/nurbs.py
index 756f78611..8a17bb7a4 100644
--- a/cadquery/occ_impl/nurbs.py
+++ b/cadquery/occ_impl/nurbs.py
@@ -582,7 +582,17 @@ def nbCurveDer(
     return out
 
 
-def nbSurface():
+def nbSurface(
+    u: Array,
+    v: Array,
+    uorder: int,
+    vorder: int,
+    uknots: Array,
+    vknots: Array,
+    pts: Array,
+    uperiodic: bool = False,
+    vperiodic: bool = False,
+):
 
     pass
 

From 9f98643e682c4b4199065e31585be795554f061d Mon Sep 17 00:00:00 2001
From: adam-urbanczyk <13981538+adam-urbanczyk@users.noreply.github.com>
Date: Wed, 16 Jul 2025 23:47:49 +0200
Subject: [PATCH 39/55] Start with surface evaluation

---
 cadquery/occ_impl/nurbs.py | 204 ++++++++++++++++++++++++++++++++++---
 1 file changed, 190 insertions(+), 14 deletions(-)

diff --git a/cadquery/occ_impl/nurbs.py b/cadquery/occ_impl/nurbs.py
index 8a17bb7a4..6c9c3f97f 100644
--- a/cadquery/occ_impl/nurbs.py
+++ b/cadquery/occ_impl/nurbs.py
@@ -4,10 +4,10 @@
 
 from numba import njit as _njit
 
-from typing import NamedTuple, Optional, Tuple, List, Union, cast, Type
+from typing import NamedTuple, Optional, Tuple, List, Union, cast
 
 from numpy.typing import NDArray
-from numpy import linspace, ndarray, dtype
+from numpy import linspace, ndarray
 
 from casadi import ldl, ldl_solve
 
@@ -19,12 +19,13 @@
 )
 from OCP.gp import gp_Pnt
 from OCP.BRepBuilderAPI import BRepBuilderAPI_MakeEdge, BRepBuilderAPI_MakeFace
+from OCP.Geom import Geom_BSplineSurface
 
 from .shapes import Face, Edge
 
-from multimethod import multidispatch, parametric
+from multimethod import multidispatch
 
-njit = _njit(cache=False, error_model="numpy", fastmath=True, parallel=False)
+njit = _njit(cache=True, error_model="numpy", fastmath=True, nogil=True, parallel=False)
 
 njiti = _njit(
     cache=True, inline="always", error_model="numpy", fastmath=True, parallel=False
@@ -151,13 +152,13 @@ def fromEdge(cls, e: Edge):
 
         return cls(pts, knots, order, periodic)
 
-    def __call__(self, us: NDArray) -> NDArray:
+    def __call__(self, us: Array) -> Array:
 
-        raise NotImplementedError()
+        return nbCurve(np.atleast_1d(us), self.order, self.knots, self.pts)
 
-    def der(self, us: NDArray) -> NDArray:
+    def der(self, us: NDArray, dorder: int) -> NDArray:
 
-        raise NotImplementedError()
+        return nbCurveDer(np.atleast_1d(us), self.order, self.knots, self.pts)
 
 
 class Surface(NamedTuple):
@@ -207,6 +208,64 @@ def face(self, tol: float = 1e-3) -> Face:
 
         return Face(BRepBuilderAPI_MakeFace(self.surface(), tol).Shape())
 
+    @classmethod
+    def fromFace(cls, f: Face):
+        """
+        Construct a surface from a face.
+        """
+
+        assert (
+            f.geomType() == "BSPLINE"
+        ), "B-spline geometry required, try converting first."
+
+        g = cast(Geom_BSplineSurface, f._geomAdaptor())
+
+        uknots = np.array(list(g.UKnotSequence()))
+        vknots = np.array(list(g.VKnotSequence()))
+
+        tmp = []
+        for i in range(1, g.NbUPoles() + 1):
+            tmp.append(
+                [
+                    [g.Pole(i, j).X(), g.Pole(i, j).Y(), g.Pole(i, j).Z(),]
+                    for j in range(1, g.NbVPoles() + 1)
+                ]
+            )
+
+        pts = np.array(tmp)
+
+        uorder = g.UDegree()
+        vorder = g.VDegree()
+
+        uperiodic = g.IsUPeriodic()
+        vperiodic = g.IsVPeriodic()
+
+        return cls(pts, uknots, vknots, uorder, vorder, uperiodic, vperiodic)
+
+    def __call__(self, u: Array, v: Array) -> Array:
+        """
+        Evaluate surface at (u,v) points.
+        """
+
+        return nbSurface(
+            np.atleast_1d(u),
+            np.atleast_1d(v),
+            self.uorder,
+            self.vorder,
+            self.uknots,
+            self.vknots,
+            self.pts,
+            self.uperiodic,
+            self.vperiodic,
+        )
+
+    def der(self, u: Array, v: Array, dorder: int) -> Array:
+        """
+        Evaluate surface and derivatives at (u,v) points.
+        """
+
+        raise NotImplementedError
+
 
 #%% basis functions
 
@@ -451,7 +510,7 @@ def nbCurve(
     pts : Array
         Control points.
     periodic : bool, optional
-        Peridocity flag. The default is False.
+        Periodicity flag. The default is False.
 
     Returns
     -------
@@ -526,7 +585,7 @@ def nbCurveDer(
     pts : Array
         Control points.
     periodic : bool, optional
-        Peridocity flag. The default is False.
+        Periodicity flag. The default is False.
 
 
     Returns
@@ -582,6 +641,7 @@ def nbCurveDer(
     return out
 
 
+@njit
 def nbSurface(
     u: Array,
     v: Array,
@@ -592,14 +652,130 @@ def nbSurface(
     pts: Array,
     uperiodic: bool = False,
     vperiodic: bool = False,
-):
+) -> Array:
+    """
+    NURBS book A3.5 with modifications to handle periodicity.
+
+    Parameters
+    ----------
+    u : Array
+        U parameter values.
+    v : Array
+        V parameter values.
+    uorder : int
+        B-spline u order.
+    vorder : int
+        B-spline v order.
+    uknots : Array
+        U knot vector..
+    vknots : Array
+        V knot vector..
+    pts : Array
+        Control points.
+    uperiodic : bool, optional
+        U periodicity flag. The default is False.
+    vperiodic : bool, optional
+        V periodicity flag. The default is False.
+
+    Returns
+    -------
+    Array
+        Surface values.
+
+    """
+
+    # number of control points
+    nub = pts.shape[0]
+    nvb = pts.shape[1]
+
+    # handle periodicity
+    if uperiodic:
+        uperiod = uknots[-1] - uknots[0]
+        u_ = u % uperiod
+        uknots_ext = extendKnots(uorder, uknots)
+        minspanu = 0
+        maxspanu = len(uknots) - 1
+        deltaspanu = uorder - 1
+    else:
+        u_ = u
+        uknots_ext = uknots
+        minspanu = None
+        maxspanu = None
+        deltaspanu = 0
+
+    if vperiodic:
+        vperiod = vknots[-1] - vknots[0]
+        v_ = v % vperiod
+        vknots_ext = extendKnots(vorder, vknots)
+        minspanv = 0
+        maxspanv = len(vknots) - 1
+        deltaspanv = vorder - 1
+    else:
+        v_ = v
+        vknots_ext = vknots
+        minspanv = None
+        maxspanv = None
+        deltaspanv = 0
+
+    # number of param values
+    nu = np.size(u)
+
+    # chunck sizes
+    un = uorder + 1
+    vn = vorder + 1
+
+    # temp chunck storage
+    utemp = np.zeros(un)
+    vtemp = np.zeros(vn)
+
+    # initialize
+    out = np.zeros((nu, 3))
+
+    for i in range(nu):
+        ui = u_[i]
+        vi = v_[i]
+
+        # find span
+        uspan = nbFindSpan(ui, uorder, uknots, minspanu, maxspanu) + deltaspanu
+        vspan = nbFindSpan(vi, vorder, vknots, minspanv, maxspanv) + deltaspanv
+
+        # evaluate chunk
+        nbBasis(uspan, ui, uorder, uknots_ext, utemp)
+        nbBasis(vspan, vi, vorder, vknots_ext, vtemp)
+
+        uind = uspan - uorder
+        temp = np.empty(3)
+
+        # multiply by ctrl points: Nu.T*P*Nv
+        for j in range(vorder + 1):
 
-    pass
+            temp[:] = 0.0
+            vind = vspan - vorder + j
 
+            # calculate Nu.T*P
+            for k in range(uorder + 1):
+                temp += utemp[k] * pts[(uind + k) % nub, vind % nvb, :]
 
-def nbSurfaceDer():
+            # multiple by Nv
+            out[i, :] += vtemp[j] * temp
+
+    return out
+
+
+def nbSurfaceDer(
+    u: Array,
+    v: Array,
+    uorder: int,
+    vorder: int,
+    dorder: int,
+    uknots: Array,
+    vknots: Array,
+    pts: Array,
+    uperiodic: bool = False,
+    vperiodic: bool = False,
+) -> Array:
 
-    pass
+    raise NotImplementedError
 
 
 #%% matrices

From a9ce0c70a521b586326b347b8af699c4027d57c9 Mon Sep 17 00:00:00 2001
From: adam-urbanczyk <13981538+adam-urbanczyk@users.noreply.github.com>
Date: Thu, 17 Jul 2025 13:01:13 +0200
Subject: [PATCH 40/55] Add surface derivatives

---
 cadquery/occ_impl/nurbs.py | 140 +++++++++++++++++++++++++++++++++++--
 1 file changed, 136 insertions(+), 4 deletions(-)

diff --git a/cadquery/occ_impl/nurbs.py b/cadquery/occ_impl/nurbs.py
index 6c9c3f97f..f2389e845 100644
--- a/cadquery/occ_impl/nurbs.py
+++ b/cadquery/occ_impl/nurbs.py
@@ -154,11 +154,15 @@ def fromEdge(cls, e: Edge):
 
     def __call__(self, us: Array) -> Array:
 
-        return nbCurve(np.atleast_1d(us), self.order, self.knots, self.pts)
+        return nbCurve(
+            np.atleast_1d(us), self.order, self.knots, self.pts, self.periodic
+        )
 
     def der(self, us: NDArray, dorder: int) -> NDArray:
 
-        return nbCurveDer(np.atleast_1d(us), self.order, self.knots, self.pts)
+        return nbCurveDer(
+            np.atleast_1d(us), self.order, self.knots, self.pts, self.periodic
+        )
 
 
 class Surface(NamedTuple):
@@ -264,7 +268,18 @@ def der(self, u: Array, v: Array, dorder: int) -> Array:
         Evaluate surface and derivatives at (u,v) points.
         """
 
-        raise NotImplementedError
+        return nbSurfaceDer(
+            np.atleast_1d(u),
+            np.atleast_1d(v),
+            self.uorder,
+            self.vorder,
+            dorder,
+            self.uknots,
+            self.vknots,
+            self.pts,
+            self.uperiodic,
+            self.vperiodic,
+        )
 
 
 #%% basis functions
@@ -762,6 +777,7 @@ def nbSurface(
     return out
 
 
+@njit
 def nbSurfaceDer(
     u: Array,
     v: Array,
@@ -774,8 +790,124 @@ def nbSurfaceDer(
     uperiodic: bool = False,
     vperiodic: bool = False,
 ) -> Array:
+    """
+    NURBS book A3.6 with modifications to handle periodicity.
+
+    Parameters
+    ----------
+    u : Array
+        U parameter values.
+    v : Array
+        V parameter values.
+    uorder : int
+        B-spline u order.
+    vorder : int
+        B-spline v order.
+    dorder : int
+        Maximum derivative order.
+    uknots : Array
+        U knot vector..
+    vknots : Array
+        V knot vector..
+    pts : Array
+        Control points.
+    uperiodic : bool, optional
+        U periodicity flag. The default is False.
+    vperiodic : bool, optional
+        V periodicity flag. The default is False.
+
+    Returns
+    -------
+    Array
+        Surface and derivative values.
+
+    """
+
+    # max derivative orders
+    du = min(dorder, uorder)
+    dv = min(dorder, vorder)
+
+    # number of control points
+    nub = pts.shape[0]
+    nvb = pts.shape[1]
+
+    # handle periodicity
+    if uperiodic:
+        uperiod = uknots[-1] - uknots[0]
+        u_ = u % uperiod
+        uknots_ext = extendKnots(uorder, uknots)
+        minspanu = 0
+        maxspanu = len(uknots) - 1
+        deltaspanu = uorder - 1
+    else:
+        u_ = u
+        uknots_ext = uknots
+        minspanu = None
+        maxspanu = None
+        deltaspanu = 0
+
+    if vperiodic:
+        vperiod = vknots[-1] - vknots[0]
+        v_ = v % vperiod
+        vknots_ext = extendKnots(vorder, vknots)
+        minspanv = 0
+        maxspanv = len(vknots) - 1
+        deltaspanv = vorder - 1
+    else:
+        v_ = v
+        vknots_ext = vknots
+        minspanv = None
+        maxspanv = None
+        deltaspanv = 0
+
+    # number of param values
+    nu = np.size(u)
+
+    # chunck sizes
+    un = uorder + 1
+    vn = vorder + 1
+
+    # temp chunck storage
+
+    utemp = np.zeros((du + 1, un))
+    vtemp = np.zeros((dv + 1, vn))
+
+    # initialize
+    out = np.zeros((nu, du + 1, dv + 1, 3))
+
+    for i in range(nu):
+        ui = u_[i]
+        vi = v_[i]
+
+        # find span
+        uspan = nbFindSpan(ui, uorder, uknots, minspanu, maxspanu) + deltaspanu
+        vspan = nbFindSpan(vi, vorder, vknots, minspanv, maxspanv) + deltaspanv
+
+        # evaluate chunk
+        nbBasisDer(uspan, ui, uorder, du, uknots_ext, utemp)
+        nbBasisDer(vspan, vi, vorder, dv, vknots_ext, vtemp)
 
-    raise NotImplementedError
+        for k in range(du + 1):
+
+            temp = np.zeros((vorder + 1, 3))
+
+            # Nu.T^(k)*pts
+            for s in range(vorder + 1):
+                for r in range(uorder + 1):
+                    temp[s, :] += (
+                        utemp[k, r]
+                        * pts[(uspan - uorder + r) % nub, (vspan - vorder + s) % nvb, :]
+                    )
+
+            # ramaining derivative orders: dk + du <= dorder
+            dd = min(dorder - k, dv)
+
+            # .. * Nv^(l)
+            for l in range(dd + 1):
+                for s in range(vorder + 1):
+                    out[i, k, l, :] += vtemp[l, s] * temp[s, :]
+
+    return out
 
 
 #%% matrices

From a3ca1e298ba2109c8438c85994dfa1bf6e291557 Mon Sep 17 00:00:00 2001
From: adam-urbanczyk <13981538+adam-urbanczyk@users.noreply.github.com>
Date: Thu, 17 Jul 2025 16:11:22 +0200
Subject: [PATCH 41/55] Roundtrip fixes

---
 cadquery/occ_impl/nurbs.py |  9 +++++++++
 tests/test_nurbs.py        | 13 +++++++++++++
 2 files changed, 22 insertions(+)

diff --git a/cadquery/occ_impl/nurbs.py b/cadquery/occ_impl/nurbs.py
index f2389e845..f6b36c46b 100644
--- a/cadquery/occ_impl/nurbs.py
+++ b/cadquery/occ_impl/nurbs.py
@@ -150,6 +150,9 @@ def fromEdge(cls, e: Edge):
         order = g.Degree()
         periodic = g.IsPeriodic()
 
+        if periodic:
+            knots = knots[order:-order]
+
         return cls(pts, knots, order, periodic)
 
     def __call__(self, us: Array) -> Array:
@@ -244,6 +247,12 @@ def fromFace(cls, f: Face):
         uperiodic = g.IsUPeriodic()
         vperiodic = g.IsVPeriodic()
 
+        if uperiodic:
+            uknots = uknots[uorder:-uorder]
+
+        if vperiodic:
+            vknots = vknots[vorder:-vorder]
+
         return cls(pts, uknots, vknots, uorder, vorder, uperiodic, vperiodic)
 
     def __call__(self, u: Array, v: Array) -> Array:
diff --git a/tests/test_nurbs.py b/tests/test_nurbs.py
index 20308831c..2ccab0c13 100644
--- a/tests/test_nurbs.py
+++ b/tests/test_nurbs.py
@@ -133,6 +133,12 @@ def test_curve():
 
     assert e2.isValid()
 
+    # check roundtrip
+    crv3 = Curve.fromEdge(e2)
+
+    assert np.allclose(crv2.knots, crv3.knots)
+    assert np.allclose(crv2.pts, crv3.pts)
+
 
 def test_surface():
 
@@ -147,6 +153,13 @@ def test_surface():
     assert f.isValid()
     assert f.Area() == approx(1)
 
+    # roundtrip
+    srf2 = Surface.fromFace(f)
+
+    assert np.allclose(srf.uknots, srf2.uknots)
+    assert np.allclose(srf.vknots, srf2.vknots)
+    assert np.allclose(srf.pts, srf2.pts)
+
 
 def test_approximate():
 

From b12d67789fc6f1e88c33861d009876b65fc43dba Mon Sep 17 00:00:00 2001
From: adam-urbanczyk <13981538+adam-urbanczyk@users.noreply.github.com>
Date: Fri, 18 Jul 2025 09:04:36 +0200
Subject: [PATCH 42/55] Rough version of reparametrize

---
 cadquery/occ_impl/nurbs.py | 41 +++++++++++++++++++++++++++++++++++++-
 1 file changed, 40 insertions(+), 1 deletion(-)

diff --git a/cadquery/occ_impl/nurbs.py b/cadquery/occ_impl/nurbs.py
index f6b36c46b..94360c7a6 100644
--- a/cadquery/occ_impl/nurbs.py
+++ b/cadquery/occ_impl/nurbs.py
@@ -164,7 +164,7 @@ def __call__(self, us: Array) -> Array:
     def der(self, us: NDArray, dorder: int) -> NDArray:
 
         return nbCurveDer(
-            np.atleast_1d(us), self.order, self.knots, self.pts, self.periodic
+            np.atleast_1d(us), self.order, dorder, self.knots, self.pts, self.periodic
         )
 
 
@@ -1603,6 +1603,45 @@ def loft(
     return rv
 
 
+def reparametrize(
+    *curves: Curve, n: int = 100, w1: float = 1, w2: float = 1e-2
+) -> List[Curve]:
+
+    from scipy.optimize import fmin_l_bfgs_b
+
+    n_curves = len(curves)
+
+    u0 = np.tile(np.linspace(0, 1, n, False), n_curves)
+
+    def cost(u: Array) -> float:
+
+        rv1 = 0
+        us = np.split(u, n_curves)
+
+        pts = []
+
+        for i, ui in enumerate(us):
+            pts.append(curves[i](ui))
+
+            # parametric distance between points on the same curve
+            rv1 += np.sum((ui[:-1] - ui[1:]) ** 2) + np.sum((u[0] + 1 - u[-1]) ** 2)
+
+        rv2 = 0
+
+        for p1, p2 in zip(pts, pts[1:]):
+
+            # geometric distance between points on adjecent curves
+            rv2 += np.sum((p1 - p2) ** 2)
+
+        return w1 * rv1 + w2 * rv2
+
+    usol, _, _ = fmin_l_bfgs_b(cost, u0, approx_grad=True)
+
+    us = np.split(usol, n_curves)
+
+    return periodicApproximate([crv(u) for crv, u in zip(curves, us)], lam=1e-3)
+
+
 #%% for removal?
 @njit
 def findSpan(v, knots):

From 03716b792d7cc332ba110137a7f522b73cc6a9fc Mon Sep 17 00:00:00 2001
From: adam-urbanczyk <13981538+adam-urbanczyk@users.noreply.github.com>
Date: Fri, 25 Jul 2025 10:29:09 +0200
Subject: [PATCH 43/55] Adding 2D design matrix

---
 cadquery/occ_impl/nurbs.py | 96 ++++++++++++++++++++++++++++++++------
 tests/test_nurbs.py        | 23 ++++++++-
 2 files changed, 104 insertions(+), 15 deletions(-)

diff --git a/cadquery/occ_impl/nurbs.py b/cadquery/occ_impl/nurbs.py
index 94360c7a6..566bb6b26 100644
--- a/cadquery/occ_impl/nurbs.py
+++ b/cadquery/occ_impl/nurbs.py
@@ -294,6 +294,33 @@ def der(self, u: Array, v: Array, dorder: int) -> Array:
 #%% basis functions
 
 
+@njiti
+def _preprocess(
+    u: Array, order: int, knots: Array, periodic: float
+) -> Tuple[Array, Array, Optional[int], Optional[int], int]:
+    """
+    Helper for handling peridocity. This function extends the knot vector,
+    wraps the parameters and calculates the delta span.
+    """
+
+    # handle periodicity
+    if periodic:
+        period = knots[-1] - knots[0]
+        u_ = u % period
+        knots_ext = extendKnots(order, knots)
+        minspan = 0
+        maxspan = len(knots) - 1
+        deltaspan = order - 1
+    else:
+        u_ = u
+        knots_ext = knots
+        minspan = None
+        maxspan = None
+        deltaspan = 0
+
+    return u_, knots_ext, minspan, maxspan, deltaspan
+
+
 @njiti
 def extendKnots(order: int, knots: Array) -> Array:
     """
@@ -546,20 +573,7 @@ def nbCurve(
     # number of control points
     nb = pts.shape[0]
 
-    # handle periodicity
-    if periodic:
-        period = knots[-1] - knots[0]
-        u_ = u % period
-        knots_ext = extendKnots(order, knots)
-        minspan = 0
-        maxspan = len(knots) - 1
-        deltaspan = order - 1
-    else:
-        u_ = u
-        knots_ext = knots
-        minspan = None
-        maxspan = None
-        deltaspan = 0
+    u_, knots_ext, minspan, maxspan, deltaspan = _preprocess(u, order, knots, periodic)
 
     # number of param values
     nu = np.size(u)
@@ -962,6 +976,60 @@ def designMatrix(u: Array, order: int, knots: Array) -> COO:
     return rv
 
 
+# @njit
+def designMatrix2D(
+    uv: Array, uorder: int, vorder: int, uknots: Array, vknots: Array
+) -> COO:
+    """
+    Create a sparse tensor product design matrix.
+    """
+
+    # number of param values
+    ni = uv.shape[0]
+
+    # chunck size
+    nu = uorder + 1
+    nv = vorder + 1
+    nj = nu * nv
+
+    # number of basis
+    nu_total = len(uknots) - uorder - 1
+    nv_total = len(vknots) - vorder - 1
+
+    # temp chunck storage
+    utemp = np.zeros(nu)
+    vtemp = np.zeros(nv)
+
+    # initialize the empty matrix
+    rv = COO(
+        i=np.empty(ni * nj, dtype=np.int64),
+        j=np.empty(ni * nj, dtype=np.int64),
+        v=np.empty(ni * nj),
+    )
+
+    # loop over param values
+    for i in range(ni):
+        ui, vi = uv[i, :]
+
+        # find the supporting span
+        uspan = nbFindSpan(ui, uorder, uknots)
+        vspan = nbFindSpan(vi, vorder, vknots)
+
+        # evaluate non-zero functions
+        nbBasis(uspan, ui, uorder, uknots, utemp)
+        nbBasis(vspan, vi, vorder, vknots, vtemp)
+
+        # update the matrix
+        rv.i[i * nj : (i + 1) * nj] = i
+        rv.j[i * nj : (i + 1) * nj] = (
+            (uspan - uorder + np.arange(nu)) * nv_total
+            + (vspan - vorder + np.arange(nv))[:, np.newaxis]
+        ).ravel()
+        rv.v[i * nj : (i + 1) * nj] = (utemp * vtemp[:, np.newaxis]).ravel()
+
+    return rv
+
+
 @njit
 def periodicDesignMatrix(u: Array, order: int, knots: Array) -> COO:
     """
diff --git a/tests/test_nurbs.py b/tests/test_nurbs.py
index 2ccab0c13..31b561abd 100644
--- a/tests/test_nurbs.py
+++ b/tests/test_nurbs.py
@@ -1,6 +1,7 @@
 from cadquery.occ_impl.nurbs import (
     designMatrix,
     periodicDesignMatrix,
+    designMatrix2D,
     nbFindSpan,
     nbBasis,
     nbBasisDer,
@@ -17,7 +18,7 @@
 import numpy as np
 import scipy.sparse as sp
 
-from pytest import approx, fixture
+from pytest import approx, fixture, mark
 
 
 @fixture
@@ -63,6 +64,26 @@ def test_periodic_dm():
     assert C.shape[1] == len(knots) - 1
 
 
+def test_dm_2d():
+
+    uknots = np.array([0, 0, 0, 0, 0.25, 0.5, 0.75, 1, 1, 1, 1])
+    uparams = np.linspace(0, 1, 100)
+    uorder = 3
+
+    vknots = np.array([0, 0, 0, 0.5, 1, 1, 1])
+    vparams = np.linspace(0, 1, 100)
+    vorder = 2
+
+    params = np.column_stack((uparams, vparams))
+
+    res = designMatrix2D(params, uorder, vorder, uknots, vknots)
+
+    C = res.coo()
+
+    assert C.shape[0] == len(uparams)
+    assert C.shape[1] == (len(uknots) - uorder - 1) * (len(vknots) - vorder - 1)
+
+
 def test_dm():
 
     knots = np.array([0, 0, 0, 0, 0.25, 0.5, 0.75, 1, 1, 1, 1])

From 0db4c093e52e7c4d2c540a57fd234606b81a67ae Mon Sep 17 00:00:00 2001
From: adam-urbanczyk <13981538+adam-urbanczyk@users.noreply.github.com>
Date: Fri, 25 Jul 2025 11:04:47 +0200
Subject: [PATCH 44/55] Use _preprocess everywhere

---
 cadquery/occ_impl/nurbs.py | 80 +++++++-------------------------------
 1 file changed, 13 insertions(+), 67 deletions(-)

diff --git a/cadquery/occ_impl/nurbs.py b/cadquery/occ_impl/nurbs.py
index 566bb6b26..2cf40b138 100644
--- a/cadquery/occ_impl/nurbs.py
+++ b/cadquery/occ_impl/nurbs.py
@@ -636,19 +636,7 @@ def nbCurveDer(
     nb = pts.shape[0]
 
     # handle periodicity
-    if periodic:
-        period = knots[-1] - knots[0]
-        u_ = u % period
-        knots_ext = extendKnots(order, knots)
-        minspan = 0
-        maxspan = len(knots) - 1
-        deltaspan = order - 1
-    else:
-        u_ = u
-        knots_ext = knots
-        minspan = None
-        maxspan = None
-        deltaspan = 0
+    u_, knots_ext, minspan, maxspan, deltaspan = _preprocess(u, order, knots, periodic)
 
     # number of param values
     nu = np.size(u)
@@ -727,33 +715,12 @@ def nbSurface(
     nvb = pts.shape[1]
 
     # handle periodicity
-    if uperiodic:
-        uperiod = uknots[-1] - uknots[0]
-        u_ = u % uperiod
-        uknots_ext = extendKnots(uorder, uknots)
-        minspanu = 0
-        maxspanu = len(uknots) - 1
-        deltaspanu = uorder - 1
-    else:
-        u_ = u
-        uknots_ext = uknots
-        minspanu = None
-        maxspanu = None
-        deltaspanu = 0
-
-    if vperiodic:
-        vperiod = vknots[-1] - vknots[0]
-        v_ = v % vperiod
-        vknots_ext = extendKnots(vorder, vknots)
-        minspanv = 0
-        maxspanv = len(vknots) - 1
-        deltaspanv = vorder - 1
-    else:
-        v_ = v
-        vknots_ext = vknots
-        minspanv = None
-        maxspanv = None
-        deltaspanv = 0
+    u_, uknots_ext, minspanu, maxspanu, deltaspanu = _preprocess(
+        u, uorder, uknots, uperiodic
+    )
+    v_, vknots_ext, minspanv, maxspanv, deltaspanv = _preprocess(
+        v, vorder, vknots, vperiodic
+    )
 
     # number of param values
     nu = np.size(u)
@@ -855,33 +822,12 @@ def nbSurfaceDer(
     nvb = pts.shape[1]
 
     # handle periodicity
-    if uperiodic:
-        uperiod = uknots[-1] - uknots[0]
-        u_ = u % uperiod
-        uknots_ext = extendKnots(uorder, uknots)
-        minspanu = 0
-        maxspanu = len(uknots) - 1
-        deltaspanu = uorder - 1
-    else:
-        u_ = u
-        uknots_ext = uknots
-        minspanu = None
-        maxspanu = None
-        deltaspanu = 0
-
-    if vperiodic:
-        vperiod = vknots[-1] - vknots[0]
-        v_ = v % vperiod
-        vknots_ext = extendKnots(vorder, vknots)
-        minspanv = 0
-        maxspanv = len(vknots) - 1
-        deltaspanv = vorder - 1
-    else:
-        v_ = v
-        vknots_ext = vknots
-        minspanv = None
-        maxspanv = None
-        deltaspanv = 0
+    u_, uknots_ext, minspanu, maxspanu, deltaspanu = _preprocess(
+        u, uorder, uknots, uperiodic
+    )
+    v_, vknots_ext, minspanv, maxspanv, deltaspanv = _preprocess(
+        v, vorder, vknots, vperiodic
+    )
 
     # number of param values
     nu = np.size(u)

From a90756f5f6d0c1c31b4b7d5e491f15b7c707abbf Mon Sep 17 00:00:00 2001
From: adam-urbanczyk <13981538+adam-urbanczyk@users.noreply.github.com>
Date: Fri, 25 Jul 2025 12:39:12 +0200
Subject: [PATCH 45/55] Reparametrize fixes

---
 cadquery/occ_impl/nurbs.py | 14 +++++++++-----
 1 file changed, 9 insertions(+), 5 deletions(-)

diff --git a/cadquery/occ_impl/nurbs.py b/cadquery/occ_impl/nurbs.py
index 2cf40b138..9eccca1a2 100644
--- a/cadquery/occ_impl/nurbs.py
+++ b/cadquery/occ_impl/nurbs.py
@@ -1618,14 +1618,18 @@ def loft(
 
 
 def reparametrize(
-    *curves: Curve, n: int = 100, w1: float = 1, w2: float = 1e-2
+    *curves: Curve, n: int = 100, w1: float = 1, w2: float = 1e-1
 ) -> List[Curve]:
 
     from scipy.optimize import fmin_l_bfgs_b
 
     n_curves = len(curves)
 
-    u0 = np.tile(np.linspace(0, 1, n, False), n_curves)
+    u0_0 = np.linspace(0, 1, n, False)
+    u0 = np.tile(u0_0, n_curves)
+
+    # scaling for the second cost term
+    scale = n * np.linalg.norm(curves[0](u0[0]) - curves[1](u0[n]))
 
     def cost(u: Array) -> float:
 
@@ -1638,14 +1642,14 @@ def cost(u: Array) -> float:
             pts.append(curves[i](ui))
 
             # parametric distance between points on the same curve
-            rv1 += np.sum((ui[:-1] - ui[1:]) ** 2) + np.sum((u[0] + 1 - u[-1]) ** 2)
+            rv1 += np.sum((ui[:-1] - ui[1:]) ** 2) + np.sum((ui[0] + 1 - ui[-1]) ** 2)
 
         rv2 = 0
 
         for p1, p2 in zip(pts, pts[1:]):
 
             # geometric distance between points on adjecent curves
-            rv2 += np.sum((p1 - p2) ** 2)
+            rv2 += np.sum(((p1 - p2) / scale) ** 2)
 
         return w1 * rv1 + w2 * rv2
 
@@ -1653,7 +1657,7 @@ def cost(u: Array) -> float:
 
     us = np.split(usol, n_curves)
 
-    return periodicApproximate([crv(u) for crv, u in zip(curves, us)], lam=1e-3)
+    return periodicApproximate([crv(u) for crv, u in zip(curves, us)], knots=n, lam=0)
 
 
 #%% for removal?

From 3461b4ac2aa5a6135df0068a1d60e566b7c63f20 Mon Sep 17 00:00:00 2001
From: adam-urbanczyk <13981538+adam-urbanczyk@users.noreply.github.com>
Date: Fri, 25 Jul 2025 14:50:19 +0200
Subject: [PATCH 46/55] Unify periodic design matrix handling

---
 cadquery/occ_impl/nurbs.py | 101 +++++++++++++++----------------------
 1 file changed, 42 insertions(+), 59 deletions(-)

diff --git a/cadquery/occ_impl/nurbs.py b/cadquery/occ_impl/nurbs.py
index 9eccca1a2..0be0cd373 100644
--- a/cadquery/occ_impl/nurbs.py
+++ b/cadquery/occ_impl/nurbs.py
@@ -314,8 +314,8 @@ def _preprocess(
     else:
         u_ = u
         knots_ext = knots
-        minspan = None
-        maxspan = None
+        minspan = order
+        maxspan = knots.shape[0] - order - 1
         deltaspan = 0
 
     return u_, knots_ext, minspan, maxspan, deltaspan
@@ -883,13 +883,23 @@ def nbSurfaceDer(
 
 
 @njit
-def designMatrix(u: Array, order: int, knots: Array) -> COO:
+def designMatrix(u: Array, order: int, knots: Array, periodic: bool = False) -> COO:
     """
-    Create a sparse design matrix.
+    Create a sparse (possibly periodic) design matrix.
     """
 
+    # extend the knots
+    knots_ext = np.concat(
+        (knots[-order:-1] - knots[-1], knots, knots[-1] + knots[1:order])
+    )
+
+    u_, knots_ext, minspan, maxspan, deltaspan = _preprocess(u, order, knots, periodic)
+
     # number of param values
-    nu = np.size(u)
+    nu = len(u)
+
+    # number of basis functions
+    nb = maxspan
 
     # chunck size
     n = order + 1
@@ -906,17 +916,19 @@ def designMatrix(u: Array, order: int, knots: Array) -> COO:
 
     # loop over param values
     for i in range(nu):
-        ui = u[i]
+        ui = u_[i]
 
         # find the supporting span
-        span = nbFindSpan(ui, order, knots)
+        span = nbFindSpan(ui, order, knots, minspan, maxspan) + deltaspan
 
         # evaluate non-zero functions
-        nbBasis(span, ui, order, knots, temp)
+        nbBasis(span, ui, order, knots_ext, temp)
 
         # update the matrix
         rv.i[i * n : (i + 1) * n] = i
-        rv.j[i * n : (i + 1) * n] = span - order + np.arange(n)
+        rv.j[i * n : (i + 1) * n] = (
+            span - order + np.arange(n)
+        ) % nb  # NB: this is due to peridicity
         rv.v[i * n : (i + 1) * n] = temp
 
     return rv
@@ -924,12 +936,25 @@ def designMatrix(u: Array, order: int, knots: Array) -> COO:
 
 # @njit
 def designMatrix2D(
-    uv: Array, uorder: int, vorder: int, uknots: Array, vknots: Array
+    uv: Array,
+    uorder: int,
+    vorder: int,
+    uknots: Array,
+    vknots: Array,
+    uperiodic: bool = False,
+    vperiodic: bool = False,
 ) -> COO:
     """
     Create a sparse tensor product design matrix.
     """
 
+    u_, uknots_ext, minspanu, maxspanu, deltaspanu = _preprocess(
+        uv[:, 0], uorder, uknots, uperiodic
+    )
+    v_, vknots_ext, minspanv, maxspanv, deltaspanv = _preprocess(
+        uv[:, 1], vorder, vknots, vperiodic
+    )
+
     # number of param values
     ni = uv.shape[0]
 
@@ -939,8 +964,7 @@ def designMatrix2D(
     nj = nu * nv
 
     # number of basis
-    nu_total = len(uknots) - uorder - 1
-    nv_total = len(vknots) - vorder - 1
+    nv_total = maxspanv
 
     # temp chunck storage
     utemp = np.zeros(nu)
@@ -955,15 +979,15 @@ def designMatrix2D(
 
     # loop over param values
     for i in range(ni):
-        ui, vi = uv[i, :]
+        ui, vi = u_[i], v_[i]
 
         # find the supporting span
-        uspan = nbFindSpan(ui, uorder, uknots)
-        vspan = nbFindSpan(vi, vorder, vknots)
+        uspan = nbFindSpan(ui, uorder, uknots, minspanu, maxspanu) + deltaspanu
+        vspan = nbFindSpan(vi, vorder, vknots, minspanv, maxspanv) + deltaspanv
 
         # evaluate non-zero functions
-        nbBasis(uspan, ui, uorder, uknots, utemp)
-        nbBasis(vspan, vi, vorder, vknots, vtemp)
+        nbBasis(uspan, ui, uorder, uknots_ext, utemp)
+        nbBasis(vspan, vi, vorder, vknots_ext, vtemp)
 
         # update the matrix
         rv.i[i * nj : (i + 1) * nj] = i
@@ -982,48 +1006,7 @@ def periodicDesignMatrix(u: Array, order: int, knots: Array) -> COO:
     Create a sparse periodic design matrix.
     """
 
-    # extend the knots
-    knots_ext = np.concat(
-        (knots[-order:-1] - knots[-1], knots, knots[-1] + knots[1:order])
-    )
-
-    # number of param values
-    nu = len(u)
-
-    # number of basis functions
-    nb = len(knots) - 1
-
-    # chunck size
-    n = order + 1
-
-    # temp chunck storage
-    temp = np.zeros(n)
-
-    # initialize the empty matrix
-    rv = COO(
-        i=np.empty(n * nu, dtype=np.int64),
-        j=np.empty(n * nu, dtype=np.int64),
-        v=np.empty(n * nu),
-    )
-
-    # loop over param values
-    for i in range(nu):
-        ui = u[i]
-
-        # find the supporting span
-        span = nbFindSpan(ui, order, knots, 0, nb) + order - 1
-
-        # evaluate non-zero functions
-        nbBasis(span, ui, order, knots_ext, temp)
-
-        # update the matrix
-        rv.i[i * n : (i + 1) * n] = i
-        rv.j[i * n : (i + 1) * n] = (
-            span - order + np.arange(n)
-        ) % nb  # NB: this is due to peridicity
-        rv.v[i * n : (i + 1) * n] = temp
-
-    return rv
+    return designMatrix(u, order, knots, periodic=True)
 
 
 @njit

From 05b5fe38a5cf4f3d1f18b1aabc3ef7d3e97ad97c Mon Sep 17 00:00:00 2001
From: adam-urbanczyk <13981538+adam-urbanczyk@users.noreply.github.com>
Date: Fri, 25 Jul 2025 16:18:43 +0200
Subject: [PATCH 47/55] 2D dm fix

---
 cadquery/occ_impl/nurbs.py | 5 +++--
 1 file changed, 3 insertions(+), 2 deletions(-)

diff --git a/cadquery/occ_impl/nurbs.py b/cadquery/occ_impl/nurbs.py
index 0be0cd373..61fe863a3 100644
--- a/cadquery/occ_impl/nurbs.py
+++ b/cadquery/occ_impl/nurbs.py
@@ -964,6 +964,7 @@ def designMatrix2D(
     nj = nu * nv
 
     # number of basis
+    nu_total = maxspanu
     nv_total = maxspanv
 
     # temp chunck storage
@@ -992,8 +993,8 @@ def designMatrix2D(
         # update the matrix
         rv.i[i * nj : (i + 1) * nj] = i
         rv.j[i * nj : (i + 1) * nj] = (
-            (uspan - uorder + np.arange(nu)) * nv_total
-            + (vspan - vorder + np.arange(nv))[:, np.newaxis]
+            ((uspan - uorder + np.arange(nu)) % nu_total) * nv_total
+            + ((vspan - vorder + np.arange(nv)) % nv_total)[:, np.newaxis]
         ).ravel()
         rv.v[i * nj : (i + 1) * nj] = (utemp * vtemp[:, np.newaxis]).ravel()
 

From 8b5082c05ecfdcb5fc9eea5015f03fd5f3010d3d Mon Sep 17 00:00:00 2001
From: adam-urbanczyk <13981538+adam-urbanczyk@users.noreply.github.com>
Date: Fri, 1 Aug 2025 12:44:10 +0200
Subject: [PATCH 48/55] Add gradient to reparametrize

---
 cadquery/occ_impl/nurbs.py | 56 ++++++++++++++++++++++++++++++++++++--
 1 file changed, 54 insertions(+), 2 deletions(-)

diff --git a/cadquery/occ_impl/nurbs.py b/cadquery/occ_impl/nurbs.py
index 61fe863a3..3462c4c65 100644
--- a/cadquery/occ_impl/nurbs.py
+++ b/cadquery/occ_impl/nurbs.py
@@ -19,7 +19,6 @@
 )
 from OCP.gp import gp_Pnt
 from OCP.BRepBuilderAPI import BRepBuilderAPI_MakeEdge, BRepBuilderAPI_MakeFace
-from OCP.Geom import Geom_BSplineSurface
 
 from .shapes import Face, Edge
 
@@ -1623,6 +1622,8 @@ def cost(u: Array) -> float:
         pts = []
 
         for i, ui in enumerate(us):
+
+            # evaluate
             pts.append(curves[i](ui))
 
             # parametric distance between points on the same curve
@@ -1637,7 +1638,58 @@ def cost(u: Array) -> float:
 
         return w1 * rv1 + w2 * rv2
 
-    usol, _, _ = fmin_l_bfgs_b(cost, u0, approx_grad=True)
+    def grad(u: Array) -> Array:
+
+        rv1 = np.zeros_like(u)
+        us = np.split(u, n_curves)
+
+        pts = []
+        tgts = []
+
+        for i, ui in enumerate(us):
+
+            # evaluate up to 1st derivative
+            tmp = curves[i].der(ui, 1)
+
+            pts.append(tmp[:, 0, :].squeeze())
+            tgts.append(tmp[:, 1, :].squeeze())
+
+            # parametric distance between points on the same curve
+            delta = np.roll(ui, -1) - ui
+            delta[-1] += 1
+            delta *= -2
+            delta -= np.roll(delta, 1)
+
+            rv1[i * n : (i + 1) * n] = delta
+
+        rv2 = np.zeros_like(u)
+
+        for i, _ in enumerate(us):
+            # geometric distance between points on adjecent curves
+
+            # first profile
+            if i == 0:
+                p1, p2, t = pts[i], pts[i + 1], tgts[i]
+
+                rv2[i * n : (i + 1) * n] = (2 / scale ** 2 * (p1 - p2) * t).sum(1)
+
+            # middle profile
+            elif i + 1 < n_curves:
+                p1, p2, t = pts[i], pts[i + 1], tgts[i]
+                p0 = pts[i - 1]
+
+                rv2[i * n : (i + 1) * n] = (2 / scale ** 2 * (p1 - p2) * t).sum(1)
+                rv2[i * n : (i + 1) * n] += (-2 / scale ** 2 * (p0 - p1) * t).sum(1)
+
+            # last profile
+            else:
+                p1, p2, t = pts[i - 1], pts[i], tgts[i]
+
+                rv2[i * n : (i + 1) * n] = (-2 / scale ** 2 * (p1 - p2) * t).sum(1)
+
+        return w1 * rv1 + w2 * rv2
+
+    usol, _, _ = fmin_l_bfgs_b(cost, u0, grad)
 
     us = np.split(usol, n_curves)
 

From 64a2eaa4a6886dacb462719517036c097a3087b9 Mon Sep 17 00:00:00 2001
From: adam-urbanczyk <13981538+adam-urbanczyk@users.noreply.github.com>
Date: Fri, 1 Aug 2025 14:52:35 +0200
Subject: [PATCH 49/55] Add test for reparametrize

---
 tests/test_nurbs.py | 31 +++++++++++++++++++++++++++++++
 1 file changed, 31 insertions(+)

diff --git a/tests/test_nurbs.py b/tests/test_nurbs.py
index 31b561abd..19931463a 100644
--- a/tests/test_nurbs.py
+++ b/tests/test_nurbs.py
@@ -11,6 +11,7 @@
     periodicApproximate,
     periodicLoft,
     loft,
+    reparametrize,
 )
 
 from cadquery.func import circle
@@ -50,6 +51,18 @@ def trimmed_circles() -> list[Curve]:
     return cs
 
 
+@fixture
+def rotated_circles() -> list[Curve]:
+
+    pts1 = np.array([v.toTuple() for v in circle(1).sample(100)[0]])
+    pts2 = np.array([v.toTuple() for v in circle(1).moved(z=1, rz=90).sample(100)[0]])
+
+    c1 = periodicApproximate(pts1)
+    c2 = periodicApproximate(pts2)
+
+    return [c1, c2]
+
+
 def test_periodic_dm():
 
     knots = np.linspace(0, 1, 5)
@@ -245,3 +258,21 @@ def test_loft(circles, trimmed_circles):
     surf2 = loft(*trimmed_circles)
 
     assert surf2.face().isValid()
+
+
+def test_reparametrize(rotated_circles):
+
+    c1, c2 = rotated_circles
+
+    # this surface will be twisted
+    surf = loft(c1, c2, order=2, lam=1e-6)
+
+    # this should adjust the paramatrizations
+    c1r, c2r = reparametrize(c1, c2)
+
+    # resulting loft should not be twisted
+    surfr = loft(c1r, c2r, order=2, lam=1e-6)
+
+    # assert that the surface is indeed not twisted
+    assert surfr.face().Area() == approx(2 * np.pi, 1e-3)
+    assert surfr.face().Area() >= 1.01 * surf.face().Area()

From 57af33b6dd5288b1c5ad9034d42056a63e78a184 Mon Sep 17 00:00:00 2001
From: adam-urbanczyk <13981538+adam-urbanczyk@users.noreply.github.com>
Date: Fri, 8 Aug 2025 13:51:15 +0200
Subject: [PATCH 50/55] Enable oversampling in reparametrize

---
 cadquery/occ_impl/nurbs.py | 22 ++++++++++++----------
 1 file changed, 12 insertions(+), 10 deletions(-)

diff --git a/cadquery/occ_impl/nurbs.py b/cadquery/occ_impl/nurbs.py
index 3462c4c65..8363918ca 100644
--- a/cadquery/occ_impl/nurbs.py
+++ b/cadquery/occ_impl/nurbs.py
@@ -1,4 +1,4 @@
-#%% imports
+# %% imports
 import numpy as np
 import scipy.sparse as sp
 
@@ -31,7 +31,7 @@
 )
 
 
-#%% internal helpers
+# %% internal helpers
 
 
 def _colPtsArray(pts: NDArray) -> TColgp_Array1OfPnt:
@@ -79,7 +79,7 @@ def _colIntArray(knots: NDArray) -> TColStd_Array1OfInteger:
     return rv
 
 
-#%% vocabulary types
+# %% vocabulary types
 
 Array = ndarray  # NDArray[np.floating]
 ArrayI = ndarray  # NDArray[np.int_]
@@ -290,7 +290,7 @@ def der(self, u: Array, v: Array, dorder: int) -> Array:
         )
 
 
-#%% basis functions
+# %% basis functions
 
 
 @njiti
@@ -539,7 +539,7 @@ def nbBasisDer(i: int, u: float, order: int, dorder: int, knots: Array, out: Arr
         r *= order - k
 
 
-#%% evaluation
+# %% evaluation
 
 
 @njit
@@ -878,7 +878,7 @@ def nbSurfaceDer(
     return out
 
 
-#%% matrices
+# %% matrices
 
 
 @njit
@@ -1249,7 +1249,7 @@ def discretePenalty(us: Array, order: int, splineorder: int = 3) -> COO:
     return rv
 
 
-#%% construction
+# %% construction
 
 
 @multidispatch
@@ -1601,7 +1601,7 @@ def loft(
 
 
 def reparametrize(
-    *curves: Curve, n: int = 100, w1: float = 1, w2: float = 1e-1
+    *curves: Curve, n: int = 100, knots: int = 100, w1: float = 1, w2: float = 1e-1
 ) -> List[Curve]:
 
     from scipy.optimize import fmin_l_bfgs_b
@@ -1693,10 +1693,12 @@ def grad(u: Array) -> Array:
 
     us = np.split(usol, n_curves)
 
-    return periodicApproximate([crv(u) for crv, u in zip(curves, us)], knots=n, lam=0)
+    return periodicApproximate(
+        [crv(u) for crv, u in zip(curves, us)], knots=knots, lam=0
+    )
 
 
-#%% for removal?
+# %% for removal?
 @njit
 def findSpan(v, knots):
 

From 820efcc309edd1515ea1caebee93ca8c15d90309 Mon Sep 17 00:00:00 2001
From: adam-urbanczyk <13981538+adam-urbanczyk@users.noreply.github.com>
Date: Fri, 8 Aug 2025 16:05:15 +0200
Subject: [PATCH 51/55] Adding approximate 2D

---
 cadquery/occ_impl/nurbs.py | 41 +++++++++++++++++++++++++++++++++++++-
 1 file changed, 40 insertions(+), 1 deletion(-)

diff --git a/cadquery/occ_impl/nurbs.py b/cadquery/occ_impl/nurbs.py
index 8363918ca..36429e8a5 100644
--- a/cadquery/occ_impl/nurbs.py
+++ b/cadquery/occ_impl/nurbs.py
@@ -933,7 +933,7 @@ def designMatrix(u: Array, order: int, knots: Array, periodic: bool = False) ->
     return rv
 
 
-# @njit
+@njit
 def designMatrix2D(
     uv: Array,
     uorder: int,
@@ -1532,6 +1532,45 @@ def approximate(
     return rv
 
 
+def approximate2D(
+    data: Array,
+    uv: Array,
+    uorder: int,
+    vorder: int,
+    uknots: int | Array = 50,
+    vknots: int | Array = 50,
+    uperiodic: bool = False,
+    vperiodic: bool = False,
+) -> Surface:
+    """
+    Simple 2D surface approximation (without any penalty).
+    """
+
+    # process the knots
+    uknots_ = uknots if isinstance(uknots, Array) else np.linspace(0, 1, uknots)
+    vknots_ = vknots if isinstance(vknots, Array) else np.linspace(0, 1, vknots)
+
+    # create the desing matrix
+    C = designMatrix2D(uv, uorder, vorder, uknots_, vknots_, uperiodic, vperiodic).csc()
+
+    # solve normal equations
+    D, L, P = ldl(C.T @ C, False)
+    pts = ldl_solve(C.T @ data, D, L, P).toarray()
+
+    # construt the result
+    rv = Surface(
+        pts.reshape((len(uknots_) - int(uperiodic), len(vknots_) - int(vperiodic), 3)),
+        uknots_,
+        vknots_,
+        uorder,
+        vorder,
+        uperiodic,
+        vperiodic,
+    )
+
+    return rv
+
+
 def periodicLoft(*curves: Curve, order: int = 3) -> Surface:
 
     nknots: int = len(curves) + 1

From 93e0e998f718ff7f90214365f14a7442468eb0ae Mon Sep 17 00:00:00 2001
From: adam-urbanczyk <13981538+adam-urbanczyk@users.noreply.github.com>
Date: Mon, 11 Aug 2025 07:47:13 +0200
Subject: [PATCH 52/55] Add penalty matrix 2D

---
 cadquery/occ_impl/nurbs.py | 84 ++++++++++++++++++++++++++++++++++++--
 1 file changed, 80 insertions(+), 4 deletions(-)

diff --git a/cadquery/occ_impl/nurbs.py b/cadquery/occ_impl/nurbs.py
index 36429e8a5..4a3fa6a29 100644
--- a/cadquery/occ_impl/nurbs.py
+++ b/cadquery/occ_impl/nurbs.py
@@ -888,10 +888,6 @@ def designMatrix(u: Array, order: int, knots: Array, periodic: bool = False) ->
     """
 
     # extend the knots
-    knots_ext = np.concat(
-        (knots[-order:-1] - knots[-1], knots, knots[-1] + knots[1:order])
-    )
-
     u_, knots_ext, minspan, maxspan, deltaspan = _preprocess(u, order, knots, periodic)
 
     # number of param values
@@ -947,6 +943,7 @@ def designMatrix2D(
     Create a sparse tensor product design matrix.
     """
 
+    # extend the knots and preprocess
     u_, uknots_ext, minspanu, maxspanu, deltaspanu = _preprocess(
         uv[:, 0], uorder, uknots, uperiodic
     )
@@ -1249,6 +1246,85 @@ def discretePenalty(us: Array, order: int, splineorder: int = 3) -> COO:
     return rv
 
 
+@njit
+def penaltyMatrix2D(
+    uv: Array,
+    uorder: int,
+    vorder: int,
+    dorder: int,
+    uknots: Array,
+    vknots: Array,
+    uperiodic: bool = False,
+    vperiodic: bool = False,
+) -> list[COO]:
+    """
+    Create sparse tensor product 2D derivative matrices.
+    """
+
+    # extend the knots and preprocess
+    u_, uknots_ext, minspanu, maxspanu, deltaspanu = _preprocess(
+        uv[:, 0], uorder, uknots, uperiodic
+    )
+    v_, vknots_ext, minspanv, maxspanv, deltaspanv = _preprocess(
+        uv[:, 1], vorder, vknots, vperiodic
+    )
+
+    # number of param values
+    ni = uv.shape[0]
+
+    # chunck size
+    nu = uorder + 1
+    nv = vorder + 1
+    nj = nu * nv
+
+    # number of basis
+    nu_total = maxspanu
+    nv_total = maxspanv
+
+    # temp chunck storage
+    utemp = np.zeros((dorder + 1, nu))
+    vtemp = np.zeros((dorder + 1, nv))
+
+    # initialize the emptry matrices
+    rv = []
+    for i in range(dorder + 1):
+        rv.append(
+            COO(
+                i=np.empty(ni * nj, dtype=np.int64),
+                j=np.empty(ni * nj, dtype=np.int64),
+                v=np.empty(ni * nj),
+            )
+        )
+
+    # loop over param values
+    for i in range(ni):
+        ui, vi = u_[i], v_[i]
+
+        # find the supporting span
+        uspan = nbFindSpan(ui, uorder, uknots, minspanu, maxspanu) + deltaspanu
+        vspan = nbFindSpan(vi, vorder, vknots, minspanv, maxspanv) + deltaspanv
+
+        # evaluate non-zero functions
+        nbBasisDer(uspan, ui, uorder, dorder, uknots_ext, utemp)
+        nbBasisDer(vspan, vi, vorder, dorder, vknots_ext, vtemp)
+
+        # update the matrices - iterate over all derivative paris
+        for dv in range(dorder + 1):
+
+            du = dorder - dv  # NB: du + dv == dorder
+
+            rv[dv].i[i * nj : (i + 1) * nj] = i
+            rv[dv].j[i * nj : (i + 1) * nj] = (
+                ((uspan - uorder + np.arange(nu)) % nu_total) * nv_total
+                + ((vspan - vorder + np.arange(nv)) % nv_total)[:, np.newaxis]
+            ).ravel()
+            rv[dv].v[i * nj : (i + 1) * nj] = (
+                utemp[du, :] * vtemp[dv, :, np.newaxis]
+            ).ravel()
+
+    return rv
+
+
 # %% construction
 
 

From 5b66418952599cb51ce0987b5959ceb01c590b61 Mon Sep 17 00:00:00 2001
From: adam-urbanczyk <13981538+adam-urbanczyk@users.noreply.github.com>
Date: Wed, 13 Aug 2025 08:02:21 +0200
Subject: [PATCH 53/55] Add penalties to approximate2D

---
 cadquery/occ_impl/nurbs.py | 32 +++++++++++++++++++++++++++++++-
 1 file changed, 31 insertions(+), 1 deletion(-)

diff --git a/cadquery/occ_impl/nurbs.py b/cadquery/occ_impl/nurbs.py
index 4a3fa6a29..1ae55f184 100644
--- a/cadquery/occ_impl/nurbs.py
+++ b/cadquery/occ_impl/nurbs.py
@@ -6,6 +6,8 @@
 
 from typing import NamedTuple, Optional, Tuple, List, Union, cast
 
+from math import comb
+
 from numpy.typing import NDArray
 from numpy import linspace, ndarray
 
@@ -1617,6 +1619,8 @@ def approximate2D(
     vknots: int | Array = 50,
     uperiodic: bool = False,
     vperiodic: bool = False,
+    penalty: int = 3,
+    lam: float = 0,
 ) -> Surface:
     """
     Simple 2D surface approximation (without any penalty).
@@ -1629,8 +1633,34 @@ def approximate2D(
     # create the desing matrix
     C = designMatrix2D(uv, uorder, vorder, uknots_, vknots_, uperiodic, vperiodic).csc()
 
+    # handle penalties if requested
+    if lam:
+        # construct the penalty grid
+        Up, Vp = np.meshgrid(
+            np.linspace(uknots_[0], uknots_[-1], 2 * len(uknots_) * uorder),
+            np.linspace(vknots_[0], vknots_[-1], 2 * len(vknots_) * vorder),
+        )
+        up = Up.ravel()
+        vp = Vp.ravel()
+        uvp = np.column_stack((up, vp))
+
+        # construct the derivative matrices
+        penalties = penaltyMatrix2D(
+            uvp, uorder, vorder, penalty, uknots_, vknots_, uperiodic, vperiodic,
+        )
+
+        # augment the design matrix
+        tmp = [comb(penalty, i) * penalties[i].csc() for i in range(penalty + 1)]
+        Lu = uknots_[-1] - uknots_[0]  # v lenght of the parametric domain
+        Lv = vknots_[-1] - vknots_[0]  # u lenght of the parametric domain
+        P = Lu * Lv / len(up) * sp.vstack(tmp)
+
+        CtC = C.T @ C + lam * P.T @ P
+    else:
+        CtC = C.T @ C
+
     # solve normal equations
-    D, L, P = ldl(C.T @ C, False)
+    D, L, P = ldl(CtC, False)
     pts = ldl_solve(C.T @ data, D, L, P).toarray()
 
     # construt the result

From 281fe44e30d5cacc175ee333fc96a8b54e7f9ddd Mon Sep 17 00:00:00 2001
From: adam-urbanczyk <13981538+adam-urbanczyk@users.noreply.github.com>
Date: Fri, 15 Aug 2025 15:52:37 +0200
Subject: [PATCH 54/55] Add fairing, normals and offset

---
 cadquery/occ_impl/nurbs.py | 169 ++++++++++++++++++++++++++++++++-----
 1 file changed, 149 insertions(+), 20 deletions(-)

diff --git a/cadquery/occ_impl/nurbs.py b/cadquery/occ_impl/nurbs.py
index 1ae55f184..cc8793e33 100644
--- a/cadquery/occ_impl/nurbs.py
+++ b/cadquery/occ_impl/nurbs.py
@@ -291,6 +291,21 @@ def der(self, u: Array, v: Array, dorder: int) -> Array:
             self.vperiodic,
         )
 
+    def normal(self, u: Array, v: Array) -> Tuple[Array, Array]:
+        """
+        Evaluate surface normals.
+        """
+
+        ders = self.der(u, v, 1)
+
+        du = ders[:, 1, 0, :].squeeze()
+        dv = ders[:, 0, 1, :].squeeze()
+
+        rv = np.atleast_2d(np.cross(du, dv))
+        rv /= np.linalg.norm(rv, axis=1)[:, np.newaxis]
+
+        return rv, ders[:, 0, 0, :].squeeze()
+
 
 # %% basis functions
 
@@ -933,7 +948,8 @@ def designMatrix(u: Array, order: int, knots: Array, periodic: bool = False) ->
 
 @njit
 def designMatrix2D(
-    uv: Array,
+    u: Array,
+    v: Array,
     uorder: int,
     vorder: int,
     uknots: Array,
@@ -947,14 +963,14 @@ def designMatrix2D(
 
     # extend the knots and preprocess
     u_, uknots_ext, minspanu, maxspanu, deltaspanu = _preprocess(
-        uv[:, 0], uorder, uknots, uperiodic
+        u, uorder, uknots, uperiodic
     )
     v_, vknots_ext, minspanv, maxspanv, deltaspanv = _preprocess(
-        uv[:, 1], vorder, vknots, vperiodic
+        v, vorder, vknots, vperiodic
     )
 
     # number of param values
-    ni = uv.shape[0]
+    ni = len(u)
 
     # chunck size
     nu = uorder + 1
@@ -1250,7 +1266,8 @@ def discretePenalty(us: Array, order: int, splineorder: int = 3) -> COO:
 
 @njit
 def penaltyMatrix2D(
-    uv: Array,
+    u: Array,
+    v: Array,
     uorder: int,
     vorder: int,
     dorder: int,
@@ -1265,14 +1282,14 @@ def penaltyMatrix2D(
 
     # extend the knots and preprocess
     u_, uknots_ext, minspanu, maxspanu, deltaspanu = _preprocess(
-        uv[:, 0], uorder, uknots, uperiodic
+        u, uorder, uknots, uperiodic
     )
     v_, vknots_ext, minspanv, maxspanv, deltaspanv = _preprocess(
-        uv[:, 1], vorder, vknots, vperiodic
+        v, vorder, vknots, vperiodic
     )
 
     # number of param values
-    ni = uv.shape[0]
+    ni = len(u)
 
     # chunck size
     nu = uorder + 1
@@ -1327,6 +1344,32 @@ def penaltyMatrix2D(
     return rv
 
 
+def uniformGrid(
+    uknots: Array,
+    vknots: Array,
+    uorder: int,
+    vorder: int,
+    uperiodic: bool,
+    vperiodic: bool,
+) -> Tuple[Array, Array]:
+    """
+    Create a uniform grid for evaluating penalties.
+    """
+
+    Up, Vp = np.meshgrid(
+        np.linspace(
+            uknots[0], uknots[-1], 2 * len(uknots) * uorder, endpoint=not uperiodic
+        ),
+        np.linspace(
+            vknots[0], vknots[-1], 2 * len(vknots) * vorder, endpoint=not vperiodic
+        ),
+    )
+    up = Up.ravel()
+    vp = Vp.ravel()
+
+    return up, vp
+
+
 # %% construction
 
 
@@ -1342,8 +1385,9 @@ def periodicApproximate(
 
     npts = data.shape[0]
 
-    # parametrize the points
-    us = linspace(0, 1, npts, endpoint=False)
+    # parametrize the points if needed
+    if us is None:
+        us = linspace(0, 1, npts, endpoint=False)
 
     # construct the knot vector
     if isinstance(knots, int):
@@ -1612,7 +1656,8 @@ def approximate(
 
 def approximate2D(
     data: Array,
-    uv: Array,
+    u: Array,
+    v: Array,
     uorder: int,
     vorder: int,
     uknots: int | Array = 50,
@@ -1631,22 +1676,18 @@ def approximate2D(
     vknots_ = vknots if isinstance(vknots, Array) else np.linspace(0, 1, vknots)
 
     # create the desing matrix
-    C = designMatrix2D(uv, uorder, vorder, uknots_, vknots_, uperiodic, vperiodic).csc()
+    C = designMatrix2D(
+        u, v, uorder, vorder, uknots_, vknots_, uperiodic, vperiodic
+    ).csc()
 
     # handle penalties if requested
     if lam:
         # construct the penalty grid
-        Up, Vp = np.meshgrid(
-            np.linspace(uknots_[0], uknots_[-1], 2 * len(uknots_) * uorder),
-            np.linspace(vknots_[0], vknots_[-1], 2 * len(vknots_) * vorder),
-        )
-        up = Up.ravel()
-        vp = Vp.ravel()
-        uvp = np.column_stack((up, vp))
+        up, vp = uniformGrid(uknots_, vknots_, uorder, vorder, uperiodic, vperiodic)
 
         # construct the derivative matrices
         penalties = penaltyMatrix2D(
-            uvp, uorder, vorder, penalty, uknots_, vknots_, uperiodic, vperiodic,
+            up, vp, uorder, vorder, penalty, uknots_, vknots_, uperiodic, vperiodic,
         )
 
         # augment the design matrix
@@ -1677,6 +1718,60 @@ def approximate2D(
     return rv
 
 
+def fairPenalty(surf: Surface, penalty: int, lam: float) -> Surface:
+    """
+    Penalty-based surface fairing.
+    """
+
+    uknots = surf.uknots
+    vknots = surf.vknots
+    pts = surf.pts.reshape((-1, 3))
+
+    # generate penalty grid
+    up, vp = uniformGrid(
+        uknots, vknots, surf.uorder, surf.vorder, surf.uperiodic, surf.vperiodic
+    )
+
+    # generate penalty matrix
+    penalties = penaltyMatrix2D(
+        up,
+        vp,
+        surf.uorder,
+        surf.vorder,
+        penalty,
+        surf.uknots,
+        surf.vknots,
+        surf.uperiodic,
+        surf.vperiodic,
+    )
+
+    tmp = [comb(penalty, i) * penalties[i].csc() for i in range(penalty + 1)]
+    Lu = uknots[-1] - uknots[0]  # v lenght of the parametric domain
+    Lv = vknots[-1] - vknots[0]  # u lenght of the parametric domain
+    P = Lu * Lv / len(up) * sp.vstack(tmp)
+
+    # form and solve normal equations
+    CtC = sp.identity(pts.shape[0]) + lam * P.T @ P
+
+    D, L, P = ldl(CtC, False)
+    pts_new = ldl_solve(pts, D, L, P).toarray()
+
+    # construt the result
+    rv = Surface(
+        pts_new.reshape(
+            (len(uknots) - int(surf.uperiodic), len(vknots) - int(surf.vperiodic), 3)
+        ),
+        uknots,
+        vknots,
+        surf.uorder,
+        surf.vorder,
+        surf.uperiodic,
+        surf.vperiodic,
+    )
+
+    return rv
+
+
 def periodicLoft(*curves: Curve, order: int = 3) -> Surface:
 
     nknots: int = len(curves) + 1
@@ -1843,6 +1938,40 @@ def grad(u: Array) -> Array:
     )
 
 
+def offset(surf: Surface, d: float, lam: float = 1e-3) -> Surface:
+    """
+    Simple approximate offset.
+    """
+
+    # construct the knot grid
+    U, V = np.meshgrid(
+        np.linspace(surf.uknots[0], surf.uknots[-1], surf.uorder * len(surf.uknots)),
+        np.linspace(surf.vknots[0], surf.vknots[-1], surf.vorder * len(surf.uknots)),
+    )
+
+    us = U.ravel()
+    vs = V.ravel()
+
+    # evaluate the normals
+    ns, pts = surf.normal(us, vs)
+
+    # move the control points
+    pts += d * ns
+
+    return approximate2D(
+        pts,
+        us,
+        vs,
+        surf.uorder,
+        surf.vorder,
+        surf.uknots,
+        surf.vknots,
+        surf.uperiodic,
+        surf.vperiodic,
+        lam=lam,
+    )
+
+
 # %% for removal?
 @njit
 def findSpan(v, knots):

From fa52065bea92dffc39bd9c720928d4ea40feab36 Mon Sep 17 00:00:00 2001
From: adam-urbanczyk <13981538+adam-urbanczyk@users.noreply.github.com>
Date: Fri, 15 Aug 2025 16:47:14 +0200
Subject: [PATCH 55/55] fix test

---
 tests/test_nurbs.py | 4 +---
 1 file changed, 1 insertion(+), 3 deletions(-)

diff --git a/tests/test_nurbs.py b/tests/test_nurbs.py
index 19931463a..9f57a8a46 100644
--- a/tests/test_nurbs.py
+++ b/tests/test_nurbs.py
@@ -87,9 +87,7 @@ def test_dm_2d():
     vparams = np.linspace(0, 1, 100)
     vorder = 2
 
-    params = np.column_stack((uparams, vparams))
-
-    res = designMatrix2D(params, uorder, vorder, uknots, vknots)
+    res = designMatrix2D(uparams, vparams, uorder, vorder, uknots, vknots)
 
     C = res.coo()