Repository: EPFL-Center-for-Imaging/splinebox

Issue:
Rendering splines as tube
Spline box should include an easy way to render a spline in 2D or 3D as a tube.

The method of the spline class should take the `r`, `shape`, `value` as parameters.
The radius `r` of the spline should be either a float or a function of the spline parameter `t` for a tube with varying radius.
`shape` specifies the  number of pixels of the output image.

Possible approaches:

1. For a constant radius one could just create a binary skeleton of the spline and compute the euclidean distance transform. This will cause some problems at the end.
2. Construct a volume mesh with tetrahydra using cylinders and compute which pixels are inside the mesh using the Jacobian of the mesh. A similar approach is implemented in [spam](https://www.spam-project.dev/docs/spam.mesh.html#submodules). Using the bounding box of the mesh pixels can be excluded to speed up the computation.
3. Construct a surface mesh. If the mesh is constructed correctly it should be possible to use the normal of the closest face to determine if a point is on the inside or on the outside.

Discussion:
@faymanns: Implemented in #38 

Solution patch:
diff --git a/examples/plot_frame.py b/examples/plot_frame.py
new file mode 100644
index 0000000..0847291
--- /dev/null
+++ b/examples/plot_frame.py
@@ -0,0 +1,97 @@
+"""
+Moving Frames
+=============
+
+`Moving frames`_ provide a local coordinate system along a curve. This system can be used to analyze the space around points on the curve, such as extracting orthogonal images from a volume. An example of this application is the Dendrite-Centric Coordinate System.
+
+SplineBox implements two types of moving frames:
+
+1. `Frenet-Serret Frame`_: Defined by the curve's tangent, normal, and binormal vectors. To compute all three vectors, the curve cannot have any straight segments and no inflection points. Additionally, the frame may rotate around the curve.
+2. `Bishop Frame`_: A twist-free alternative that eliminates rotations around the curve and is defined on straight segments and at inflections points [Bishop1975]_.
+
+.. _Moving frames: https://en.wikipedia.org/wiki/Moving_frame
+.. _Frenet-Serret Frame: https://en.wikipedia.org/wiki/Frenet-Serret_formulas
+.. _Bishop Frame: https://www.jstor.org/stable/2319846
+"""
+
+import numpy as np
+import pyvista
+import splinebox
+
+# %%
+# 1. Create a Spline
+# ^^^^^^^^^^^^^^^^^^
+
+spline = splinebox.Spline(M=4, basis_function=splinebox.B3(), closed=False)
+spline.control_points = np.array(
+    [
+        [0.0, 0.0, 0.0],
+        [1.0, 0.0, 0.0],
+        [1.0, 1.0, 0.0],
+        [1.0, 1.0, 1.0],
+        [0.0, 1.0, 1.0],
+        [0.0, 0.0, 0.0],
+    ]
+)
+
+t = np.linspace(0, spline.M - 1, spline.M * 3)
+
+# %%
+# 2. Compute the Frenet-Serret Frame
+# ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+frenet_frame = spline.moving_frame(t, method="frenet")
+
+# %%
+# 3. Compute the Bishop Frame
+# ^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+bishop_frame = spline.moving_frame(t, method="bishop")
+
+# %%
+# 4. Visualize the Frames
+# ^^^^^^^^^^^^^^^^^^^^^^^
+# We can visualize the Frenet-Serret and Bishop frames using PyVista.
+# The following code plots the two frames side-by-side to highlight their differences.
+# The Frenet-Serret frame (left) visibly twists along the curve, while the Bishop frame (right) avoids this twist.
+
+# Create a PyVista mesh for the curve
+spline_mesh = pyvista.MultipleLines(points=spline.eval(t))
+
+# Add vectors to the mesh for visualization
+spline_mesh["frenet0"] = frenet_frame[:, 0] * 0.2
+spline_mesh["frenet1"] = frenet_frame[:, 1] * 0.2
+spline_mesh["frenet2"] = frenet_frame[:, 2] * 0.2
+spline_mesh["bishop0"] = bishop_frame[:, 0] * 0.2
+spline_mesh["bishop1"] = bishop_frame[:, 1] * 0.2
+spline_mesh["bishop2"] = bishop_frame[:, 2] * 0.2
+
+# Initialize a PyVista plotter
+plotter = pyvista.Plotter(shape=(1, 2), border=False)
+
+# Plot the Frenet-Serret frame
+plotter.subplot(0, 0)
+spline_mesh.set_active_vectors("frenet0")
+plotter.add_mesh(spline_mesh.arrows, lighting=False, color="blue")
+spline_mesh.set_active_vectors("frenet1")
+plotter.add_mesh(spline_mesh.arrows, lighting=False, color="red")
+spline_mesh.set_active_vectors("frenet2")
+plotter.add_mesh(spline_mesh.arrows, lighting=False, color="green")
+plotter.add_mesh(spline_mesh, line_width=10, color="black")
+
+# Plot the Bishop frame
+plotter.subplot(0, 1)
+spline_mesh.set_active_vectors("bishop0")
+plotter.add_mesh(spline_mesh.arrows, lighting=False, color="blue")
+spline_mesh.set_active_vectors("bishop1")
+plotter.add_mesh(spline_mesh.arrows, lighting=False, color="red")
+spline_mesh.set_active_vectors("bishop2")
+plotter.add_mesh(spline_mesh.arrows, lighting=False, color="green")
+plotter.add_mesh(spline_mesh, line_width=10, color="black", copy_mesh=True)
+
+plotter.link_views()
+plotter.show()
+
+# %%
+# If we want the Bishop frame to start in a specific orientation, we can specify
+# an :code:`initial_vector` in :meth:`splinebox.spline_curves.Spline.moving_frame()`.
diff --git a/examples/plot_mesh.py b/examples/plot_mesh.py
new file mode 100644
index 0000000..c19c115
--- /dev/null
+++ b/examples/plot_mesh.py
@@ -0,0 +1,171 @@
+"""
+Spline to mesh
+==============
+This guide demonstrates how to convert a 3D spline curve into various types of meshes using SplineBox.
+"""
+
+# sphinx_gallery_thumbnail_number = 4
+
+import numpy as np
+import pyvista
+import splinebox
+
+# %%
+# 1. Constructing a Spline
+# ------------------------
+# We begin by creating a circular spline.
+M = 4
+spline = splinebox.Spline(M=M, basis_function=splinebox.Exponential(M), closed=True)
+spline.knots = np.array([[0, 0, 1], [0, 1, 0], [0, 0, -1], [0, -1, 0]])
+
+# %%
+# 2. Generating a Mesh Without Radius
+# -----------------------------------
+# The :code:`step_t` parameter determines the granularity of the resulting mesh, corresponding to the step size in the spline parameter space (t).
+# Setting the radius to None or 0 results in a line mesh.
+
+# Generate a simple line mesh
+points, connectivity = spline.mesh(step_t=0.1, radius=None)
+
+# Prepend the number of points in each element (2 for a line) for PyVista
+connectivity = np.hstack((np.full((connectivity.shape[0], 1), 2), connectivity))
+
+# Create and plot the PyVista mesh
+mesh = pyvista.PolyData(points, lines=connectivity)
+mesh.plot()
+
+# %%
+# 3. Mesh with a Fixed Radius
+# ---------------------------
+# Here, we generate a surface mesh (a "tube") using a fixed radius.
+# We employ the Frenet-Serret frame to avoid selecting an initial vector.
+
+# Generate a surface mesh with a fixed radius
+points, connectivity = spline.mesh(step_t=0.1, radius=0.2, frame="frenet")
+
+# Prepend the number of points in each element (3 for triangles) for PyVista
+connectivity = np.hstack((np.full((connectivity.shape[0], 1), 3), connectivity))
+
+# Create and plot the PyVista mesh
+mesh = pyvista.PolyData(points, faces=connectivity)
+mesh.plot(show_edges=True)
+
+# %%
+# 3. Mesh with an Elliptical Cross-Section
+# ----------------------------------------
+# You can define a custom cross-section shape by specifying the radius as a function of the spline parameter (:code:`t`) and the polar angle (:code:`phi`).
+# Example 1: Elliptical Cross-Section
+# ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+
+def elliptical_radius(t, phi):
+    a = 0.1
+    b = 0.05
+    phi = np.deg2rad(phi)
+    r = (a * b) / np.sqrt((b * np.cos(phi)) ** 2 + (a * np.sin(phi)) ** 2)
+    return r
+
+
+points, connectivity = spline.mesh(step_t=0.1, step_angle=36, radius=elliptical_radius, frame="frenet")
+
+connectivity = np.hstack((np.full((connectivity.shape[0], 1), 3), connectivity))
+
+mesh = pyvista.PolyData(points, faces=connectivity)
+mesh.plot(show_edges=True)
+
+# %%
+# Example 2: Varying Radius Along the Spline
+# ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+
+def radius(t, phi):
+    return 0.1 + 0.03 * np.sin(t / spline.M * 16 * np.pi)
+
+
+points, connectivity = spline.mesh(step_t=0.1, step_angle=36, radius=radius, frame="frenet")
+
+connectivity = np.hstack((np.full((connectivity.shape[0], 1), 3), connectivity))
+
+mesh = pyvista.PolyData(points, faces=connectivity)
+mesh.plot(show_edges=True)
+
+# %%
+# 5. Bishop Frame for Mesh Generation
+# -----------------------------------
+# The Frenet-Serret frame is not defined on straight segments and at inflections point.
+# In those cases, we can use the Bishop frame instead. Another advantage of the Bishop frame
+# is that it does not twist around the spline.
+#
+# Correcting Twists with the Bishop Frame
+# ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+# Create a spline with for which the Frenet frame twists
+spline = splinebox.Spline(M=M, basis_function=splinebox.B3(), closed=False)
+spline.control_points = np.array(
+    [
+        [0.0, 0.0, 0.0],
+        [2.0, 0.0, 0.0],
+        [2.0, 2.0, 0.0],
+        [2.0, 2.0, 2.0],
+        [0.0, 2.0, 2.0],
+        [0.0, 0.0, 0.0],
+    ]
+)
+
+points, connectivity = spline.mesh(step_t=0.1, step_angle=36, radius=elliptical_radius, frame="frenet")
+
+connectivity = np.hstack((np.full((connectivity.shape[0], 1), 3), connectivity))
+
+mesh = pyvista.PolyData(points, faces=connectivity)
+mesh.plot(show_edges=True)
+
+# %%
+# You can clearly see how the ellipse twists around the spline.
+# The Bishop frame eliminates this twist.
+
+points, connectivity = spline.mesh(step_t=0.1, step_angle=36, radius=elliptical_radius, frame="bishop")
+
+connectivity = np.hstack((np.full((connectivity.shape[0], 1), 3), connectivity))
+
+mesh = pyvista.PolyData(points, faces=connectivity)
+mesh.plot(show_edges=True)
+
+# %%
+# Changing the initial vector rotates the frame and therefore the ellipse.
+
+initial_vector = np.array([0.5, -0.5, 1])
+
+points, connectivity = spline.mesh(
+    step_t=0.1, step_angle=36, radius=elliptical_radius, frame="bishop", initial_vector=initial_vector
+)
+
+connectivity = np.hstack((np.full((connectivity.shape[0], 1), 3), connectivity))
+
+mesh = pyvista.PolyData(points, faces=connectivity)
+mesh.plot(show_edges=True)
+
+# %%
+# 6. Volume Mesh
+# --------------
+# Finally, you can generate a volumetric mesh by setting the :code:`mesh_type` to :code:`"volume"`.
+
+points, connectivity = spline.mesh(
+    radius=radius, step_t=0.5, step_angle=72, initial_vector=initial_vector, mesh_type="volume"
+)
+
+connectivity = np.hstack((np.full((connectivity.shape[0], 1), 4), connectivity))
+cell_types = np.full(len(connectivity), fill_value=pyvista.CellType.TETRA, dtype=np.uint8)
+
+mesh = pyvista.UnstructuredGrid(connectivity, cell_types, points)
+mesh.plot(show_edges=True)
+
+# %%
+# For a better understanding of the volume mesh we can explode it.
+# This allows us to see the individual tetrahedra.
+mesh.explode(factor=0.5).plot(show_edges=True)
+
+# %%
+# Tips
+# ----
+# * Save meshes for visualization in ParaView using :code:`mesh.save("mesh.vtk")`.
+# * Surface meshes are open by default. Use the :code:`cap_ends` keyword in :meth:`splinebox.spline_curves.Spline.mesh()` to close them.
diff --git a/src/splinebox/spline_curves.py b/src/splinebox/spline_curves.py
index 3634978..d200d2f 100644
--- a/src/splinebox/spline_curves.py
+++ b/src/splinebox/spline_curves.py
@@ -641,7 +641,7 @@ def normal(self, t, frame="bishop", initial_vector=None):
             normals = (np.array([[0, -1], [1, 0]]) @ first_deriv.T).T
             normals /= np.linalg.norm(normals, axis=1)[:, np.newaxis]
         elif self.control_points.shape[1] == 3:
-            frame = self.moving_frame(t, kind=frame, initial_vector=initial_vector)
+            frame = self.moving_frame(t, method=frame, initial_vector=initial_vector)
             normals = frame[:, 1:]
         else:
             raise RuntimeError(
@@ -649,35 +649,67 @@ def normal(self, t, frame="bishop", initial_vector=None):
             )
         return normals
 
-    def moving_frame(self, t, kind="frenet", initial_vector=None):
+    def moving_frame(self, t, method="frenet", initial_vector=None):
         """
-        This function defines two of the `moving frames` on the spline curve,
-        the Frenet-Serre frame and the Bishop frame. It returns an orthonormal
-        basis for each t. The Bishop frame [Bishop1975]_ is constructed such that
-        it does not twist around the curve, i.e. it has zero torsion.
+        Compute a `moving frame`_ (local orthonormal coordinate system) along the spline.
+
+        This method computes either the Frenet-Serret frame or the Bishop frame for
+        the spline. A moving frame consists of three orthonormal basis vectors at
+        each point on the curve. The Frenet-Serret frame is derived from the curve's
+        derivatives but may twist around the curve. The Bishop frame eliminates
+        this twist, providing a zero-torsion alternative.
 
         Parameters
         ----------
         t : np.array
-            The parameters values of the spline for which the frame should be
-            evaluated.
-        kind : string
-            The type of frame. Can be "frenet" or "bishop".
-        initial_vector : np.array or None
-            The initial vector for the Bishop frame. It has to be orthogonal
-            to the tangent at :code:`t[0]` and defines the orientation of the basis
-            at :code:`t[0]`. This orientation is then propagated along the spline
-            since the Bishop frame does not allow the basis to twist around
-            the curve.
+            A 1D array of parameter values at which to evaluate the frame. These
+            correspond to positions along the spline.
+        method : str, optional
+            The type of moving frame to compute. Options are:
+
+            - "frenet": The classical Frenet-Serret frame, based on tangent, normal, and binormal vectors.
+            - "bishop": A twist-free frame that requires an initial orientation.
+
+            Default is "frenet".
+        initial_vector : np.array or None, optional
+            For the Bishop frame, an initial vector that is orthogonal to the tangent
+            vector at `t[0]`. This vector determines the initial orientation of the
+            basis, which is propagated along the curve without twisting. If None,
+            the method computes a suitable initial vector automatically. This
+            parameter is ignored when :code:`method="frenet"`.
 
         Returns
         -------
         frame : np.array
-            A numpy array with 3 dimensions. The first dimension corresponds
-            to t, the second to the 3 basis vectors and the last one are the
-            components of each basis vector.
-
-        .. _moving frames: https://en.wikipedia.org/wiki/Moving_frame
+            A 3D numpy array with shape `(len(t), 3, 3)`. The dimensions are:
+
+            - The first axis corresponds to the parameter values in `t`.
+            - The second axis contains the three basis vectors at each `t`:
+              [tangent, normal, binormal] for "frenet" or equivalent vectors for "bishop".
+            - The third axis contains the components of each basis vector in 3D space.
+
+        Raises
+        ------
+        RuntimeError
+            If the spline is not defined in 3D or if the Frenet frame cannot be
+            computed due to inflection points, straight segments, or undefined
+            tangent/normal vectors.
+        ValueError
+            If the initial vector for the Bishop frame is not orthogonal to the
+            tangent at `t[0]`, or if an invalid `method` is specified.
+
+        Notes
+        -----
+        - The Frenet frame is not defined at points where the curve has zero curvature,
+          such as straight segments or inflection points. In these cases, the Bishop
+          frame is recommended.
+        - For closed curves, check for discontinuities of the Bishop frame.
+
+        References
+        ----------
+        .. [1] Bishop, R. L. (1975). "There is More than One Way to Frame a Curve."
+               American Mathematical Monthly, 82(3), 246-251.
+        .. _moving frame: https://en.wikipedia.org/wiki/Moving_frame
         """
         self._check_control_points()
         if self.control_points.ndim != 2 or self.control_points.shape[1] != 3:
@@ -687,7 +719,7 @@ def moving_frame(self, t, kind="frenet", initial_vector=None):
         frame = np.zeros((len(t), 3, 3))
         frame[:, 0] = first_derivative / np.linalg.norm(first_derivative, axis=-1)[:, np.newaxis]
 
-        if kind == "frenet":
+        if method == "frenet":
             second_derivative = self.eval(t, derivative=2)
             frame[:, 2] = np.cross(first_derivative, second_derivative)
             norm_binormal = np.linalg.norm(frame[:, 2], axis=-1)[:, np.newaxis]
@@ -701,9 +733,17 @@ def moving_frame(self, t, kind="frenet", initial_vector=None):
                 )
             frame[:, 2] /= norm_binormal
             frame[:, 1] = np.cross(frame[:, 2], frame[:, 0])
