From 371acfbebe6327f7fad81832c246a044435b9fff Mon Sep 17 00:00:00 2001 From: adam-urbanczyk <13981538+adam-urbanczyk@users.noreply.github.com> Date: Fri, 20 Mar 2026 21:01:01 +0100 Subject: [PATCH 01/22] Implement direct visualization of Curve/Surface --- cadquery/vis.py | 18 ++++++++++++++---- 1 file changed, 14 insertions(+), 4 deletions(-) diff --git a/cadquery/vis.py b/cadquery/vis.py index f46252964..77c731856 100644 --- a/cadquery/vis.py +++ b/cadquery/vis.py @@ -11,8 +11,9 @@ Edge, ) from .occ_impl.assembly import _loc2vtk, toVTKAssy +from .occ_impl.nurbs import Curve, Surface -from typing import Union, Any, List, Tuple, Iterable, cast, Optional +from typing import Union, List, Tuple, Iterable, cast, Optional from OCP.TopoDS import TopoDS_Shape from OCP.Geom import Geom_BSplineSurface @@ -30,7 +31,6 @@ vtkRenderWindow, vtkWindowToImageFilter, vtkRenderer, - vtkPropCollection, ) from vtkmodules.vtkCommonCore import vtkPoints from vtkmodules.vtkCommonDataModel import vtkCellArray, vtkPolyData @@ -102,6 +102,10 @@ def _split_showables( for el in objs: if instance_of(el, ShapeLike): rv_s.append(el) + if isinstance(el, Curve): + rv_s.append(el.edge()) + if isinstance(el, Surface): + rv_s.append(el.face()) elif isinstance(el, Vector): rv_v.append(el) elif isinstance(el, Location): @@ -195,14 +199,20 @@ def _to_vtk_shapes( def ctrlPts( - s: Union[Face, Edge], + s: Union[Face, Edge, Surface, Curve], size: float = DEFAULT_CTRL_PT_SIZE, color: str = DEFAULT_CTRL_PT_COLOR, ) -> vtkActor: """ - Convert Edge or Face to a vtkActor representing control points. + Convert Edge, Face, Surface or Curve to a vtkActor representing control points. """ + # handle Surface, Curve first + if isinstance(s, Surface): + return ctrlPts(s.face(), size, color) + elif isinstance(s, Curve): + return ctrlPts(s.edge(), size, color) + rv = vtkActor() mapper = vtkPolyDataMapper() From d7e43c20b8d4aeb4165631e2d861dc388996a6b3 Mon Sep 17 00:00:00 2001 From: adam-urbanczyk <13981538+adam-urbanczyk@users.noreply.github.com> Date: Fri, 20 Mar 2026 21:01:22 +0100 Subject: [PATCH 02/22] Fig visibility fixes --- cadquery/fig.py | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/cadquery/fig.py b/cadquery/fig.py index 7a1037851..d5b67f9b5 100644 --- a/cadquery/fig.py +++ b/cadquery/fig.py @@ -472,6 +472,12 @@ def onVisibility(self, event: dict): for act in actors: act.SetVisibility(event["visible"]) + # synchronize state by hand + for act in self.state.actors: + if act["id"] == event["id"]: + act["visible"] = event["visible"] + + self._update_state("actors") self.view.update() def onSelection(self, event: list[str]): From 78656a607b64a586810f5ed75610d7046821514a Mon Sep 17 00:00:00 2001 From: adam-urbanczyk <13981538+adam-urbanczyk@users.noreply.github.com> Date: Fri, 20 Mar 2026 21:02:23 +0100 Subject: [PATCH 03/22] Fix matrix construction and start working on isolines --- cadquery/occ_impl/nurbs.py | 183 ++++++++++++++++++++++++++++++++++++- 1 file changed, 179 insertions(+), 4 deletions(-) diff --git a/cadquery/occ_impl/nurbs.py b/cadquery/occ_impl/nurbs.py index da520f81b..0dcbc933a 100644 --- a/cadquery/occ_impl/nurbs.py +++ b/cadquery/occ_impl/nurbs.py @@ -4,7 +4,7 @@ from numba import njit as _njit -from typing import NamedTuple, Optional, Tuple, List, Union, cast +from typing import NamedTuple, Optional, Tuple, List, Union, cast, Literal from math import comb @@ -95,10 +95,15 @@ class COO(NamedTuple): i: ArrayI j: ArrayI v: Array + shape: tuple[int, int] def coo(self): - return sp.coo_matrix((self.v, (self.i, self.j))) + return sp.coo_matrix((self.v, (self.i, self.j)), shape=self.shape) + + def cooarr(self): + + return sp.coo_array((self.v, (self.i, self.j)), shape=self.shape) def csc(self): @@ -289,6 +294,29 @@ def normal(self, u: Array, v: Array) -> Tuple[Array, Array]: return rv, ders[:, 0, 0, :].squeeze() + def isoline(self, param: float, dir: Literal["u", "v"] = "u") -> Curve: + """ + Extract an isoline. + """ + + if dir == "u": + knots = self.vknots + knots_orig = self.uknots + periodic = self.uperiodic + order = self.uorder + ix = 0 + else: + knots = self.uknots + knots_orig = self.vknots + periodic = self.vperiodic + order = self.vorder + ix = 1 + + C = designMatrix(np.atleast_1d(param), order, knots_orig, periodic).cooarr() + pts = C.tensordot(self.pts, axes=((1,), (ix,))).squeeze() + + return Curve(pts, knots, order, periodic) + # %% basis functions @@ -906,6 +934,7 @@ def designMatrix(u: Array, order: int, knots: Array, periodic: bool = False) -> i=np.empty(n * nu, dtype=np.int64), j=np.empty(n * nu, dtype=np.int64), v=np.empty(n * nu), + shape=(nu, len(knots) - 1 + (not periodic) * order), ) # loop over param values @@ -972,6 +1001,11 @@ def designMatrix2D( i=np.empty(ni * nj, dtype=np.int64), j=np.empty(ni * nj, dtype=np.int64), v=np.empty(ni * nj), + shape=( + ni, + (len(uknots) - 1 + (not uperiodic) * uorder) + * (len(vknots) - 1 + (not vperiodic) * vorder), + ), ) # loop over param values @@ -1030,6 +1064,7 @@ def derMatrix(u: Array, order: int, dorder: int, knots: Array) -> list[COO]: i=np.empty(n * nu, dtype=np.int64), j=np.empty(n * nu, dtype=np.int64), v=np.empty(n * nu), + shape=(nu, len(knots) - 1 + order), ) ) @@ -1084,6 +1119,7 @@ def periodicDerMatrix(u: Array, order: int, dorder: int, knots: Array) -> list[C i=np.empty(n * nu, dtype=np.int64), j=np.empty(n * nu, dtype=np.int64), v=np.empty(n * nu), + shape=(nu, len(knots) - 1), ) ) @@ -1127,6 +1163,7 @@ def periodicDiscretePenalty(us: Array, order: int) -> COO: i=np.empty(nb * ne, dtype=np.int64), j=np.empty(nb * ne, dtype=np.int64), v=np.empty(nb * ne), + shape=(nb, nb), ) if order == 1: @@ -1157,7 +1194,7 @@ def periodicDiscretePenalty(us: Array, order: int) -> COO: @njit -def discretePenalty(us: Array, order: int, splineorder: int = 3) -> COO: +def discretePenalty(us: Array, order: int) -> COO: if order not in (1, 2): raise ValueError( @@ -1175,6 +1212,7 @@ def discretePenalty(us: Array, order: int, splineorder: int = 3) -> COO: i=np.empty(nb * ne, dtype=np.int64), j=np.empty(nb * ne, dtype=np.int64), v=np.empty(nb * ne), + shape=(nb, nb), ) if order == 1: @@ -1294,6 +1332,11 @@ def penaltyMatrix2D( i=np.empty(ni * nj, dtype=np.int64), j=np.empty(ni * nj, dtype=np.int64), v=np.empty(ni * nj), + shape=( + ni, + (len(uknots) - 1 + (not uperiodic) * uorder) + * (len(vknots) - 1 + (not vperiodic) * vorder), + ), ) ) @@ -1352,6 +1395,60 @@ def uniformGrid( return up, vp +@njit +def uIsoMatrix(surf: Surface, u: float) -> COO: + """ + Create a matrix that extracts the requested u isoline control points from flattened surf control points. + """ + + nu = surf.pts.shape[0] + nv = surf.pts.shape[1] + + block = designMatrix(np.atleast_1d(u), surf.uorder, surf.uknots, surf.uperiodic) + + shape = (1, nu * nv) + + n = len(block.i) + N = nv * n # total number of elements + i = np.empty(N, dtype=np.int64) + j = np.empty(N, dtype=np.int64) + v = np.empty(N) + + for ix in range(nv): + i[ix * n : (ix + 1) * n] = 0 + j[ix * n : (ix + 1) * n] = block.j * nv + v[ix * n : (ix + 1) * n] = block.v + + return COO(i, j, v, shape) + + +@njit +def vIsoMatrix(surf: Surface, v: float) -> COO: + """ + Create a matrix that extracts the requested v isoline control points from flattened surf control points. + """ + + nu = surf.pts.shape[0] + nv = surf.pts.shape[1] + + block = designMatrix(np.atleast_1d(v), surf.vorder, surf.vknots, surf.vperiodic) + + shape = (1, nu * nv) + + n = len(block.i) + N = nu * n # total number of elements + i = np.empty(N, dtype=np.int64) + j = np.empty(N, dtype=np.int64) + vals = np.empty(N) + + for ix in range(nu): + i[ix * n : (ix + 1) * n] = 0 + j[ix * n : (ix + 1) * n] = block.j + ix * nv + vals[ix * n : (ix + 1) * n] = block.v + + return COO(i, j, vals, shape) + + # %% construction @@ -1661,7 +1758,7 @@ def approximate2D( lam: float = 0, ) -> Surface: """ - Simple 2D surface approximation (without any penalty). + Simple 2D surface approximation. """ # process the knots @@ -1711,6 +1808,84 @@ def approximate2D( return rv +def constrainedApproximate2D( + As: tuple[Optional[COO], Optional[COO], Optional[COO]], + bs: tuple[Optional[Array], Optional[Array], Optional[Array]], + data: Array, + u: Array, + v: Array, + uorder: int, + vorder: int, + uknots: int | Array = 50, + vknots: int | Array = 50, + uperiodic: bool = False, + vperiodic: bool = False, + penalty: int = 3, + lam: float = 0, +) -> Surface: + """ + 2D surface approximation with additional optional constraints per (x,y,z) component. + """ + + # 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 design matrix + C = designMatrix2D( + u, v, uorder, vorder, uknots_, vknots_, uperiodic, vperiodic + ).csc() + + # handle penalties if requested + if lam: + # construct the penalty grid + up, vp = uniformGrid(uknots_, vknots_, uorder, vorder, uperiodic, vperiodic) + + # construct the derivative matrices + penalties = penaltyMatrix2D( + up, vp, 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 length of the parametric domain + Lv = vknots_[-1] - vknots_[0] # u length 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 the augmented normal equations + pts = [] + for i, (A, b) in enumerate(zip(As, bs)): + if A and b: + LHS = sp.bmat(((CtC, A.coo().T,), (A.coo(), None))) + RHS = np.stack((C.T @ data[:, i], b)) + else: + LHS = CtC + RHS = C.T @ data[:, i] + + D, L, P = ldl(LHS, False) + sol = ldl_solve(RHS, D, L, P).toarray() + pts.append(sol[: CtC.shape[0]]) + + # construct the result + rv = Surface( + np.stack(pts, axis=1).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 c94da159466e3dd277a4f30a478551efc63cf6e3 Mon Sep 17 00:00:00 2001 From: adam-urbanczyk <13981538+adam-urbanczyk@users.noreply.github.com> Date: Fri, 20 Mar 2026 23:40:49 +0100 Subject: [PATCH 04/22] Fix isoline matrix creation --- cadquery/occ_impl/nurbs.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/cadquery/occ_impl/nurbs.py b/cadquery/occ_impl/nurbs.py index 0dcbc933a..ba2f4797b 100644 --- a/cadquery/occ_impl/nurbs.py +++ b/cadquery/occ_impl/nurbs.py @@ -1406,7 +1406,7 @@ def uIsoMatrix(surf: Surface, u: float) -> COO: block = designMatrix(np.atleast_1d(u), surf.uorder, surf.uknots, surf.uperiodic) - shape = (1, nu * nv) + shape = (nv, nu * nv) n = len(block.i) N = nv * n # total number of elements @@ -1415,8 +1415,8 @@ def uIsoMatrix(surf: Surface, u: float) -> COO: v = np.empty(N) for ix in range(nv): - i[ix * n : (ix + 1) * n] = 0 - j[ix * n : (ix + 1) * n] = block.j * nv + i[ix * n : (ix + 1) * n] = ix + j[ix * n : (ix + 1) * n] = block.j * nv + ix v[ix * n : (ix + 1) * n] = block.v return COO(i, j, v, shape) @@ -1433,7 +1433,7 @@ def vIsoMatrix(surf: Surface, v: float) -> COO: block = designMatrix(np.atleast_1d(v), surf.vorder, surf.vknots, surf.vperiodic) - shape = (1, nu * nv) + shape = (nu, nu * nv) n = len(block.i) N = nu * n # total number of elements @@ -1442,7 +1442,7 @@ def vIsoMatrix(surf: Surface, v: float) -> COO: vals = np.empty(N) for ix in range(nu): - i[ix * n : (ix + 1) * n] = 0 + i[ix * n : (ix + 1) * n] = ix j[ix * n : (ix + 1) * n] = block.j + ix * nv vals[ix * n : (ix + 1) * n] = block.v From 5195db10830dcdfd19cb19a0e20e65f9cfe2b15a Mon Sep 17 00:00:00 2001 From: adam-urbanczyk <13981538+adam-urbanczyk@users.noreply.github.com> Date: Mon, 23 Mar 2026 08:42:30 +0100 Subject: [PATCH 05/22] Fix wrong dimensions --- cadquery/occ_impl/nurbs.py | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/cadquery/occ_impl/nurbs.py b/cadquery/occ_impl/nurbs.py index ba2f4797b..46590bd91 100644 --- a/cadquery/occ_impl/nurbs.py +++ b/cadquery/occ_impl/nurbs.py @@ -934,7 +934,7 @@ def designMatrix(u: Array, order: int, knots: Array, periodic: bool = False) -> i=np.empty(n * nu, dtype=np.int64), j=np.empty(n * nu, dtype=np.int64), v=np.empty(n * nu), - shape=(nu, len(knots) - 1 + (not periodic) * order), + shape=(nu, len(knots) - 1 - (not periodic) * order), ) # loop over param values @@ -1003,8 +1003,8 @@ def designMatrix2D( v=np.empty(ni * nj), shape=( ni, - (len(uknots) - 1 + (not uperiodic) * uorder) - * (len(vknots) - 1 + (not vperiodic) * vorder), + (len(uknots) - 1 - (not uperiodic) * uorder) + * (len(vknots) - 1 - (not vperiodic) * vorder), ), ) @@ -1064,7 +1064,7 @@ def derMatrix(u: Array, order: int, dorder: int, knots: Array) -> list[COO]: i=np.empty(n * nu, dtype=np.int64), j=np.empty(n * nu, dtype=np.int64), v=np.empty(n * nu), - shape=(nu, len(knots) - 1 + order), + shape=(nu, len(knots) - 1 - order), ) ) @@ -1334,8 +1334,8 @@ def penaltyMatrix2D( v=np.empty(ni * nj), shape=( ni, - (len(uknots) - 1 + (not uperiodic) * uorder) - * (len(vknots) - 1 + (not vperiodic) * vorder), + (len(uknots) - 1 - (not uperiodic) * uorder) + * (len(vknots) - 1 - (not vperiodic) * vorder), ), ) ) @@ -1616,7 +1616,7 @@ def approximate( # discrete + exact derivatives if penalty > order: Pexact = derMatrix(up, order, order - 1, knots_)[-1].csc() - Pdiscrete = discretePenalty(up, penalty - order, order).csc() + Pdiscrete = discretePenalty(up, penalty - order).csc() P = Pdiscrete @ Pexact @@ -1695,7 +1695,7 @@ def approximate( # discrete + exact derivatives if penalty > order: Pexact = derMatrix(up, order, order - 1, knots_)[-1].csc() - Pdiscrete = discretePenalty(up, penalty - order, order).csc() + Pdiscrete = discretePenalty(up, penalty - order).csc() P = Pdiscrete @ Pexact From a50683d619ae73d970e5dff1cab0512cb1c35984 Mon Sep 17 00:00:00 2001 From: adam-urbanczyk <13981538+adam-urbanczyk@users.noreply.github.com> Date: Mon, 23 Mar 2026 18:40:49 +0100 Subject: [PATCH 06/22] Add some docstrings + spelling --- cadquery/occ_impl/nurbs.py | 62 ++++++++++++++++++++++++++++---------- 1 file changed, 46 insertions(+), 16 deletions(-) diff --git a/cadquery/occ_impl/nurbs.py b/cadquery/occ_impl/nurbs.py index 46590bd91..17402cd41 100644 --- a/cadquery/occ_impl/nurbs.py +++ b/cadquery/occ_impl/nurbs.py @@ -185,6 +185,9 @@ class Surface(NamedTuple): vperiodic: bool def surface(self) -> Geom_BSplineSurface: + """ + Convert to a OCCT Surface object. + """ unique_knots, mults_arr = np.unique(self.uknots, return_counts=True) uknots = _colRealArray(unique_knots) @@ -207,6 +210,9 @@ def surface(self) -> Geom_BSplineSurface: ) def face(self, tol: float = 1e-3) -> Face: + """ + Convert to a Face object. + """ return Face(BRepBuilderAPI_MakeFace(self.surface(), tol).Shape()) @@ -856,7 +862,7 @@ def nbSurfaceDer( v, vorder, vknots, vperiodic ) - # number of param values + # number of parameter values nu = np.size(u) # chunk sizes @@ -917,7 +923,7 @@ def designMatrix(u: Array, order: int, knots: Array, periodic: bool = False) -> # extend the knots u_, knots_ext, minspan, maxspan, deltaspan = _preprocess(u, order, knots, periodic) - # number of param values + # number of parameter values nu = len(u) # number of basis functions @@ -937,7 +943,7 @@ def designMatrix(u: Array, order: int, knots: Array, periodic: bool = False) -> shape=(nu, len(knots) - 1 - (not periodic) * order), ) - # loop over param values + # loop over parameter values for i in range(nu): ui = u_[i] @@ -972,7 +978,7 @@ def designMatrix2D( Create a sparse tensor product design matrix. """ - # extend the knots and preprocess + # extend the knots and pre-process u_, uknots_ext, minspanu, maxspanu, deltaspanu = _preprocess( u, uorder, uknots, uperiodic ) @@ -980,7 +986,7 @@ def designMatrix2D( v, vorder, vknots, vperiodic ) - # number of param values + # number of parameter values ni = len(u) # chunk size @@ -1008,7 +1014,7 @@ def designMatrix2D( ), ) - # loop over param values + # loop over parameter values for i in range(ni): ui, vi = u_[i], v_[i] @@ -1046,7 +1052,7 @@ 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 + # number of parameter values nu = np.size(u) # chunk size @@ -1068,7 +1074,7 @@ def derMatrix(u: Array, order: int, dorder: int, knots: Array) -> list[COO]: ) ) - # loop over param values + # loop over parameter values for i in range(nu): ui = u[i] @@ -1098,7 +1104,7 @@ def periodicDerMatrix(u: Array, order: int, dorder: int, knots: Array) -> list[C (knots[-order:-1] - knots[-1], knots, knots[-1] + knots[1:order]) ) - # number of param values + # number of parameter values nu = len(u) # number of basis functions @@ -1123,7 +1129,7 @@ def periodicDerMatrix(u: Array, order: int, dorder: int, knots: Array) -> list[C ) ) - # loop over param values + # loop over parameter values for i in range(nu): ui = u[i] @@ -1138,7 +1144,7 @@ def periodicDerMatrix(u: Array, order: int, dorder: int, knots: Array) -> list[C 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 + ) % nb # NB: this is due to periodicity rv[di].v[i * n : (i + 1) * n] = temp[di, :] return rv @@ -1146,6 +1152,9 @@ def periodicDerMatrix(u: Array, order: int, dorder: int, knots: Array) -> list[C @njit def periodicDiscretePenalty(us: Array, order: int) -> COO: + """ + Create a periodic discrete penalty matrix for approximating higher order penalties. + """ if order not in (1, 2): raise ValueError( @@ -1195,6 +1204,9 @@ def periodicDiscretePenalty(us: Array, order: int) -> COO: @njit def discretePenalty(us: Array, order: int) -> COO: + """ + Create a discrete penalty matrix for approximating higher order penalties. + """ if order not in (1, 2): raise ValueError( @@ -1300,7 +1312,7 @@ def penaltyMatrix2D( Create sparse tensor product 2D derivative matrices. """ - # extend the knots and preprocess + # extend the knots and pre-process u_, uknots_ext, minspanu, maxspanu, deltaspanu = _preprocess( u, uorder, uknots, uperiodic ) @@ -1308,7 +1320,7 @@ def penaltyMatrix2D( v, vorder, vknots, vperiodic ) - # number of param values + # number of parameter values ni = len(u) # chunk size @@ -1324,7 +1336,7 @@ def penaltyMatrix2D( utemp = np.zeros((dorder + 1, nu)) vtemp = np.zeros((dorder + 1, nv)) - # initialize the emptry matrices + # initialize the empty matrices rv = [] for i in range(dorder + 1): rv.append( @@ -1340,7 +1352,7 @@ def penaltyMatrix2D( ) ) - # loop over param values + # loop over parameter values for i in range(ni): ui, vi = u_[i], v_[i] @@ -1472,6 +1484,9 @@ def periodicApproximate( penalty: int = 4, lam: float = 0, ) -> Curve: + """ + Approximate a periodic curve. + """ npts = data.shape[0] @@ -1529,6 +1544,9 @@ def periodicApproximate( penalty: int = 4, lam: float = 0, ) -> List[Curve]: + """ + Approximate multiple periodic curves using same parametrization. + """ rv = [] @@ -1589,6 +1607,9 @@ def approximate( lam: float = 0, tangents: Optional[Tuple[Array, Array]] = None, ) -> Curve: + """ + Approximate a curve. + """ npts = data.shape[0] @@ -1666,6 +1687,9 @@ def approximate( lam: float = 0, tangents: Optional[Union[Tuple[Array, Array], List[Tuple[Array, Array]]]] = None, ) -> List[Curve]: + """ + Approximate multiple curves using the same parametrization. + """ rv = [] @@ -1758,7 +1782,7 @@ def approximate2D( lam: float = 0, ) -> Surface: """ - Simple 2D surface approximation. + Approximate a 2D surface. """ # process the knots @@ -1887,6 +1911,9 @@ def constrainedApproximate2D( def periodicLoft(*curves: Curve, order: int = 3) -> Surface: + """ + Periodic loft from curves. + """ nknots: int = len(curves) + 1 @@ -1917,6 +1944,9 @@ def loft( penalty: int = 4, tangents: Optional[List[Tuple[Array, Array]]] = None, ) -> Surface: + """ + Regular (non-periodic) loft form curves. + """ nknots: int = len(curves) From 44a54a6891f03f9e2105fc5d9c234e75ac9c5bd1 Mon Sep 17 00:00:00 2001 From: adam-urbanczyk <13981538+adam-urbanczyk@users.noreply.github.com> Date: Fri, 27 Mar 2026 14:51:29 +0100 Subject: [PATCH 07/22] Preliminary solution for column shift --- cadquery/occ_impl/nurbs.py | 111 +++++++++++++++++++++++++++++++++---- 1 file changed, 101 insertions(+), 10 deletions(-) diff --git a/cadquery/occ_impl/nurbs.py b/cadquery/occ_impl/nurbs.py index 17402cd41..9512a9386 100644 --- a/cadquery/occ_impl/nurbs.py +++ b/cadquery/occ_impl/nurbs.py @@ -4,7 +4,7 @@ from numba import njit as _njit -from typing import NamedTuple, Optional, Tuple, List, Union, cast, Literal +from typing import NamedTuple, Optional, Tuple, List, Union, cast, Literal, Sequence from math import comb @@ -23,8 +23,8 @@ from OCP.BRepBuilderAPI import BRepBuilderAPI_MakeEdge, BRepBuilderAPI_MakeFace from .shapes import Face, Edge - -from multimethod import multidispatch +from .geom import Vector +from ..utils import multidispatch njit = _njit(cache=True, error_model="numpy", fastmath=True, nogil=True, parallel=False) @@ -81,6 +81,16 @@ def _colIntArray(knots: NDArray) -> TColStd_Array1OfInteger: return rv +def vec2array(vecs: Sequence[Vector]) -> NDArray: + + return np.array([v.toTuple() for v in vecs]) + + +def array2vec(data: NDArray) -> list[Vector]: + + return [Vector(*row) for row in data] + + # %% vocabulary types Array = ndarray # NDArray[np.floating] @@ -332,7 +342,7 @@ 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, + Helper for handling periodicity. This function extends the knot vector, wraps the parameters and calculates the delta span. """ @@ -1029,8 +1039,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)) % nu_total) * nv_total - + ((vspan - vorder + np.arange(nv)) % nv_total)[:, np.newaxis] + ((uspan - uorder + 1 + np.arange(nu)) % nu_total) * nv_total + + ((vspan - vorder + 1 + np.arange(nv)) % nv_total)[:, np.newaxis] ).ravel() rv.v[i * nj : (i + 1) * nj] = (utemp * vtemp[:, np.newaxis]).ravel() @@ -1535,7 +1545,7 @@ def periodicApproximate( return rv -@periodicApproximate.register +@multidispatch def periodicApproximate( data: List[Array], us: Optional[Array] = None, @@ -1677,7 +1687,7 @@ def approximate( return rv -@approximate.register +@multidispatch def approximate( data: List[Array], us: Optional[Array] = None, @@ -1832,6 +1842,7 @@ def approximate2D( return rv +@multidispatch def constrainedApproximate2D( As: tuple[Optional[COO], Optional[COO], Optional[COO]], bs: tuple[Optional[Array], Optional[Array], Optional[Array]], @@ -1883,9 +1894,9 @@ def constrainedApproximate2D( # solve the augmented normal equations pts = [] for i, (A, b) in enumerate(zip(As, bs)): - if A and b: + if A is not None and b is not None: LHS = sp.bmat(((CtC, A.coo().T,), (A.coo(), None))) - RHS = np.stack((C.T @ data[:, i], b)) + RHS = np.concatenate((C.T @ data[:, i], b)) else: LHS = CtC RHS = C.T @ data[:, i] @@ -1910,6 +1921,86 @@ def constrainedApproximate2D( return rv +@multidispatch +def constrainedApproximate2D( + A: COO, + b: Array, + data: Array, + u: Array, + v: Array, + uorder: int, + vorder: int, + uknots: int | Array = 50, + vknots: int | Array = 50, + uperiodic: bool = False, + vperiodic: bool = False, + penalty: int = 3, + lam: float = 0, +) -> Surface: + """ + 2D surface approximation with additional constraint for all components (x,y,z) combined. + """ + + # 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 design matrix + C = designMatrix2D( + u, v, uorder, vorder, uknots_, vknots_, uperiodic, vperiodic + ).csc() + + # handle penalties if requested + if lam: + # construct the penalty grid + up, vp = uniformGrid(uknots_, vknots_, uorder, vorder, uperiodic, vperiodic) + + # construct the derivative matrices + penalties = penaltyMatrix2D( + up, vp, 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 length of the parametric domain + Lv = vknots_[-1] - vknots_[0] # u length 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 + + # assemble one large block diagonal matrix + CtC_xyz = sp.block_diag((CtC, CtC, CtC)) + + # solve the augmented normal equations + LHS = sp.bmat(((CtC_xyz, A.coo().T,), (A.coo(), None))) + RHS = np.stack((CtC_xyz.T @ data.flatten(order="F"), b)) + + D, L, P = ldl(LHS, False) + sol = ldl_solve(RHS, D, L, P).toarray() + pts = sol[: CtC_xyz.shape[0]] + + # construct the result + rv = Surface( + pts.reshape((-1, 3), order="F").reshape( + ( + len(uknots_) - int(uperiodic), + len(vknots_) - int(vperiodic), + 3, + ) # FIXME: check if this reshape is correct + ), + uknots_, + vknots_, + uorder, + vorder, + uperiodic, + vperiodic, + ) + + return rv + + def periodicLoft(*curves: Curve, order: int = 3) -> Surface: """ Periodic loft from curves. From e59cbf1ab8d49166219c110674d4fce211bcbad8 Mon Sep 17 00:00:00 2001 From: adam-urbanczyk <13981538+adam-urbanczyk@users.noreply.github.com> Date: Sun, 29 Mar 2026 20:52:21 +0200 Subject: [PATCH 08/22] Fix off by one bug for periodic surfaces and curves --- cadquery/occ_impl/nurbs.py | 19 ++++++++++++------- tests/test_nurbs.py | 25 +++++++++++++++++++++++-- 2 files changed, 35 insertions(+), 9 deletions(-) diff --git a/cadquery/occ_impl/nurbs.py b/cadquery/occ_impl/nurbs.py index 9512a9386..524d214aa 100644 --- a/cadquery/occ_impl/nurbs.py +++ b/cadquery/occ_impl/nurbs.py @@ -339,7 +339,7 @@ def isoline(self, param: float, dir: Literal["u", "v"] = "u") -> Curve: @njiti def _preprocess( - u: Array, order: int, knots: Array, periodic: float + u: Array, order: int, knots: Array, periodic: bool ) -> Tuple[Array, Array, Optional[int], Optional[int], int]: """ Helper for handling periodicity. This function extends the knot vector, @@ -966,7 +966,7 @@ def designMatrix(u: Array, order: int, knots: Array, periodic: bool = False) -> # update the matrix rv.i[i * n : (i + 1) * n] = i rv.j[i * n : (i + 1) * n] = ( - span - order + np.arange(n) + span - order + int(periodic) + np.arange(n) ) % nb # NB: this is due to periodicity rv.v[i * n : (i + 1) * n] = temp @@ -1039,8 +1039,10 @@ def designMatrix2D( # update the matrix rv.i[i * nj : (i + 1) * nj] = i rv.j[i * nj : (i + 1) * nj] = ( - ((uspan - uorder + 1 + np.arange(nu)) % nu_total) * nv_total - + ((vspan - vorder + 1 + np.arange(nv)) % nv_total)[:, np.newaxis] + ((uspan - uorder + int(uperiodic) + np.arange(nu)) % nu_total) * nv_total + + ((vspan - vorder + int(vperiodic) + np.arange(nv)) % nv_total)[ + :, np.newaxis + ] ).ravel() rv.v[i * nj : (i + 1) * nj] = (utemp * vtemp[:, np.newaxis]).ravel() @@ -1153,7 +1155,7 @@ def periodicDerMatrix(u: Array, order: int, dorder: int, knots: Array) -> list[C 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) + span - order + 1 + np.arange(n) ) % nb # NB: this is due to periodicity rv[di].v[i * n : (i + 1) * n] = temp[di, :] @@ -1381,8 +1383,11 @@ def penaltyMatrix2D( 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] + ((uspan - uorder + int(uperiodic) + np.arange(nu)) % nu_total) + * nv_total + + ((vspan - vorder + int(vperiodic) + np.arange(nv)) % nv_total)[ + :, np.newaxis + ] ).ravel() rv[dv].v[i * nj : (i + 1) * nj] = ( utemp[du, :] * vtemp[dv, :, np.newaxis] diff --git a/tests/test_nurbs.py b/tests/test_nurbs.py index 5270a18ce..921335644 100644 --- a/tests/test_nurbs.py +++ b/tests/test_nurbs.py @@ -12,6 +12,7 @@ periodicApproximate, periodicLoft, loft, + array2vec, ) from cadquery.func import circle, torus @@ -277,7 +278,8 @@ def test_approximate(): def test_periodic_approximate(): - pts_ = circle(1).sample(100)[0] + circ = circle(1) + pts_ = circ.sample(100)[0] pts = np.array([list(p) for p in pts_]) crv = periodicApproximate(pts) @@ -286,6 +288,12 @@ def test_periodic_approximate(): assert e.isValid() assert e.Length() == approx(2 * np.pi) + # check params + us0 = circ.params(pts_) + us1 = e.params(pts_) + + assert np.allclose(us0, np.array(us1) * 2 * np.pi) + # multiple approximate crvs = periodicApproximate([pts, pts]) @@ -352,5 +360,18 @@ def test_approximate2D(lam, penalty): f = surf.face() + # general sanity checks assert f.isValid() - assert f.Area() == approx(t.Area(), rel=1e-4) + assert f.Area() == approx(t.Area(), rel=1e-3) + + # check the stability of the parameters + us0 = us[None, :].repeat(len(vs), 0).ravel() + vs0 = vs[:, None].repeat(len(us), 1).ravel() + + us2, vs2 = f.params(array2vec(pts)) + + delta_u = np.array(us2 - us0) + delta_v = np.array(vs2 - vs0) + + assert np.allclose(np.where(delta_u >= 1, delta_u, 0) % 1, 0) + assert np.allclose(np.where(delta_v >= 1, delta_v, 0) % 1, 0) From 32af7015e69f37e4d7f7a8aa96d87b2fe9c52a2e Mon Sep 17 00:00:00 2001 From: adam-urbanczyk <13981538+adam-urbanczyk@users.noreply.github.com> Date: Mon, 30 Mar 2026 21:26:16 +0200 Subject: [PATCH 09/22] Fix periodic surface and curve evaluation --- cadquery/occ_impl/nurbs.py | 16 ++++++--- tests/test_nurbs.py | 74 +++++++++++++++++++++++++++++++++++++- 2 files changed, 84 insertions(+), 6 deletions(-) diff --git a/cadquery/occ_impl/nurbs.py b/cadquery/occ_impl/nurbs.py index 524d214aa..35ffbf064 100644 --- a/cadquery/occ_impl/nurbs.py +++ b/cadquery/occ_impl/nurbs.py @@ -641,7 +641,7 @@ def nbCurve( # multiply by ctrl points for j in range(order + 1): - out[i, :] += temp[j] * pts[(span - order + j) % nb, :] + out[i, :] += temp[j] * pts[(span - order + int(periodic) + j) % nb, :] return out @@ -705,7 +705,9 @@ def nbCurveDer( # 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, :] + out[i, k, :] += ( + temp[k, j] * pts[(span - order + int(periodic) + j) % nb, :] + ) return out @@ -791,14 +793,14 @@ def nbSurface( nbBasis(uspan, ui, uorder, uknots_ext, utemp) nbBasis(vspan, vi, vorder, vknots_ext, vtemp) - uind = uspan - uorder + uind = uspan - uorder + int(uperiodic) temp = np.empty(3) # multiply by ctrl points: Nu.T*P*Nv for j in range(vorder + 1): temp[:] = 0.0 - vind = vspan - vorder + j + vind = vspan - vorder + int(vperiodic) + j # calculate Nu.T*P for k in range(uorder + 1): @@ -907,7 +909,11 @@ def nbSurfaceDer( for r in range(uorder + 1): temp[s, :] += ( utemp[k, r] - * pts[(uspan - uorder + r) % nub, (vspan - vorder + s) % nvb, :] + * pts[ + (uspan - uorder + int(uperiodic) + r) % nub, + (vspan - vorder + int(vperiodic) + s) % nvb, + :, + ] ) # remaining derivative orders: dk + du <= dorder diff --git a/tests/test_nurbs.py b/tests/test_nurbs.py index 921335644..ca0092a42 100644 --- a/tests/test_nurbs.py +++ b/tests/test_nurbs.py @@ -13,9 +13,10 @@ periodicLoft, loft, array2vec, + vec2array, ) -from cadquery.func import circle, torus +from cadquery.func import circle, torus, ellipse, spline, plane from cadquery import Vector from OCP.gp import gp_Pnt, gp_Vec @@ -375,3 +376,74 @@ def test_approximate2D(lam, penalty): assert np.allclose(np.where(delta_u >= 1, delta_u, 0) % 1, 0) assert np.allclose(np.where(delta_v >= 1, delta_v, 0) % 1, 0) + + +EDGES = ( + periodicApproximate(vec2array(ellipse(2, 1).sample(100)[0])).edge(), + spline([(0, 0), (1, 0)], tgts=((0, 1), (1, 0))), +) + +PARAMS = np.array((0, 0.1, 0.5)) + + +@mark.parametrize("e", EDGES) +def test_curve_position(e): + + crv = Curve.fromEdge(e) + + for u in PARAMS: + assert np.allclose(np.array(e.positionAt(u, mode="param").toTuple()), crv(u)) + + +@mark.parametrize("e", EDGES) +def test_curve_tangents(e): + + crv = Curve.fromEdge(e) + + for u in PARAMS: + tgt = crv.der(u, 1)[0, 1, :] + tgt /= np.linalg.norm(tgt) + + assert np.allclose(np.array(e.tangentAt(u, mode="param").toTuple()), tgt) + + +def torus_surf(): + + t = torus(5, 1).face() + + # double periodic surface + us = np.linspace(0, 1, endpoint=False) + vs = np.linspace(0, 1, endpoint=False) + + pts = np.array( + [t.positionAt(u * 2 * np.pi, v * 2 * np.pi).toTuple() for v in vs for u in us] + ) + + surf = approximate2D( + pts, + us[None, :].repeat(len(vs), 0).ravel(), + vs[:, None].repeat(len(us), 1).ravel(), + 3, + 3, + 50, + 50, + uperiodic=True, + vperiodic=True, + penalty=3, + lam=1e-6, + ) + + return surf.face() + + +FACES = (torus_surf(), plane(1, 1).toNURBS()) + + +@mark.parametrize("f", FACES) +def test_surface_positions(f): + + surf = Surface.fromFace(f) + + for u in PARAMS: + for v in PARAMS: + assert np.allclose(np.array(f.positionAt(u, v).toTuple()), surf(u, v)) From b75a130c58ca767507d0598b912db97e4bab9f3f Mon Sep 17 00:00:00 2001 From: adam-urbanczyk <13981538+adam-urbanczyk@users.noreply.github.com> Date: Tue, 31 Mar 2026 08:58:28 +0200 Subject: [PATCH 10/22] Add test for tangents --- cadquery/occ_impl/shapes.py | 16 ++++++++++++++++ tests/test_nurbs.py | 17 +++++++++++++++++ 2 files changed, 33 insertions(+) diff --git a/cadquery/occ_impl/shapes.py b/cadquery/occ_impl/shapes.py index 7ee861614..e28cc752d 100644 --- a/cadquery/occ_impl/shapes.py +++ b/cadquery/occ_impl/shapes.py @@ -3392,6 +3392,22 @@ def normals( return rv_n, rv_p + def tangentAt(self, u: Real, v: Real) -> Tuple[Vector, Vector, Vector]: + """ + Computes tangent vectors at the desired location in the u,v parameter space. + + :returns: a vector representing the tangent directions and the position + :param u: the u parametric location to compute at. + :param v: the v parametric location to compute at. + """ + + p = gp_Pnt() + du = gp_Vec() + dv = gp_Vec() + BRepAdaptor_Surface(self.wrapped).D1(u, v, p, du, dv) + + return Vector(du).normalized(), Vector(dv).normalized(), Vector(p) + def Center(self) -> Vector: Properties = GProp_GProps() diff --git a/tests/test_nurbs.py b/tests/test_nurbs.py index ca0092a42..d73d67853 100644 --- a/tests/test_nurbs.py +++ b/tests/test_nurbs.py @@ -447,3 +447,20 @@ def test_surface_positions(f): for u in PARAMS: for v in PARAMS: assert np.allclose(np.array(f.positionAt(u, v).toTuple()), surf(u, v)) + + +@mark.parametrize("f", FACES) +def test_surface_tangents(f): + + surf = Surface.fromFace(f) + + for u in PARAMS: + for v in PARAMS: + dun, dvn, p = f.tangentAt(u, v) + der = surf.der(u, v, 1).squeeze() + du = der[1, 0, :] + dv = der[0, 1, :] + + assert np.allclose(np.array(dun.toTuple()), du / np.linalg.norm(du)) + assert np.allclose(np.array(dvn.toTuple()), dv / np.linalg.norm(dv)) + assert np.allclose(np.array(p.toTuple()), der[0, 0, :]) From 660ac168162cc8aed0934c6199c5d6005043282b Mon Sep 17 00:00:00 2001 From: adam-urbanczyk <13981538+adam-urbanczyk@users.noreply.github.com> Date: Tue, 31 Mar 2026 20:05:12 +0200 Subject: [PATCH 11/22] More tests --- tests/test_nurbs.py | 62 ++++++++++++++++++++++++++++++++++++--------- tests/test_vis.py | 15 ++++++++++- 2 files changed, 64 insertions(+), 13 deletions(-) diff --git a/tests/test_nurbs.py b/tests/test_nurbs.py index d73d67853..ce372bccd 100644 --- a/tests/test_nurbs.py +++ b/tests/test_nurbs.py @@ -407,7 +407,8 @@ def test_curve_tangents(e): assert np.allclose(np.array(e.tangentAt(u, mode="param").toTuple()), tgt) -def torus_surf(): +@fixture +def torus_face(): t = torus(5, 1).face() @@ -436,23 +437,35 @@ def torus_surf(): return surf.face() -FACES = (torus_surf(), plane(1, 1).toNURBS()) +@fixture +def torus_surf(torus_face): + + return Surface.fromFace(torus_face) + + +@fixture +def plane_face(): + + return plane(1, 1).toNURBS() + + +FACES = (torus_face, plane_face) -@mark.parametrize("f", FACES) -def test_surface_positions(f): +@fixture(params=FACES) +def test_surface_positions(request): - surf = Surface.fromFace(f) + surf = Surface.fromFace(request.getfixturevalue(request.param)) for u in PARAMS: for v in PARAMS: - assert np.allclose(np.array(f.positionAt(u, v).toTuple()), surf(u, v)) + assert np.allclose(f.positionAt(u, v).toTuple(), surf(u, v)) -@mark.parametrize("f", FACES) -def test_surface_tangents(f): +@fixture(params=FACES) +def test_surface_tangents(request): - surf = Surface.fromFace(f) + surf = Surface.fromFace(request.getfixturevalue(request.param)) for u in PARAMS: for v in PARAMS: @@ -461,6 +474,31 @@ def test_surface_tangents(f): du = der[1, 0, :] dv = der[0, 1, :] - assert np.allclose(np.array(dun.toTuple()), du / np.linalg.norm(du)) - assert np.allclose(np.array(dvn.toTuple()), dv / np.linalg.norm(dv)) - assert np.allclose(np.array(p.toTuple()), der[0, 0, :]) + assert np.allclose(dun.toTuple(), du / np.linalg.norm(du)) + assert np.allclose(dvn.toTuple(), dv / np.linalg.norm(dv)) + assert np.allclose(p.toTuple(), der[0, 0, :]) + + +@mark.parametrize("isoparam", PARAMS) +@mark.parametrize("u", PARAMS) +def test_isolines(torus_surf, isoparam, u): + + uiso = torus_surf.isoline(isoparam) + viso = torus_surf.isoline(isoparam, "v") + + assert isinstance(uiso, Curve) + assert isinstance(viso, Curve) + + pt_u = uiso(u) + pt_v = viso(u) + + # ref + f = torus_surf.face() + uiso_ref = f.isoline(isoparam, "u") + viso_ref = f.isoline(isoparam, "v") + + pt_u_ref = uiso_ref.positionAt(u, mode="param") + pt_v_ref = viso_ref.positionAt(u, mode="param") + + assert np.allclose(pt_u_ref.toTuple(), pt_u) + assert np.allclose(pt_v_ref.toTuple(), pt_v) diff --git a/tests/test_vis.py b/tests/test_vis.py index 9104dd1a0..b59367d73 100644 --- a/tests/test_vis.py +++ b/tests/test_vis.py @@ -1,6 +1,7 @@ from cadquery import Workplane, Assembly, Sketch, Location, Vector -from cadquery.func import circle, sweep, spline, plane, torus +from cadquery.func import circle, sweep, spline, plane, torus, segment from cadquery.vis import show, show_object, vtkAxesActor, ctrlPts, style +from cadquery.occ_impl.nurbs import Surface, Curve import cadquery.vis as vis @@ -139,6 +140,18 @@ def test_show(wp, assy, sk, patch_vtk): show(vtkAxesActor(), [vtkAnnotatedCubeActor()]) +def test_show_nurbs(patch_vtk): + + # smoke tests for Surface and Curve + surf = Surface.fromFace(plane(1, 1).toNURBS()) + crv = Curve.fromEdge(segment((0, 0), (1, 1)).toNURBS()) + + show(surf) + show(crv) + show(ctrlPts(surf)) + show(ctrlPts(crv)) + + def test_screenshot(wp, patch_vtk): # smoke test for now From c7b0bc081ac2d590d6b59ccf7b655ae5c2fb7662 Mon Sep 17 00:00:00 2001 From: adam-urbanczyk <13981538+adam-urbanczyk@users.noreply.github.com> Date: Wed, 1 Apr 2026 08:43:52 +0200 Subject: [PATCH 12/22] Rework tests --- tests/test_nurbs.py | 16 +++++++++------- 1 file changed, 9 insertions(+), 7 deletions(-) diff --git a/tests/test_nurbs.py b/tests/test_nurbs.py index ce372bccd..863fd558a 100644 --- a/tests/test_nurbs.py +++ b/tests/test_nurbs.py @@ -449,23 +449,25 @@ def plane_face(): return plane(1, 1).toNURBS() -FACES = (torus_face, plane_face) +FACES = ("torus_face", "plane_face") -@fixture(params=FACES) -def test_surface_positions(request): +@mark.parametrize("face", FACES) +def test_surface_positions(face, request): - surf = Surface.fromFace(request.getfixturevalue(request.param)) + f = request.getfixturevalue(face) + surf = Surface.fromFace(f) for u in PARAMS: for v in PARAMS: assert np.allclose(f.positionAt(u, v).toTuple(), surf(u, v)) -@fixture(params=FACES) -def test_surface_tangents(request): +@mark.parametrize("face", FACES) +def test_surface_tangents(face, request): - surf = Surface.fromFace(request.getfixturevalue(request.param)) + f = request.getfixturevalue(face) + surf = Surface.fromFace(f) for u in PARAMS: for v in PARAMS: From 03c5c26b0010ad2c430e472fa7d86f25bdb82ef4 Mon Sep 17 00:00:00 2001 From: adam-urbanczyk <13981538+adam-urbanczyk@users.noreply.github.com> Date: Wed, 1 Apr 2026 19:55:01 +0200 Subject: [PATCH 13/22] More tests --- cadquery/occ_impl/shapes.py | 2 +- tests/test_nurbs.py | 89 ++++++++++++++++++++++++++++++++++--- 2 files changed, 85 insertions(+), 6 deletions(-) diff --git a/cadquery/occ_impl/shapes.py b/cadquery/occ_impl/shapes.py index e28cc752d..d24860011 100644 --- a/cadquery/occ_impl/shapes.py +++ b/cadquery/occ_impl/shapes.py @@ -3396,7 +3396,7 @@ def tangentAt(self, u: Real, v: Real) -> Tuple[Vector, Vector, Vector]: """ Computes tangent vectors at the desired location in the u,v parameter space. - :returns: a vector representing the tangent directions and the position + :returns: vectors representing the tangent directions and the position :param u: the u parametric location to compute at. :param v: the v parametric location to compute at. """ diff --git a/tests/test_nurbs.py b/tests/test_nurbs.py index 863fd558a..b30d2382f 100644 --- a/tests/test_nurbs.py +++ b/tests/test_nurbs.py @@ -13,7 +13,10 @@ periodicLoft, loft, array2vec, + uIsoMatrix, + vIsoMatrix, vec2array, + constrainedApproximate2D, ) from cadquery.func import circle, torus, ellipse, spline, plane @@ -277,23 +280,32 @@ def test_approximate(): assert e.Length() == approx(1) -def test_periodic_approximate(): +@mark.parametrize("penalty", (1, 2, 3, 4, 5)) +@mark.parametrize("lam", (0, 1e-9)) +def test_periodic_approximate(penalty, lam): + + EPS = 1e-6 if lam == 0 else 1e-3 circ = circle(1) pts_ = circ.sample(100)[0] pts = np.array([list(p) for p in pts_]) - crv = periodicApproximate(pts) + crv = periodicApproximate(pts, lam=lam, penalty=penalty) e = crv.edge() assert e.isValid() - assert e.Length() == approx(2 * np.pi) + assert e.Length() == approx(2 * np.pi, rel=EPS) # check params us0 = circ.params(pts_) us1 = e.params(pts_) - assert np.allclose(us0, np.array(us1) * 2 * np.pi) + # NB: I'm skipping the first and last point to not handle wrap around + assert np.allclose(us0[1:-1], np.array(us1)[1:-1] * 2 * np.pi, rtol=EPS) + # special case for wrap around + for ix in (0, -1): + delta = us0[ix] - np.array(us1[ix]) * 2 * np.pi + assert abs(delta) < EPS or abs(delta) - 2 * np.pi < EPS # multiple approximate crvs = periodicApproximate([pts, pts]) @@ -302,7 +314,7 @@ def test_periodic_approximate(): e = crv.edge() assert e.isValid() - assert e.Length() == approx(2 * np.pi) + assert e.Length() == approx(2 * np.pi, rel=EPS) def test_periodic_loft(circles, trimmed_circles): @@ -504,3 +516,70 @@ def test_isolines(torus_surf, isoparam, u): assert np.allclose(pt_u_ref.toTuple(), pt_u) assert np.allclose(pt_v_ref.toTuple(), pt_v) + + +@mark.parametrize("lam", [0, 1e-6]) +@mark.parametrize("penalty", [2, 3, 4, 5]) +def test_constrainedApproximate2D(torus_surf, lam, penalty): + + # sample the surface + us_ = np.linspace(0, 1, endpoint=False) + vs_ = np.linspace(0, 1, endpoint=False) + + pts = np.stack([torus_surf(u, v) for v in vs_ for u in us_]).squeeze() + + us = us_[None, :].repeat(len(vs_), 0).ravel() + vs = vs_[:, None].repeat(len(us_), 1).ravel() + + # add some noise to make the problem non-trivial + pts += np.random.randn(*pts.shape) / 100 + + surf = approximate2D( + pts, + us, + vs, + 3, + 3, + 25, + 25, + uperiodic=True, + vperiodic=True, + penalty=penalty, + lam=lam, + ) + + isou = surf.isoline(0.0, "u") + isov = surf.isoline(0.0, "v") + + # check that the constraints will do something + assert not np.allclose(isou.pts[:, 1], 0) + assert not np.allclose(isov.pts[:, 2], 0) + + # constraints per direction + Au = uIsoMatrix(surf, np.array(0.0)) + Av = vIsoMatrix(surf, np.array(0.0)) + by = np.zeros(Au.shape[0]) + bz = np.zeros(Av.shape[0]) + + surf = constrainedApproximate2D( + (None, Au, Av), + (None, by, bz), + pts, + us, + vs, + 3, + 3, + len(surf.uknots), + len(surf.vknots), + uperiodic=True, + vperiodic=True, + penalty=penalty, + lam=lam, + ) + + isou = surf.isoline(0.0, "u") + isov = surf.isoline(0.0, "v") + + # check the constraints + assert np.allclose(isou.pts[:, 1], 0) + assert np.allclose(isov.pts[:, 2], 0) From 9403a90e24965e4371a38f018b7a676fab3a60ed Mon Sep 17 00:00:00 2001 From: adam-urbanczyk <13981538+adam-urbanczyk@users.noreply.github.com> Date: Fri, 3 Apr 2026 09:32:03 +0200 Subject: [PATCH 14/22] Fix constrainedApproximate --- cadquery/occ_impl/nurbs.py | 68 ++++++++++++++++++-------------------- 1 file changed, 33 insertions(+), 35 deletions(-) diff --git a/cadquery/occ_impl/nurbs.py b/cadquery/occ_impl/nurbs.py index 35ffbf064..d0ef8679d 100644 --- a/cadquery/occ_impl/nurbs.py +++ b/cadquery/occ_impl/nurbs.py @@ -25,6 +25,7 @@ from .shapes import Face, Edge from .geom import Vector from ..utils import multidispatch +from ..types import Real njit = _njit(cache=True, error_model="numpy", fastmath=True, nogil=True, parallel=False) @@ -1855,8 +1856,8 @@ def approximate2D( @multidispatch def constrainedApproximate2D( - As: tuple[Optional[COO], Optional[COO], Optional[COO]], - bs: tuple[Optional[Array], Optional[Array], Optional[Array]], + As: COO, + bs: Array, data: Array, u: Array, v: Array, @@ -1867,10 +1868,10 @@ def constrainedApproximate2D( uperiodic: bool = False, vperiodic: bool = False, penalty: int = 3, - lam: float = 0, + lam: Real = 0, ) -> Surface: """ - 2D surface approximation with additional optional constraints per (x,y,z) component. + 2D surface approximation with additional constraint for all components (x,y,z) combined. """ # process the knots @@ -1902,24 +1903,22 @@ def constrainedApproximate2D( else: CtC = C.T @ C + # assemble one large block diagonal matrix + CtC_xyz = sp.block_diag((CtC, CtC, CtC)) + C_xyz = sp.block_diag((C, C, C)) + # solve the augmented normal equations - pts = [] - for i, (A, b) in enumerate(zip(As, bs)): - if A is not None and b is not None: - LHS = sp.bmat(((CtC, A.coo().T,), (A.coo(), None))) - RHS = np.concatenate((C.T @ data[:, i], b)) - else: - LHS = CtC - RHS = C.T @ data[:, i] + LHS = sp.bmat(((CtC_xyz, As.coo().T,), (As.coo(), None))) + RHS = np.concatenate((C_xyz.T @ data.flatten(order="F"), bs)) - D, L, P = ldl(LHS, False) - sol = ldl_solve(RHS, D, L, P).toarray() - pts.append(sol[: CtC.shape[0]]) + D, L, P = ldl(LHS, False) + sol = ldl_solve(RHS, D, L, P).toarray() + pts = sol[: CtC_xyz.shape[0]] # construct the result rv = Surface( - np.stack(pts, axis=1).reshape( - (len(uknots_) - int(uperiodic), len(vknots_) - int(vperiodic), 3) + pts.reshape((-1, 3), order="F").reshape( + (len(uknots_) - int(uperiodic), len(vknots_) - int(vperiodic), 3,) ), uknots_, vknots_, @@ -1934,8 +1933,8 @@ def constrainedApproximate2D( @multidispatch def constrainedApproximate2D( - A: COO, - b: Array, + As: tuple[Optional[COO], Optional[COO], Optional[COO]], + bs: tuple[Optional[Array], Optional[Array], Optional[Array]], data: Array, u: Array, v: Array, @@ -1946,10 +1945,10 @@ def constrainedApproximate2D( uperiodic: bool = False, vperiodic: bool = False, penalty: int = 3, - lam: float = 0, + lam: Real = 0, ) -> Surface: """ - 2D surface approximation with additional constraint for all components (x,y,z) combined. + 2D surface approximation with additional optional constraints per (x,y,z) component. """ # process the knots @@ -1981,25 +1980,24 @@ def constrainedApproximate2D( else: CtC = C.T @ C - # assemble one large block diagonal matrix - CtC_xyz = sp.block_diag((CtC, CtC, CtC)) - # solve the augmented normal equations - LHS = sp.bmat(((CtC_xyz, A.coo().T,), (A.coo(), None))) - RHS = np.stack((CtC_xyz.T @ data.flatten(order="F"), b)) + pts = [] + for i, (A, b) in enumerate(zip(As, bs)): + if A is not None and b is not None: + LHS = sp.bmat(((CtC, A.coo().T,), (A.coo(), None))) + RHS = np.concatenate((C.T @ data[:, i], b)) + else: + LHS = CtC + RHS = C.T @ data[:, i] - D, L, P = ldl(LHS, False) - sol = ldl_solve(RHS, D, L, P).toarray() - pts = sol[: CtC_xyz.shape[0]] + D, L, P = ldl(LHS, False) + sol = ldl_solve(RHS, D, L, P).toarray() + pts.append(sol[: CtC.shape[0]]) # construct the result rv = Surface( - pts.reshape((-1, 3), order="F").reshape( - ( - len(uknots_) - int(uperiodic), - len(vknots_) - int(vperiodic), - 3, - ) # FIXME: check if this reshape is correct + np.stack(pts, axis=1).reshape( + (len(uknots_) - int(uperiodic), len(vknots_) - int(vperiodic), 3) ), uknots_, vknots_, From 0cfafe22d0c98e9f80890c8078d4a087c7770fcc Mon Sep 17 00:00:00 2001 From: adam-urbanczyk <13981538+adam-urbanczyk@users.noreply.github.com> Date: Fri, 3 Apr 2026 09:32:23 +0200 Subject: [PATCH 15/22] Add test and optimize test vectors --- tests/test_nurbs.py | 48 ++++++++++++++++++++++++++++++++++----------- 1 file changed, 37 insertions(+), 11 deletions(-) diff --git a/tests/test_nurbs.py b/tests/test_nurbs.py index b30d2382f..3fe68a62d 100644 --- a/tests/test_nurbs.py +++ b/tests/test_nurbs.py @@ -17,6 +17,7 @@ vIsoMatrix, vec2array, constrainedApproximate2D, + COO, ) from cadquery.func import circle, torus, ellipse, spline, plane @@ -280,8 +281,9 @@ def test_approximate(): assert e.Length() == approx(1) -@mark.parametrize("penalty", (1, 2, 3, 4, 5)) -@mark.parametrize("lam", (0, 1e-9)) +@mark.parametrize( + ("penalty", "lam"), [(1, 1e-9), (2, 1e-9), (3, 1e-9), (4, 1e-9), (5, 1e-9), (2, 0)] +) def test_periodic_approximate(penalty, lam): EPS = 1e-6 if lam == 0 else 1e-3 @@ -343,8 +345,9 @@ def test_loft(circles, trimmed_circles): assert surf2.face().isValid() -@mark.parametrize("lam", [0, 1e-6]) -@mark.parametrize("penalty", [2, 3, 4, 5]) +@mark.parametrize( + ("penalty", "lam"), [(2, 1e-6), (3, 1e-6), (4, 1e-6), (5, 1e-6), (2, 0)] +) def test_approximate2D(lam, penalty): t = torus(5, 1).face() @@ -363,8 +366,8 @@ def test_approximate2D(lam, penalty): vs[:, None].repeat(len(us), 1).ravel(), 3, 3, - 50, - 50, + 15, + 15, uperiodic=True, vperiodic=True, penalty=penalty, @@ -518,8 +521,8 @@ def test_isolines(torus_surf, isoparam, u): assert np.allclose(pt_v_ref.toTuple(), pt_v) -@mark.parametrize("lam", [0, 1e-6]) -@mark.parametrize("penalty", [2, 3, 4, 5]) +@mark.parametrize("lam", [0.0, 1e-6]) +@mark.parametrize("penalty", [3]) def test_constrainedApproximate2D(torus_surf, lam, penalty): # sample the surface @@ -540,8 +543,8 @@ def test_constrainedApproximate2D(torus_surf, lam, penalty): vs, 3, 3, - 25, - 25, + 15, + 15, uperiodic=True, vperiodic=True, penalty=penalty, @@ -555,7 +558,7 @@ def test_constrainedApproximate2D(torus_surf, lam, penalty): assert not np.allclose(isou.pts[:, 1], 0) assert not np.allclose(isov.pts[:, 2], 0) - # constraints per direction + # # constraints per direction Au = uIsoMatrix(surf, np.array(0.0)) Av = vIsoMatrix(surf, np.array(0.0)) by = np.zeros(Au.shape[0]) @@ -583,3 +586,26 @@ def test_constrainedApproximate2D(torus_surf, lam, penalty): # check the constraints assert np.allclose(isou.pts[:, 1], 0) assert np.allclose(isov.pts[:, 2], 0) + + # all constraints at once + Ay = sp.bmat(((sp.coo_matrix(Au.shape), Au.coo(), sp.coo_matrix(Au.shape)),)) + Az = sp.bmat(((sp.coo_matrix(Av.shape), sp.coo_matrix(Av.shape), Av.coo()),)) + + A = sp.vstack((Ay, Az)) + b = np.concatenate((by, bz)) + + surf = constrainedApproximate2D( + COO(A.row, A.col, A.data, A.shape), + b, + pts, + us, + vs, + 3, + 3, + len(surf.uknots), + len(surf.vknots), + uperiodic=True, + vperiodic=True, + penalty=penalty, + lam=lam, + ) From 969afb6a08dd2cd6de62c9d15005028f2bcf516b Mon Sep 17 00:00:00 2001 From: adam-urbanczyk <13981538+adam-urbanczyk@users.noreply.github.com> Date: Fri, 3 Apr 2026 18:57:12 +0200 Subject: [PATCH 16/22] Optimize tests further --- tests/test_nurbs.py | 20 ++++++++------------ 1 file changed, 8 insertions(+), 12 deletions(-) diff --git a/tests/test_nurbs.py b/tests/test_nurbs.py index 3fe68a62d..ce984481d 100644 --- a/tests/test_nurbs.py +++ b/tests/test_nurbs.py @@ -281,9 +281,7 @@ def test_approximate(): assert e.Length() == approx(1) -@mark.parametrize( - ("penalty", "lam"), [(1, 1e-9), (2, 1e-9), (3, 1e-9), (4, 1e-9), (5, 1e-9), (2, 0)] -) +@mark.parametrize(("penalty", "lam"), [(3, 1e-9), (4, 1e-9), (5, 1e-9), (2, 0)]) def test_periodic_approximate(penalty, lam): EPS = 1e-6 if lam == 0 else 1e-3 @@ -345,9 +343,7 @@ def test_loft(circles, trimmed_circles): assert surf2.face().isValid() -@mark.parametrize( - ("penalty", "lam"), [(2, 1e-6), (3, 1e-6), (4, 1e-6), (5, 1e-6), (2, 0)] -) +@mark.parametrize(("penalty", "lam"), [(3, 1e-6), (4, 1e-6), (5, 1e-6), (2, 0)]) def test_approximate2D(lam, penalty): t = torus(5, 1).face() @@ -398,7 +394,7 @@ def test_approximate2D(lam, penalty): spline([(0, 0), (1, 0)], tgts=((0, 1), (1, 0))), ) -PARAMS = np.array((0, 0.1, 0.5)) +PARAMS = np.array((0, 0.1)) @mark.parametrize("e", EDGES) @@ -441,8 +437,8 @@ def torus_face(): vs[:, None].repeat(len(us), 1).ravel(), 3, 3, - 50, - 50, + 10, + 10, uperiodic=True, vperiodic=True, penalty=3, @@ -543,8 +539,8 @@ def test_constrainedApproximate2D(torus_surf, lam, penalty): vs, 3, 3, - 15, - 15, + 10, + 10, uperiodic=True, vperiodic=True, penalty=penalty, @@ -595,7 +591,7 @@ def test_constrainedApproximate2D(torus_surf, lam, penalty): b = np.concatenate((by, bz)) surf = constrainedApproximate2D( - COO(A.row, A.col, A.data, A.shape), + COO.fromSP(A), b, pts, us, From 8864a6e41a37d0cc4beeced788937b4bd7f926a7 Mon Sep 17 00:00:00 2001 From: adam-urbanczyk <13981538+adam-urbanczyk@users.noreply.github.com> Date: Fri, 3 Apr 2026 18:58:40 +0200 Subject: [PATCH 17/22] Add fromSP helper --- cadquery/occ_impl/nurbs.py | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/cadquery/occ_impl/nurbs.py b/cadquery/occ_impl/nurbs.py index d0ef8679d..035a781f0 100644 --- a/cadquery/occ_impl/nurbs.py +++ b/cadquery/occ_impl/nurbs.py @@ -124,6 +124,16 @@ def csr(self): return self.coo().tocsr() + @classmethod + def fromSP(cls, sparr: sp.sparray | sp.spmatrix): + """ + Helper for converting from any SciPy sparse array or matrix. + """ + + tmp = sparr.tocoo() + + return cls(tmp.row, tmp.col, tmp.data, tmp.shape) + class Curve(NamedTuple): """ From 778c8f2438313c14f6b724850071c1e2f52b75da Mon Sep 17 00:00:00 2001 From: AU Date: Sat, 4 Apr 2026 11:01:01 +0200 Subject: [PATCH 18/22] Pin mypy for now --- environment.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/environment.yml b/environment.yml index 0e1b50492..451362e10 100644 --- a/environment.yml +++ b/environment.yml @@ -9,7 +9,7 @@ dependencies: - pyparsing>=3.0.0 - sphinx=9.1.0 # - sphinx_rtd_theme=3.1.0 - - mypy + - mypy=1.19.* - codecov - pytest - pytest-cov From c3c38ad127ad35bec689e8faff47a33a5c89955d Mon Sep 17 00:00:00 2001 From: AU Date: Sun, 5 Apr 2026 15:19:55 +0200 Subject: [PATCH 19/22] Debugging CI failure --- environment.yml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/environment.yml b/environment.yml index 451362e10..29a72dbd9 100644 --- a/environment.yml +++ b/environment.yml @@ -9,7 +9,7 @@ dependencies: - pyparsing>=3.0.0 - sphinx=9.1.0 # - sphinx_rtd_theme=3.1.0 - - mypy=1.19.* + - mypy#=1.19.* - codecov - pytest - pytest-cov @@ -32,6 +32,6 @@ dependencies: - trame-vuetify - pip - pip: - - sphinx_rtd_theme==3.1.0 +# - sphinx_rtd_theme==3.1.0 - --editable=. - git+https://github.com/cadquery/black.git@cq From 73186ecc35099e91b831e5ce374973417ebbaf39 Mon Sep 17 00:00:00 2001 From: AU Date: Sun, 5 Apr 2026 15:37:02 +0200 Subject: [PATCH 20/22] Typo --- environment.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/environment.yml b/environment.yml index 29a72dbd9..8dd3b605b 100644 --- a/environment.yml +++ b/environment.yml @@ -9,7 +9,7 @@ dependencies: - pyparsing>=3.0.0 - sphinx=9.1.0 # - sphinx_rtd_theme=3.1.0 - - mypy#=1.19.* + - mypy #=1.19.* - codecov - pytest - pytest-cov From fbe61b6c5f72ed631f25b32956a1caf4a8d3cbf8 Mon Sep 17 00:00:00 2001 From: adam-urbanczyk <13981538+adam-urbanczyk@users.noreply.github.com> Date: Sun, 5 Apr 2026 16:22:54 +0200 Subject: [PATCH 21/22] Solve previously undetectd mypy issue --- cadquery/occ_impl/shapes.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cadquery/occ_impl/shapes.py b/cadquery/occ_impl/shapes.py index d24860011..28db281f7 100644 --- a/cadquery/occ_impl/shapes.py +++ b/cadquery/occ_impl/shapes.py @@ -617,7 +617,7 @@ def geomType(self) -> Geoms: if isinstance(tr, str): rv = tr elif tr is BRepAdaptor_Curve: - rv = geom_LUT_EDGE[tr(self.wrapped).GetType()] + rv = geom_LUT_EDGE[tr(tcast(TopoDS_Edge, self.wrapped)).GetType()] else: rv = geom_LUT_FACE[tr(self.wrapped).GetType()] From ca16bfb58af00d53b19c53924b7931db6dc11f31 Mon Sep 17 00:00:00 2001 From: AU Date: Sun, 5 Apr 2026 16:47:56 +0200 Subject: [PATCH 22/22] Remove temp changes --- environment.yml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/environment.yml b/environment.yml index 8dd3b605b..0e1b50492 100644 --- a/environment.yml +++ b/environment.yml @@ -9,7 +9,7 @@ dependencies: - pyparsing>=3.0.0 - sphinx=9.1.0 # - sphinx_rtd_theme=3.1.0 - - mypy #=1.19.* + - mypy - codecov - pytest - pytest-cov @@ -32,6 +32,6 @@ dependencies: - trame-vuetify - pip - pip: -# - sphinx_rtd_theme==3.1.0 + - sphinx_rtd_theme==3.1.0 - --editable=. - git+https://github.com/cadquery/black.git@cq