From b8e648a8180b3962134f8128acdc630236a3809c Mon Sep 17 00:00:00 2001 From: Kjetil Andre Johannessen Date: Fri, 2 Jan 2026 16:34:33 +0100 Subject: [PATCH 1/2] Fixes #168: const_par_curve snaps to nearest knot --- core/src/splipy_core/__init__.py | 4 ++-- core/src/splipy_core/__init__.pyi | 3 ++- core/src/splipy_core/core.pyx | 32 +++++++++++++++++++++++++++---- src/splipy/basis.py | 14 ++++++++++++-- src/splipy/surface.py | 3 +++ tests/surface_test.py | 17 +++++++++++++++- 6 files changed, 63 insertions(+), 10 deletions(-) diff --git a/core/src/splipy_core/__init__.py b/core/src/splipy_core/__init__.py index e4dbe14..e4208eb 100644 --- a/core/src/splipy_core/__init__.py +++ b/core/src/splipy_core/__init__.py @@ -1,5 +1,5 @@ -from ._core import evaluate, snap +from ._core import evaluate, snap_points, snap_point -__all__ = ["evaluate", "snap"] +__all__ = ["evaluate", "snap_points", "snap_point"] __version__ = "1.10.1" diff --git a/core/src/splipy_core/__init__.pyi b/core/src/splipy_core/__init__.pyi index 9c61335..b9e7cec 100644 --- a/core/src/splipy_core/__init__.pyi +++ b/core/src/splipy_core/__init__.pyi @@ -1,7 +1,8 @@ from numpy import floating, int_ from numpy.typing import NDArray -def snap(knots: NDArray[floating], eval_pts: NDArray[floating], tolerance: float) -> None: ... +def snap_point(knots: NDArray[floating], eval_pt: float, tolerance: float) -> float: ... +def snap_points(knots: NDArray[floating], eval_pts: NDArray[floating], tolerance: float) -> None: ... def evaluate( knots: NDArray[floating], order: int, diff --git a/core/src/splipy_core/core.pyx b/core/src/splipy_core/core.pyx index fc7844e..18d4d60 100644 --- a/core/src/splipy_core/core.pyx +++ b/core/src/splipy_core/core.pyx @@ -59,9 +59,9 @@ def evaluate(cnp.ndarray[cnp.float_t, ndim=1] knots_in, # wrap everything into c-type datastructures for optimized performance cdef cnp.float_t[:] knots = knots_in - cdef unsigned int n_all = len(knots) - p # number of basis functions (without periodicity) - cdef unsigned int n = len(knots) - p - (periodic+1) # number of basis functions (with periodicity) - cdef unsigned int m = len(eval_t_in) + cdef unsigned int n_all = len(knots) - p # number of basis functions (without periodicity) + cdef unsigned int n = len(knots) - p - (periodic+1) # number of basis functions (with periodicity) + cdef unsigned int m = len(eval_t_in) cdef cnp.float_t start = knots[p-1] cdef cnp.float_t end = knots[n_all] cdef cnp.float_t evalT @@ -131,7 +131,7 @@ def evaluate(cnp.ndarray[cnp.float_t, ndim=1] knots_in, @cython.boundscheck(False) @cython.wraparound(False) -def snap(cnp.ndarray[cnp.float_t, ndim=1] knots_in, +def snap_points(cnp.ndarray[cnp.float_t, ndim=1] knots_in, cnp.ndarray[cnp.float_t, ndim=1] eval_t_in, cnp.float_t tolerance): """ Snap evaluation points to knots if they are sufficiently close @@ -153,3 +153,27 @@ def snap(cnp.ndarray[cnp.float_t, ndim=1] knots_in, t[j] = knots[i] elif i > 0 and abs(knots[i-1]-t[j]) < tolerance: t[j] = knots[i-1] + +@cython.boundscheck(False) +@cython.wraparound(False) +def snap_point(cnp.ndarray[cnp.float_t, ndim=1] knots_in, + cnp.float_t eval_t_in, + cnp.float_t tolerance): + """ Snap singular evaluation point to knots if they are sufficiently + close as given in by state.state.knot_tolerance. + + :param knots_in: Knot vector + :param eval_t_in: The parametric coordinate in which to evaluate + :param tolerance_in: Knot tolerance for detecting end-point-evaluations + :return: eval_t_in if sufficiently far from knot. Otherwise closest knot + """ + cdef unsigned int i,j + cdef unsigned int n = len(knots_in) + cdef cnp.float_t t = eval_t_in + cdef cnp.float_t[:] knots = knots_in + i = my_bisect_left(knots, t, n) + if i < n and abs(knots[i]-t) < tolerance: + t = knots[i] + elif i > 0 and abs(knots[i-1]-t) < tolerance: + t = knots[i-1] + return t diff --git a/src/splipy/basis.py b/src/splipy/basis.py index 9c3e723..58bc52a 100644 --- a/src/splipy/basis.py +++ b/src/splipy/basis.py @@ -194,7 +194,7 @@ def evaluate( :rtype: numpy.array """ t_arr = np.atleast_1d(np.asarray(t, dtype=np.float64)) - splipy_core.snap(self.knots, t_arr, state.knot_tolerance) + splipy_core.snap_points(self.knots, t_arr, state.knot_tolerance) if self.order <= d: # requesting more derivatives than polynomial degree: return all zeros return np.zeros((len(t_arr), self.num_functions())) @@ -598,6 +598,16 @@ def matches(self, bspline: BSplineBasis, reverse: bool = False) -> bool: atol=state.knot_tolerance, ) + def snap_point(self, t: float) -> float: + """Snap evaluation point to knots if it is sufficiently close + as given in by state.state.knot_tolerance. + + :param t: evaluation point + :type t: float + :return: float + """ + return splipy_core.snap_point(self.knots, t, state.knot_tolerance) + def snap_points(self, t: FloatArray) -> None: """Snap evaluation points to knots if they are sufficiently close as given in by state.state.knot_tolerance. This will modify the input @@ -607,7 +617,7 @@ def snap_points(self, t: FloatArray) -> None: :type t: [float] :return: none """ - splipy_core.snap(self.knots, t, state.knot_tolerance) + splipy_core.snap_points(self.knots, t, state.knot_tolerance) # TODO(Eivind): Deprecate this later. def snap(self, t): # type: ignore[no-untyped-def] diff --git a/src/splipy/surface.py b/src/splipy/surface.py index 73c4b36..40d55c3 100644 --- a/src/splipy/surface.py +++ b/src/splipy/surface.py @@ -303,6 +303,9 @@ def const_par_curve(self, knot: Scalar, direction: Direction) -> Curve: # clone basis since we need to augment this by knot insertion b = self.bases[direction].clone() + # snap to existing knot if close enough + knot = b.snap_point(knot) + # compute mapping matrix C which is the knot insertion operator mult = b.min_continuity(knot, b.order - 1) C = np.identity(self.shape[direction]) diff --git a/tests/surface_test.py b/tests/surface_test.py index 57e5641..e44e6a2 100644 --- a/tests/surface_test.py +++ b/tests/surface_test.py @@ -6,7 +6,7 @@ import numpy as np import splipy.surface_factory as sf -from splipy import BSplineBasis, Surface +from splipy import BSplineBasis, Curve, Surface class TestSurface(unittest.TestCase): @@ -842,6 +842,21 @@ def test_const_par_crv(self): u = np.linspace(0, 1, 13) self.assertTrue(np.allclose(surf(u, 1.0).reshape(13, 2), crv(u))) + # check near internal knot (issue #168) + bu = BSplineBasis(4) + crv1 = Curve(bu, [[0, 0, 0], [1, 0, 0], [1, 1, 0], [0, 1, 0]], raw=True) + crv2 = Curve(bu, [[0, 0, 3.5], [0.2, 0, 3.5], [0.2, 0.8, 3.5], [0, 0.7, 3.5]], raw=True) + crv3 = Curve(bu, [[0.2, -0.4, 5], [0.2, 1, 5], [0.6, 0.8, 5], [-0.5, 1.7, 5]], raw=True) + crv4 = Curve(bu, [[-0.2, 0.1, 7], [0.28, 1.5, 7], [0.1, -0.8, 7], [0.5, 1.0, 7]], raw=True) + srf = sf.loft(crv1, crv2, crv3, crv4) + srf.bases[1].normalize() + srf.insert_knot(0.85, direction=1) + pt1 = srf(0, 0.85 + 1e-11) + pt2 = srf.const_par_curve(0.85 + 1e-11, direction=1)(0) + self.assertTrue(np.allclose(pt1, pt2)) + pt3 = srf.const_par_curve(0.85 - 1e-11, direction=1)(0) + self.assertTrue(np.allclose(pt1, pt3)) + def test_antiderivative(self): """Test that antiderivative inverts the derivative operation for surfaces.""" From 72313544b7f73cf4b43544a2c152e18e74fe4987 Mon Sep 17 00:00:00 2001 From: Kjetil Andre Johannessen Date: Fri, 2 Jan 2026 16:34:55 +0100 Subject: [PATCH 2/2] Added volume.const_par_surface --- src/splipy/volume.py | 40 ++++++++++++++++++++++++++--- tests/volume_test.py | 61 ++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 98 insertions(+), 3 deletions(-) diff --git a/src/splipy/volume.py b/src/splipy/volume.py index 22ddd20..0ddad5c 100644 --- a/src/splipy/volume.py +++ b/src/splipy/volume.py @@ -1,19 +1,20 @@ from __future__ import annotations +from bisect import bisect_left from typing import TYPE_CHECKING, ClassVar, cast import numpy as np from .basis import BSplineBasis from .splineobject import SplineObject -from .utils import ensure_listlike, sections +from .surface import Surface +from .utils import check_direction, ensure_listlike, sections if TYPE_CHECKING: from collections.abc import Sequence from .curve import Curve - from .surface import Surface - from .typing import ArrayLike, Direction, FloatArray + from .typing import ArrayLike, Direction, FloatArray, Scalar __all__ = ["Volume"] @@ -100,6 +101,39 @@ def faces( tuple(boundary_faces), ) + def const_par_surface(self, knot: Scalar, direction: Direction) -> Surface: + """Get a Surface representation of the parametric volume at some constant + knot value. + :param float knot: The constant knot value to sample the volume + :param int direction: The parametric direction for the constant value + :return: surface in this volume + :rtype: Surface + """ + direction = check_direction(direction, 3) + + # clone basis since we need to augment this by knot insertion + b = self.bases[direction].clone() + + # snap to existing knot if close enough + knot = b.snap_point(knot) + + # compute mapping matrix C which is the knot insertion operator + mult = b.min_continuity(knot, b.order - 1) + C = np.identity(self.shape[direction]) + for i in range(mult): + C = b.insert_knot(knot) @ C + + # at this point we have a C0 basis, find the right interpolating index + i = max(bisect_left(b.knots, knot) - 1, 0) + + # compute the controlpoints and return Surface + cp = np.tensordot(C[i, :], self.controlpoints, axes=(0, direction)) + cp = cp.transpose(1, 0, 2) + + # return surface with the two remaining directions + remaining_dirs = [d for d in range(3) if d != direction] + return Surface(self.bases[remaining_dirs[0]], self.bases[remaining_dirs[1]], cp, self.rational) + def volume(self) -> float: """Computes the volume of the object in geometric space""" diff --git a/tests/volume_test.py b/tests/volume_test.py index 30b2c47..68f6b1e 100644 --- a/tests/volume_test.py +++ b/tests/volume_test.py @@ -752,6 +752,67 @@ def test_volume(self): v[:, :, 1, 0:2] = 0.0 # squeeze top together, creating a pyramid self.assertAlmostEqual(v.volume(), 1.0 / 3) + def test_const_par_surface(self): + # test with a unit cube + vol = Volume() + + # test u-direction (direction 0) + surf = vol.const_par_surface(0.5, 0) + print(vol) + print(surf) + self.assertEqual(surf.pardim, 2) # should be a surface + self.assertEqual(surf.dimension, 3) # 3D geometry + # evaluate at the center of the resulting surface + pt = surf(0.1, 0.2) + self.assertAlmostEqual(pt[0], 0.5) # u is fixed at 0.5 + self.assertAlmostEqual(pt[1], 0.1) # v varies + self.assertAlmostEqual(pt[2], 0.2) # w varies + + # test v-direction (direction 1) + surf = vol.const_par_surface(0.25, 1) + self.assertEqual(surf.pardim, 2) + pt = surf(0.3, 0.4) + self.assertAlmostEqual(pt[0], 0.3) # u varies + self.assertAlmostEqual(pt[1], 0.25) # v is fixed at 0.25 + self.assertAlmostEqual(pt[2], 0.4) # w varies + + # test w-direction (direction 2) + surf = vol.const_par_surface(0.75, 2) + self.assertEqual(surf.pardim, 2) + pt = surf(0.9, 0.7) + self.assertAlmostEqual(pt[0], 0.9) # u varies + self.assertAlmostEqual(pt[1], 0.7) # v varies + self.assertAlmostEqual(pt[2], 0.75) # w is fixed at 0.75 + + # test with string directions + surf = vol.const_par_surface(0.3, "u") + pt = surf(0.5, 0.5) + self.assertAlmostEqual(pt[0], 0.3) + + surf = vol.const_par_surface(0.7, "v") + pt = surf(0.5, 0.5) + self.assertAlmostEqual(pt[1], 0.7) + + surf = vol.const_par_surface(0.9, "w") + pt = surf(0.5, 0.5) + self.assertAlmostEqual(pt[2], 0.9) + + b1 = BSplineBasis(2) + b2 = BSplineBasis(3) + b3 = BSplineBasis(4) + vol = Volume(b1, b2, b3) + surf = vol.const_par_surface(0.5, 0) + self.assertEqual(surf.bases[0].order, 3) + self.assertEqual(surf.bases[1].order, 4) + + surf = vol.const_par_surface(0.3, 1) + self.assertEqual(surf.bases[0].order, 2) + self.assertEqual(surf.bases[1].order, 4) + + surf = vol.const_par_surface(0.8, 2) + self.assertEqual(surf.bases[0].order, 2) + self.assertEqual(surf.bases[1].order, 3) + def test_operators(self): v = Volume() v.raise_order(1, 1, 2)