-        elif kind == "bishop":
+        elif method == "bishop":
             if initial_vector is None:
-                raise ValueError("The bishop frame requires the keyword argument initial_vector.")
+                tangent = frame[0, 0]
+                # Try to do the same as for the Frenet frame
+                initial_vector = np.cross(np.cross(tangent, self.eval(t[0], derivative=2)), tangent)
+                if np.isclose(np.linalg.norm(initial_vector), 0) or np.any(np.isnan(initial_vector)):
+                    initial_vector = np.zeros(3)
+                    max_axis = np.argmax(np.abs(tangent))
+                    other_axis = (max_axis + 1) % 3
+                    initial_vector[max_axis] = tangent[other_axis]
+                    initial_vector[other_axis] = -tangent[max_axis]
             initial_vector /= np.linalg.norm(initial_vector)
             if not np.isclose(np.dot(frame[0, 0], initial_vector), 0):
                 raise ValueError("The initial vector has to be orthogonal to the tangent at t[0].")
@@ -729,7 +769,7 @@ def moving_frame(self, t, kind="frenet", initial_vector=None):
                         + n * np.dot(n, frame[i - 1, 2]) * (1 - np.cos(phi))
                     )
         else:
-            raise ValueError(f"Unkown kind of frame {kind}.")
+            raise ValueError(f"Unkown method '{method}' for moving frame.")
         return frame
 
     def eval(self, t, derivative=0):
