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]): diff --git a/cadquery/occ_impl/nurbs.py b/cadquery/occ_impl/nurbs.py index da520f81b..035a781f0 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, Sequence from math import comb @@ -23,8 +23,9 @@ 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 +from ..types import Real njit = _njit(cache=True, error_model="numpy", fastmath=True, nogil=True, parallel=False) @@ -81,6 +82,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] @@ -95,10 +106,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): @@ -108,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): """ @@ -180,6 +206,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) @@ -202,6 +231,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()) @@ -289,16 +321,39 @@ 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 @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 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. """ @@ -597,7 +652,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 @@ -661,7 +716,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 @@ -747,14 +804,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): @@ -828,7 +885,7 @@ def nbSurfaceDer( v, vorder, vknots, vperiodic ) - # number of param values + # number of parameter values nu = np.size(u) # chunk sizes @@ -863,7 +920,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 @@ -889,7 +950,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 @@ -906,9 +967,10 @@ 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 + # loop over parameter values for i in range(nu): ui = u_[i] @@ -921,7 +983,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 @@ -943,7 +1005,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 ) @@ -951,7 +1013,7 @@ def designMatrix2D( v, vorder, vknots, vperiodic ) - # number of param values + # number of parameter values ni = len(u) # chunk size @@ -972,9 +1034,14 @@ 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 + # loop over parameter values for i in range(ni): ui, vi = u_[i], v_[i] @@ -989,8 +1056,10 @@ 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 + 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() @@ -1012,7 +1081,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 @@ -1030,10 +1099,11 @@ 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), ) ) - # loop over param values + # loop over parameter values for i in range(nu): ui = u[i] @@ -1063,7 +1133,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 @@ -1084,10 +1154,11 @@ 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), ) ) - # loop over param values + # loop over parameter values for i in range(nu): ui = u[i] @@ -1101,8 +1172,8 @@ 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) - ) % nb # NB: this is due to peridicity + span - order + 1 + np.arange(n) + ) % nb # NB: this is due to periodicity rv[di].v[i * n : (i + 1) * n] = temp[di, :] return rv @@ -1110,6 +1181,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( @@ -1127,6 +1201,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 +1232,10 @@ 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: + """ + Create a discrete penalty matrix for approximating higher order penalties. + """ if order not in (1, 2): raise ValueError( @@ -1175,6 +1253,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: @@ -1262,7 +1341,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 ) @@ -1270,7 +1349,7 @@ def penaltyMatrix2D( v, vorder, vknots, vperiodic ) - # number of param values + # number of parameter values ni = len(u) # chunk size @@ -1286,7 +1365,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( @@ -1294,10 +1373,15 @@ 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), + ), ) ) - # loop over param values + # loop over parameter values for i in range(ni): ui, vi = u_[i], v_[i] @@ -1316,8 +1400,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] @@ -1352,6 +1439,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 = (nv, 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] = 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) + + +@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 = (nu, 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] = ix + 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 @@ -1375,6 +1516,9 @@ def periodicApproximate( penalty: int = 4, lam: float = 0, ) -> Curve: + """ + Approximate a periodic curve. + """ npts = data.shape[0] @@ -1423,7 +1567,7 @@ def periodicApproximate( return rv -@periodicApproximate.register +@multidispatch def periodicApproximate( data: List[Array], us: Optional[Array] = None, @@ -1432,6 +1576,9 @@ def periodicApproximate( penalty: int = 4, lam: float = 0, ) -> List[Curve]: + """ + Approximate multiple periodic curves using same parametrization. + """ rv = [] @@ -1492,6 +1639,9 @@ def approximate( lam: float = 0, tangents: Optional[Tuple[Array, Array]] = None, ) -> Curve: + """ + Approximate a curve. + """ npts = data.shape[0] @@ -1519,7 +1669,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 @@ -1559,7 +1709,7 @@ def approximate( return rv -@approximate.register +@multidispatch def approximate( data: List[Array], us: Optional[Array] = None, @@ -1569,6 +1719,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 = [] @@ -1598,7 +1751,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 @@ -1661,7 +1814,7 @@ def approximate2D( lam: float = 0, ) -> Surface: """ - Simple 2D surface approximation (without any penalty). + Approximate a 2D surface. """ # process the knots @@ -1711,7 +1864,166 @@ def approximate2D( return rv +@multidispatch +def constrainedApproximate2D( + As: COO, + bs: 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: Real = 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)) + C_xyz = sp.block_diag((C, C, C)) + + # solve the augmented normal equations + 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 = 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,) + ), + uknots_, + vknots_, + uorder, + vorder, + uperiodic, + vperiodic, + ) + + return rv + + +@multidispatch +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: Real = 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 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.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: + """ + Periodic loft from curves. + """ nknots: int = len(curves) + 1 @@ -1742,6 +2054,9 @@ def loft( penalty: int = 4, tangents: Optional[List[Tuple[Array, Array]]] = None, ) -> Surface: + """ + Regular (non-periodic) loft form curves. + """ nknots: int = len(curves) diff --git a/cadquery/occ_impl/shapes.py b/cadquery/occ_impl/shapes.py index 7ee861614..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()] @@ -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: 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. + """ + + 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/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() diff --git a/tests/test_nurbs.py b/tests/test_nurbs.py index 5270a18ce..ce984481d 100644 --- a/tests/test_nurbs.py +++ b/tests/test_nurbs.py @@ -12,9 +12,15 @@ periodicApproximate, periodicLoft, loft, + array2vec, + uIsoMatrix, + vIsoMatrix, + vec2array, + constrainedApproximate2D, + COO, ) -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 @@ -275,16 +281,31 @@ def test_approximate(): assert e.Length() == approx(1) -def test_periodic_approximate(): +@mark.parametrize(("penalty", "lam"), [(3, 1e-9), (4, 1e-9), (5, 1e-9), (2, 0)]) +def test_periodic_approximate(penalty, lam): - pts_ = circle(1).sample(100)[0] + 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_) + + # 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]) @@ -293,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): @@ -322,8 +343,7 @@ 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"), [(3, 1e-6), (4, 1e-6), (5, 1e-6), (2, 0)]) def test_approximate2D(lam, penalty): t = torus(5, 1).face() @@ -342,8 +362,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, @@ -352,5 +372,236 @@ 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) + + +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)) + + +@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) + + +@fixture +def torus_face(): + + 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, + 10, + 10, + uperiodic=True, + vperiodic=True, + penalty=3, + lam=1e-6, + ) + + return surf.face() + + +@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("face", FACES) +def test_surface_positions(face, request): + + 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)) + + +@mark.parametrize("face", FACES) +def test_surface_tangents(face, request): + + f = request.getfixturevalue(face) + 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(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) + + +@mark.parametrize("lam", [0.0, 1e-6]) +@mark.parametrize("penalty", [3]) +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, + 10, + 10, + 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) + + # 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.fromSP(A), + b, + pts, + us, + vs, + 3, + 3, + len(surf.uknots), + len(surf.vknots), + uperiodic=True, + vperiodic=True, + penalty=penalty, + lam=lam, + ) 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