Perform knot snap on surface.const_par_curve

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"

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,

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

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]

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])

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"""
 

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."""
 

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)