@@ -906,6 +946,231 @@ def _distance(t):
 
         return min_distance
 
+    def mesh(
+        self,
+        radius=None,
+        step_t=0.1,
+        step_angle=5,
+        mesh_type="surface",
+        cap_ends=False,
+        frame="bishop",
+        initial_vector=None,
+    ):
+        """
+        Create a 3D mesh around the spline curve.
+
+        This method generates a mesh surrounding the spline, with the distance from
+        the spline controlled by the `radius` parameter. The mesh can be either a
+        surface mesh with triangular cells or a volume mesh with tetrahedral cells.
+        The orientation of the mesh is determined using either the Frenet-Serret or
+        Bishop frame.
+
+        Parameters
+        ----------
+        radius : float or callable, optional
+            Defines the distance from the spline to the mesh surface. Can be:
+
+            - A float for a constant radius.
+            - A callable function that takes the spline parameter `t` and the polar
+              angle in the normal plane as arguments and returns a float.
+
+            If None, the mesh will follow the spline without an offset. Default is None.
+        step_t : float, optional
+            Step size for the spline parameter `t`, controlling the resolution along
+            the curve. Smaller values increase mesh resolution. Default is 0.1.
+        step_angle : float, optional
+            Step size for the polar angle (in degrees), controlling the resolution
+            around the spline. Smaller values increase mesh resolution. Default is 5.
+        mesh_type : str, optional
+            Specifies the type of mesh to create:
+
+            - "surface": A surface mesh with triangular cells.
+            - "volume": A volume mesh with tetrahedral cells.
+
+            Default is "surface".
+        cap_ends : bool, optional
+            If True, the ends of an open surface mesh are capped with orthogonal
+            planes. Ignored for closed splines or volume meshes. Default is False.
+        frame : str, optional
+            The frame to use for orientation of the mesh:
+            - "frenet": Uses the Frenet-Serret frame.
+            - "bishop": Uses the Bishop frame, requiring an `initial_vector`.
+            See :meth:`splinebox.spline_curves.moving_frame`. Default is "bishop".
+        initial_vector : numpy array or None, optional
+            For the Bishop frame, an initial vector that defines the orientation of
+            the frame at the start of the spline (`t[0]`). This vector must be
+            orthogonal to the tangent at `t[0]`. Ignored for the Frenet frame. If
+            None, a suitable initial vector is computed automatically. Default is None.
+
+        Returns
+        -------
+        points : numpy array
+            A 2D array of shape `(N, 3)`, where `N` is the number of mesh points.
+            Each row represents the (x, y, z) coordinates of a point in the mesh.
+        connectivity : numpy array
+            A 2D array of shape `(M, K)`, where `M` is the number of elements in the
+            mesh and `K` is the number of vertices per element (3 for surface meshes
+            and 4 for volume meshes). Each row contains the indices of `points` that
+            form an element.
+
+        Raises
+        ------
+        NotImplementedError
+            If the spline is not defined in 3D, as meshes are only supported for 3D splines.
+
+        Notes
+        -----
+        - Surface meshes are useful for visualization, while volume meshes are
+          typically used in simulations and finite element analysis.
+        - For open splines, capping the ends (`cap_ends=True`) creates closed surfaces,
+          which may be useful for some applications.
+        - The Bishop frame is recommended for curves with inflection points or
+          straight segments where the Frenet frame is undefined.
+        - When radius is callable, the Bishop frame is recommend to avoid "drift" of the
+          polar angle around the curve.
+
+        Examples
+        --------
+        Create a surface mesh with constant radius:
+
+        >>> points, connectivity = spline.mesh(radius=0.5, step_t=0.1, step_angle=10, mesh_type="surface")
+
+        Create a volume mesh with variable radius:
+
+        >>> def radius_function(t, angle):
+        >>>     return 0.5 + 0.2 * np.sin(np.radians(angle))
+        >>> points, connectivity = spline.mesh(radius=radius_function, mesh_type="volume")
+        """
+        if self.control_points.ndim != 2 or self.control_points.shape[1] != 3:
+            raise NotImplementedError("Meshes are only implemented for splines in 3D.")
+        t = np.arange(0, self.M if self.closed else self.M - 1 + step_t, step_t)
+        if radius is None or radius == 0:
+            points = self.eval(t)
+            connectivity = np.stack((np.arange(len(points)), np.arange(len(points)) + 1), axis=-1)
+            if self.closed:
+                # Connect end to the beginning
+                connectivity[-1, -1] = 0
+            else:
+                connectivity = connectivity[:-1]
+        else:
+            _radius = (lambda t, phi: np.full(t.shape, radius)) if not callable(radius) else radius
+
+            phi = np.arange(0, 360, step_angle)
+
+            if mesh_type == "surface":
+                phiphi, tt = np.meshgrid(phi, t)
+                tt = tt.flatten()
+                phiphi = phiphi.flatten()
+                centers = self.eval(tt.flatten())
+                rr = _radius(tt, phiphi)
+                normals = self.normal(t, frame=frame, initial_vector=initial_vector)
+
+                n_angles = len(phi)
+                n_t = len(t)
+                normals = (
+                    np.repeat(normals[:, 0], n_angles, axis=0) * np.sin(np.deg2rad(phiphi))[:, np.newaxis]
+                    + np.repeat(normals[:, 1], n_angles, axis=0) * np.cos(np.deg2rad(phiphi))[:, np.newaxis]
+                )
+
+                points = centers + rr[:, np.newaxis] * normals
+                n_points = len(points)
+                connectivity = self._surface_mesh_connectivity(self.closed, n_angles, n_t, n_points)
+                if cap_ends and not self.closed:
+                    points = np.concatenate((centers[0].reshape((1, -1)), points, centers[-1].reshape((1, -1))), axis=0)
+                    start_connectivity = np.zeros((n_angles, 3), dtype=int)
+                    start_connectivity[:, 1] = np.arange(1, n_angles + 1)
+                    start_connectivity[:, 2] = np.roll(start_connectivity[:, 1], -1)
+                    end_connectivity = np.zeros((n_angles, 3), dtype=int)
+                    end_connectivity[:, 0] = n_points + 1
+                    end_connectivity[:, 1] = np.arange(n_points, n_points - n_angles, -1)
+                    end_connectivity[:, 2] = np.roll(end_connectivity[:, 1], -1)
+                    connectivity = np.concatenate((start_connectivity, connectivity + 1, end_connectivity))
+
+            elif mesh_type == "volume":
+                phiphi, tt = np.meshgrid(phi, t)
+                rr = _radius(tt, phiphi)
+                # Add columns for the center points
+                rr = np.hstack((np.zeros((rr.shape[0], 1)), rr))
+                tt = np.hstack((tt[:, 0][:, np.newaxis], tt))
+                phiphi = np.hstack((phiphi[:, 0][:, np.newaxis], phiphi))
+
+                tt = tt.flatten()
+                phiphi = phiphi.flatten()
+                rr = rr.flatten()
+
+                centers = self.eval(tt.flatten())
+                normals = self.normal(t, frame=frame, initial_vector=initial_vector)
+
+                n_angles = len(phi)
+                n_t = len(t)
+                normals = (
+                    np.repeat(normals[:, 0], n_angles + 1, axis=0) * np.sin(np.deg2rad(phiphi))[:, np.newaxis]
+                    + np.repeat(normals[:, 1], n_angles + 1, axis=0) * np.cos(np.deg2rad(phiphi))[:, np.newaxis]
+                )
+                points = centers + rr[:, np.newaxis] * normals
+                n_points = len(points)
+
+                connectivity = self._volume_mesh_connectivity(self.closed, n_angles, n_t, n_points)
+
+        return points, connectivity
+
+    @staticmethod
+    @numba.jit(nopython=True, nogil=True, cache=True)
+    def _surface_mesh_connectivity(closed, n_angles, n_t, n_points):  # pragma: no cover
+        if closed:
+            connectivity = np.zeros((2 * n_angles * n_t, 3), dtype=numba.int64)
+        else:
+            connectivity = np.zeros((2 * n_angles * (n_t - 1), 3), dtype=numba.int64)
+        face = 0
+        for i in range(n_t if closed else n_t - 1):
+            for j in range(n_angles):
+                connectivity[face] = [
+                    i * n_angles + j,
+                    ((i + 1) * n_angles + j) % n_points,
+                    ((i + 1) * n_angles + (j + 1) % n_angles) % n_points,
+                ]
+                face += 1
+                connectivity[face] = [
+                    i * n_angles + j,
+                    ((i + 1) * n_angles + (j + 1) % n_angles) % n_points,
+                    (i * n_angles + (j + 1) % n_angles) % n_points,
+                ]
+                face += 1
+        return connectivity
+
+    @staticmethod
+    @numba.jit(nopython=True, nogil=True, cache=True)
+    def _volume_mesh_connectivity(closed, n_angles, n_t, n_points):  # pragma: no cover
+        if closed:
+            connectivity = np.zeros((3 * n_angles * n_t, 4), dtype=numba.int64)
+        else:
+            connectivity = np.zeros((3 * n_angles * (n_t - 1), 4), dtype=numba.int64)
+        vol = 0
+        for i in range(n_t if closed else n_t - 1):
+            for j in range(1, n_angles + 1):
+                connectivity[vol] = [
+                    i * (n_angles + 1) + j,
+                    ((i + 1) * (n_angles + 1) + j) % n_points,
+                    ((i + 1) * (n_angles + 1) + 1 + j % n_angles) % n_points,
+                    (i + 1) * (n_angles + 1) % n_points,
+                ]
+                vol += 1
+                connectivity[vol] = [
+                    i * (n_angles + 1) + j,
+                    ((i + 1) * (n_angles + 1) + 1 + j % n_angles) % n_points,
+                    (i * (n_angles + 1) + 1 + j % n_angles) % n_points,
+                    i * (n_angles + 1),
+                ]
+                vol += 1
+                connectivity[vol] = [
+                    i * (n_angles + 1) + j,
+                    ((i + 1) * (n_angles + 1) + 1 + j % n_angles) % n_points,
+                    (i + 1) * (n_angles + 1) % n_points,
+                    i * (n_angles + 1),
+                ]
+                vol += 1
+        return connectivity
+
 
 class HermiteSpline(Spline):
     """


Test patch:
diff --git a/tests/conftest.py b/tests/conftest.py
index 9eb87ed..3682f7a 100644
--- a/tests/conftest.py
+++ b/tests/conftest.py
@@ -295,3 +295,41 @@ def multigenerator(request):
 @pytest.fixture(params=[2, 3, 4])
 def support(request):
     return request.param
+
+
+def radius_func(t, phi):
+    a = 0.1 + 0.03 * np.sin(t * 4 * np.pi)
+    b = 0.05
+    phi = np.deg2rad(phi)
+    r = (a * b) / np.sqrt((b * np.cos(phi)) ** 2 + (a * np.sin(phi)) ** 2)
+    return r
+
+
+@pytest.fixture(params=[None, 1, 1.5, radius_func])
+def radius(request):
+    return request.param
+
+
+@pytest.fixture(params=[0.1, 1])
+def step_t(request):
+    return request.param
+
+
+@pytest.fixture(params=[5, 25])
+def step_angle(request):
+    return request.param
+
+
+@pytest.fixture(params=["surface", "volume"])
+def mesh_type(request):
+    return request.param
+
+
+@pytest.fixture(params=[False, True])
+def cap_ends(request):
+    return request.param
+
+
+@pytest.fixture(params=["frenet", "bishop"])
+def frame(request):
+    return request.param
diff --git a/tests/test_spline_curves.py b/tests/test_spline_curves.py
index 6e2228b..4c1dcff 100644
--- a/tests/test_spline_curves.py
+++ b/tests/test_spline_curves.py
@@ -773,9 +773,9 @@ def test_moving_frame(initialized_spline_curve, not_differentiable_twice):
     else:
         if not_differentiable_twice(spline):
             with pytest.raises(RuntimeError):
-                frame = spline.moving_frame(t, kind="frenet")
+                frame = spline.moving_frame(t, method="frenet")
         else:
-            frame = spline.moving_frame(t, kind="frenet")
+            frame = spline.moving_frame(t, method="frenet")
 
             # Check that they are unit vectors
             assert np.allclose(np.linalg.norm(frame, axis=-1), 1)
@@ -798,19 +798,15 @@ def test_moving_frame(initialized_spline_curve, not_differentiable_twice):
             deriv2 -= np.sum(deriv2 * frame[:, 0], axis=-1)[:, np.newaxis] * frame[:, 0]
             assert np.allclose(np.cross(frame[:, 1], deriv2), 0)
 
-        with pytest.raises(ValueError):
-            # Missing initial_vector
-            frame = spline.moving_frame(t, kind="bishop")
-
         with pytest.raises(ValueError):
             # Non-orthogonal initial vector
-            frame = spline.moving_frame(t, kind="bishop", initial_vector=np.ones(3))
+            frame = spline.moving_frame(t, method="bishop", initial_vector=np.ones(3))
 
         initial_vector = np.zeros(3)
         tangent = spline.eval(t[0], derivative=1)
         initial_vector[1] = -tangent[2]
         initial_vector[2] = tangent[1]
-        frame = spline.moving_frame(t, kind="bishop", initial_vector=initial_vector)
+        frame = spline.moving_frame(t, method="bishop", initial_vector=initial_vector)
 
         # Check that they are unit vectors
         assert np.allclose(np.linalg.norm(frame, axis=-1), 1)
@@ -823,8 +819,74 @@ def test_moving_frame(initialized_spline_curve, not_differentiable_twice):
         # Check that T is parallel to the first derivative
         assert np.allclose(np.cross(frame[:, 0], spline.eval(t, derivative=1)), 0)
 
-        # Torsion tau = dN/ds * B(s) (dot product) should be zero in the Bishop frame
-        # dnds = np.gradient(frame[:, 1], axis=0)
-        # tau = np.sum(dnds * frame[:, 2], axis=-1)
-        # print(tau)
-        # assert np.allclose(tau, 0)
+
+def test_moving_frame_raises():
+    M = 6
+    spline = splinebox.spline_curves.Spline(M=M, basis_function=splinebox.basis_functions.B3(), closed=False)
+    t = np.linspace(0, M - 1, 100)
+
+    spline.knots = np.random.rand(6, 3)
+
+    # Wrong method
+    with pytest.raises(ValueError):
+        spline.moving_frame(t, method="madeup method")
+
+    # Vanishing tangent at the end because of padding
+    with pytest.raises(RuntimeError):
+        spline.moving_frame(t, method="frenet")
+
+    # Straight segment
+    spline.control_points = np.stack([np.arange(M + 2), np.zeros(M + 2), np.zeros(M + 2)], axis=-1)
+    with pytest.raises(RuntimeError):
+        spline.moving_frame(t, method="frenet")
+
+
+def test_mesh(
+    initialized_spline_curve, radius, step_t, step_angle, mesh_type, cap_ends, frame, not_differentiable_twice
+):
+    spline = initialized_spline_curve
+
+    if spline.control_points.ndim != 2 or spline.control_points.shape[1] != 3:
+        # mesh is only implemented for splines in 3D
+        with pytest.raises(RuntimeError):
+            spline.mesh(
+                radius=radius, step_t=step_t, step_angle=step_angle, mesh_type=mesh_type, cap_ends=cap_ends, frame=frame
+            )
+    elif not_differentiable_twice(spline) and radius is not None:
+        with pytest.raises(RuntimeError):
+            spline.mesh(
+                radius=radius, step_t=step_t, step_angle=step_angle, mesh_type=mesh_type, cap_ends=cap_ends, frame=frame
+            )
+    else:
+        points, connectivity = spline.mesh(
+            radius=radius, step_t=step_t, step_angle=step_angle, mesh_type=mesh_type, cap_ends=cap_ends, frame=frame
+        )
+
+        # Check that all points are part of at least one element
+        assert np.all(np.arange(len(points)) == np.sort(np.unique(connectivity)))
+
+        if mesh_type == "surface" and (cap_ends or spline.closed):
+            # Compute the Euler characteristic to check that the mesh is closed
+
+            # Number of verices
+            V = len(points)
+            # Number of faces
+            F = len(connectivity) if connectivity.shape[1] > 2 else 0
+            # Collect all in a set (no duplicates)
+            edges = set()
+            for i in range(len(connectivity)):
+                edges.add(tuple(sorted([connectivity[i, 0], connectivity[i, 1]])))
+                if radius is not None:
+                    edges.add(tuple(sorted([connectivity[i, 0], connectivity[i, 2]])))
+                    edges.add(tuple(sorted([connectivity[i, 1], connectivity[i, 2]])))
+            # Number of edges
+            E = len(edges)
+
+            if radius is None and not spline.closed:
+                assert V + F - E == 1
+            elif spline.closed:
+                # toroidal polyhedra have Euler characteristic 0
+                assert V + F - E == 0
+            elif cap_ends:
+                # spherical polyhedra have Euler characteristic 2
+                assert V + F - E == 2