From 976ffcfb0138278bd61e7828d651a6ef2dd0a4e1 Mon Sep 17 00:00:00 2001 From: faymanns Date: Fri, 13 Dec 2024 11:34:05 +0100 Subject: [PATCH 01/21] First version of surface mesh --- src/splinebox/spline_curves.py | 75 ++++++++++++++++++++++++++++++++++ 1 file changed, 75 insertions(+) diff --git a/src/splinebox/spline_curves.py b/src/splinebox/spline_curves.py index 3634978..2f3175e 100644 --- a/src/splinebox/spline_curves.py +++ b/src/splinebox/spline_curves.py @@ -906,6 +906,81 @@ def _distance(t): return min_distance + def mesh(self, radius=None, resolution=0.1, mesh_type="surface", angular_resolution=10, cap_ends=False): + t = np.arange(0, self.M if self.closed else self.M - 1 + resolution, resolution) + if radius is None: + vertices = self.eval(t) + connections = np.stack((np.arange(len(vertices)), np.arange(len(vertices)) + 1), axis=-1) + if self.closed: + # Connect end to the beginning + connections[-1, -1] = 0 + else: + connections = connections[:-1] + else: + _radius = (lambda t, phi: np.full((len(t),), radius)) if not callable(radius) else radius + + if mesh_type == "surface": + phi = np.arange(0, 360, angular_resolution) + phiphi, tt = np.meshgrid(phi, t) + tt = tt.flatten() + phiphi = phiphi.flatten() + centers = self.eval(tt.flatten()) + rr = _radius(tt, phiphi) + # This should be changed once the no twist framework is implemented + deriv = self.eval(t, derivative=1) + normal1 = np.zeros((len(t), 3)) + normal1[:, 1] = deriv[:, 2] + normal1[:, 2] = -deriv[:, 1] + normal2 = np.zeros((len(t), 3)) + normal2 = np.cross(deriv, normal1) + normal1 /= np.linalg.norm(normal1, axis=1)[:, np.newaxis] + normal2 /= np.linalg.norm(normal2, axis=1)[:, np.newaxis] + + n_angles = len(phi) + n_t = len(t) + normals = ( + np.repeat(normal1, n_angles, axis=0) * np.sin(np.deg2rad(phiphi))[:, np.newaxis] + + np.repeat(normal2, n_angles, axis=0) * np.cos(np.deg2rad(phiphi))[:, np.newaxis] + ) + + vertices = centers + rr[:, np.newaxis] * normals + if self.closed: + connections = np.zeros((2 * n_angles * n_t, 3), dtype=int) + else: + connections = np.zeros((2 * n_angles * (n_t - 1), 3), dtype=int) + face = 0 + n_vertices = len(vertices) + for i in range(n_t if self.closed else n_t - 1): + for j in range(n_angles): + connections[face] = [ + i * n_angles + j, + ((i + 1) * n_angles + j) % n_vertices, + ((i + 1) * n_angles + (j + 1) % n_angles) % n_vertices, + ] + face += 1 + connections[face] = [ + i * n_angles + j, + ((i + 1) * n_angles + (j + 1) % n_angles) % n_vertices, + (i * n_angles + (j + 1) % n_angles) % n_vertices, + ] + face += 1 + if cap_ends and not self.closed: + vertices = np.concatenate( + (centers[0].reshape((1, -1)), vertices, centers[-1].reshape((1, -1))), axis=0 + ) + start_connections = np.zeros((n_angles, 3), dtype=int) + start_connections[:, 1] = np.arange(1, n_angles + 1) + start_connections[:, 2] = np.roll(start_connections[:, 1], -1) + end_connections = np.zeros((n_angles, 3), dtype=int) + end_connections[:, 0] = n_vertices + end_connections[:, 1] = np.arange(n_vertices - 1, n_vertices - 1 - n_angles, -1) + end_connections[:, 2] = np.roll(end_connections[:, 1], -1) + connections = np.concatenate((start_connections, connections + 1, end_connections)) + elif mesh_type == "volume": + pass + + return vertices, connections + class HermiteSpline(Spline): """ From 10932a07f610fd1c659e02bfde59ae9a1dbabbf5 Mon Sep 17 00:00:00 2001 From: faymanns Date: Tue, 17 Dec 2024 12:39:55 +0100 Subject: [PATCH 02/21] Use frame and implement volume mesh --- src/splinebox/spline_curves.py | 72 +++++++++++++++++++++++++++------- 1 file changed, 58 insertions(+), 14 deletions(-) diff --git a/src/splinebox/spline_curves.py b/src/splinebox/spline_curves.py index 2f3175e..dbb988a 100644 --- a/src/splinebox/spline_curves.py +++ b/src/splinebox/spline_curves.py @@ -906,7 +906,16 @@ def _distance(t): return min_distance - def mesh(self, radius=None, resolution=0.1, mesh_type="surface", angular_resolution=10, cap_ends=False): + def mesh( + self, + radius=None, + resolution=0.1, + mesh_type="surface", + angular_resolution=10, + cap_ends=False, + frame="bishop", + initial_vector=None, + ): t = np.arange(0, self.M if self.closed else self.M - 1 + resolution, resolution) if radius is None: vertices = self.eval(t) @@ -919,28 +928,21 @@ def mesh(self, radius=None, resolution=0.1, mesh_type="surface", angular_resolut else: _radius = (lambda t, phi: np.full((len(t),), radius)) if not callable(radius) else radius + phi = np.arange(0, 360, angular_resolution) + if mesh_type == "surface": - phi = np.arange(0, 360, angular_resolution) phiphi, tt = np.meshgrid(phi, t) tt = tt.flatten() phiphi = phiphi.flatten() centers = self.eval(tt.flatten()) rr = _radius(tt, phiphi) - # This should be changed once the no twist framework is implemented - deriv = self.eval(t, derivative=1) - normal1 = np.zeros((len(t), 3)) - normal1[:, 1] = deriv[:, 2] - normal1[:, 2] = -deriv[:, 1] - normal2 = np.zeros((len(t), 3)) - normal2 = np.cross(deriv, normal1) - normal1 /= np.linalg.norm(normal1, axis=1)[:, np.newaxis] - normal2 /= np.linalg.norm(normal2, axis=1)[:, np.newaxis] + normals = self.normal(t, frame=frame, initial_vector=initial_vector) n_angles = len(phi) n_t = len(t) normals = ( - np.repeat(normal1, n_angles, axis=0) * np.sin(np.deg2rad(phiphi))[:, np.newaxis] - + np.repeat(normal2, n_angles, axis=0) * np.cos(np.deg2rad(phiphi))[:, np.newaxis] + 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] ) vertices = centers + rr[:, np.newaxis] * normals @@ -977,7 +979,49 @@ def mesh(self, radius=None, resolution=0.1, mesh_type="surface", angular_resolut end_connections[:, 2] = np.roll(end_connections[:, 1], -1) connections = np.concatenate((start_connections, connections + 1, end_connections)) elif mesh_type == "volume": - pass + 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] + ) + vertices = centers + rr[:, np.newaxis] * normals + if self.closed: + connections = np.zeros((2 * n_angles * n_t, 4), dtype=int) + else: + connections = np.zeros((2 * n_angles * (n_t - 1), 4), dtype=int) + vol = 0 + n_vertices = len(vertices) + for i in range(n_t if self.closed else n_t - 1): + for j in range(1, n_angles + 1): + connections[vol] = [ + i * (n_angles + 1) + j, + ((i + 1) * (n_angles + 1) + j) % n_vertices, + ((i + 1) * (n_angles + 1) + 1 + j % n_angles) % n_vertices, + (i + 1) * (n_angles + 1), + ] + vol += 1 + connections[vol] = [ + i * (n_angles + 1) + j, + ((i + 1) * (n_angles + 1) + 1 + j % n_angles) % n_vertices, + (i * (n_angles + 1) + 1 + j % n_angles) % n_vertices, + i * (n_angles + 1), + ] + vol += 1 return vertices, connections From 2475f2c2603cc543e9568189d3d6700d75b45ba2 Mon Sep 17 00:00:00 2001 From: faymanns Date: Tue, 17 Dec 2024 16:10:25 +0100 Subject: [PATCH 03/21] Add radius=0 as alternative to radius=None --- src/splinebox/spline_curves.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/splinebox/spline_curves.py b/src/splinebox/spline_curves.py index dbb988a..96ee916 100644 --- a/src/splinebox/spline_curves.py +++ b/src/splinebox/spline_curves.py @@ -917,7 +917,7 @@ def mesh( initial_vector=None, ): t = np.arange(0, self.M if self.closed else self.M - 1 + resolution, resolution) - if radius is None: + if radius is None or radius == 0: vertices = self.eval(t) connections = np.stack((np.arange(len(vertices)), np.arange(len(vertices)) + 1), axis=-1) if self.closed: From 9bfa83d2cd487f172c0400f3dd2536e9fa27aa1a Mon Sep 17 00:00:00 2001 From: faymanns Date: Tue, 17 Dec 2024 16:11:10 +0100 Subject: [PATCH 04/21] Add example with meshes --- examples/plot_mesh.py | 162 ++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 162 insertions(+) create mode 100644 examples/plot_mesh.py diff --git a/examples/plot_mesh.py b/examples/plot_mesh.py new file mode 100644 index 0000000..8b12e19 --- /dev/null +++ b/examples/plot_mesh.py @@ -0,0 +1,162 @@ +""" +Spline to mesh +============== +This example demonstrates how you can turn a spline curve in 3D into a mesh. +""" + +# sphinx_gallery_thumbnail_number = 4 + +import numpy as np +import pyvista +import splinebox + +# %% +# 1. Construct a spline +# --------------------- +# We begin by constructing 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. Generate a mesh with no radius +# --------------------------------- +# The resolution determines how fine the resulting mesh is an corresponds +# to the step size in the spline parameter space t. +# Setting the radius to :code:`None` or 0 results in points along the spline connected by lines. +vertices, connections = spline.mesh(resolution=0.1, radius=None) + +# %% +# To visulainze the "mesh" with pyvista we have to preprend the number of vertices +# that are connected by a specific connection. In this case the number of vertices +# is two since we are simply connecing a two points with a line. +# Lastly, we have to turn the vertices and connections into a pyvista mesh using the +# :code:`PolyData` class. +connections = np.hstack((np.full((connections.shape[0], 1), 2), connections)) +mesh = pyvista.PolyData(vertices, lines=connections) +mesh.plot() + +# %% +# 3. Fixed radius +# --------------- +# Next, we will use a fixed radius to generate a surface mesh ("tube"). +# We choose the Frenet-Serret frame in order to avoid having to choose an initial vector. +vertices, connections = spline.mesh(resolution=0.1, radius=0.2, frame="frenet") +# %% +# Once again, we have to prepend the number how vertices in each connection. +# Since, our mesh consists of triangles we choose 3. +connections = np.hstack((np.full((connections.shape[0], 1), 3), connections)) +mesh = pyvista.PolyData(vertices, faces=connections) +mesh.plot(show_edges=True) + +# %% +# 3. Eliptical cross section +# -------------------------- +# In order to create, a spline with a custom cross section we can specify +# the radius as a function. The function was to take the spline parameter :code:`t` +# and the angle polar angle :code:`phi` as arguments in that order. +# For this example we create an 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 + + +vertices, connections = spline.mesh(resolution=0.1, angular_resolution=36, radius=elliptical_radius, frame="frenet") +connections = np.hstack((np.full((connections.shape[0], 1), 3), connections)) +mesh = pyvista.PolyData(vertices, faces=connections) +mesh.plot(show_edges=True) + +# %% +# Or you can change the radius along the spline. + + +def radius(t, phi): + return 0.1 + 0.03 * np.sin(t / spline.M * 16 * np.pi) + + +vertices, connections = spline.mesh(resolution=0.1, angular_resolution=36, radius=radius, frame="frenet") +connections = np.hstack((np.full((connections.shape[0], 1), 3), connections)) +mesh = pyvista.PolyData(vertices, faces=connections) +mesh.plot(show_edges=True) + +# %% +# 4. Bishop frame +# --------------- +# 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. +# Let's create a spline with some unwated twist and see how the Bishop frame fixes it. +spline = splinebox.Spline(M=M, 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], + ] +) + +vertices, connections = spline.mesh(resolution=0.1, angular_resolution=36, radius=elliptical_radius, frame="frenet") +connections = np.hstack((np.full((connections.shape[0], 1), 3), connections)) +mesh = pyvista.PolyData(vertices, faces=connections) +mesh.plot(show_edges=True) + +# %% +# You can clearly see how the ellipse twists around the spline. +# The Bishop frame eliminates this twist but requires an initial orientation for the ellipse, given by an initial vector. + +initial_vector = np.zeros(3) +tangent = spline.eval(0, derivative=1) +initial_vector[1] = tangent[2] +initial_vector[2] = -tangent[1] + +vertices, connections = spline.mesh( + resolution=0.1, angular_resolution=36, radius=elliptical_radius, frame="bishop", initial_vector=initial_vector +) +connections = np.hstack((np.full((connections.shape[0], 1), 3), connections)) +mesh = pyvista.PolyData(vertices, faces=connections) +mesh.plot(show_edges=True) + +# %% +# Changing the initial vector rotates the ellipse, because the initial vector is the reference for phi=0. + +initial_vector = np.array([0.5, -0.5, 1]) + +vertices, connections = spline.mesh( + resolution=0.1, angular_resolution=36, radius=elliptical_radius, frame="bishop", initial_vector=initial_vector +) +connections = np.hstack((np.full((connections.shape[0], 1), 3), connections)) +mesh = pyvista.PolyData(vertices, faces=connections) +mesh.plot(show_edges=True) + +# %% +# 5. Volume mesh +# -------------- +# Besides a surface mesh, we can also turn the spline into a volume mesh. + +vertices, connections = spline.mesh( + radius=radius, resolution=0.5, angular_resolution=72, initial_vector=initial_vector, mesh_type="volume" +) +connections = np.hstack((np.full((connections.shape[0], 1), 4), connections)) +cell_types = np.full(len(connections), fill_value=pyvista.CellType.TETRA, dtype=np.uint8) + +mesh = pyvista.UnstructuredGrid(connections, cell_types, vertices) +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 +# ---- +# * To save the meshes for visulalisation in ParaView, you can use pyvista: :code:`mesh.save("mesh.vtk")` +# * By default the surface meshes are open at the end. You can close them using the :code:`cap_ends` keyword argument of :meth:`splinebox.spline_curves.Spline.mesh`. From cf8ad175b00d3fa5be61ff194f6c224de16118e0 Mon Sep 17 00:00:00 2001 From: faymanns Date: Tue, 17 Dec 2024 17:08:15 +0100 Subject: [PATCH 05/21] Example showing the different frames --- examples/plot_frame.py | 87 ++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 87 insertions(+) create mode 100644 examples/plot_frame.py diff --git a/examples/plot_frame.py b/examples/plot_frame.py new file mode 100644 index 0000000..892f278 --- /dev/null +++ b/examples/plot_frame.py @@ -0,0 +1,87 @@ +""" +Moving Frames +============= + +`Moving frames`_ can be used to to implment a local coordinate system on a curve. +This coordinate system can them be used to interogate the space around a point on the spline, +for instance orthogonal image can be extracted from a volume as in the :ref:`Dendrite-Centric Coordinate System`. +Splinebox implements two different frames, the `Frenet-Serret`_ and the `Bishop`_ [Bishop1975]_ frame. +The Bishop frame has the advantage that it does not twist around the curve. + +.. _Moving frames: https://en.wikipedia.org/wiki/Moving_frame +.. _Frenet-Serret: https://en.wikipedia.org/wiki/Frenet-Serret_formulas +.. _Bishop: https://www.jstor.org/stable/2319846 +""" + +import numpy as np +import pyvista +import splinebox + +# %% +# Let's construct a spline and visualize the two frames. + +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) + +# %% +# First, we will compute the Frenet-Serret frame. +frenet_frame = spline.moving_frame(t, kind="frenet") +print(frenet_frame.shape) + +# %% +# The returned array contains three vectors that form an orthonormal basis +# for each t. +# To compute the Bishop frame we need an initial vector that fixes the roation +# around the curve. This vector has to be orthogonal to the tangent at :code:`t[0]`. + +tangent = spline.eval(t[0], derivative=1) +initial_vector = np.zeros(3) +initial_vector[1] = -tangent[2] +initial_vector[2] = tangent[1] + +bishop_frame = spline.moving_frame(t, kind="bishop", initial_vector=initial_vector) + +# %% +# To visualize the frames we can use a pyvista mesh. +# In the plot below it is clearly visible how the Frenet-Serret frame on the left +# twists around the curve. + +spline_mesh = pyvista.MultipleLines(points=spline.eval(t)) + +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 + +plotter = pyvista.Plotter(shape=(1, 2), border=False) +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") +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() From f8433c18228b92ac5cf8534e812f532ee89e1ebf Mon Sep 17 00:00:00 2001 From: faymanns Date: Thu, 19 Dec 2024 11:27:24 +0100 Subject: [PATCH 06/21] Fix end cap of surface mesh --- src/splinebox/spline_curves.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/splinebox/spline_curves.py b/src/splinebox/spline_curves.py index 96ee916..3234f4c 100644 --- a/src/splinebox/spline_curves.py +++ b/src/splinebox/spline_curves.py @@ -974,8 +974,8 @@ def mesh( start_connections[:, 1] = np.arange(1, n_angles + 1) start_connections[:, 2] = np.roll(start_connections[:, 1], -1) end_connections = np.zeros((n_angles, 3), dtype=int) - end_connections[:, 0] = n_vertices - end_connections[:, 1] = np.arange(n_vertices - 1, n_vertices - 1 - n_angles, -1) + end_connections[:, 0] = n_vertices + 1 + end_connections[:, 1] = np.arange(n_vertices, n_vertices - n_angles, -1) end_connections[:, 2] = np.roll(end_connections[:, 1], -1) connections = np.concatenate((start_connections, connections + 1, end_connections)) elif mesh_type == "volume": From bc2b52028a3cf988310c529b6d9fc2f19a4139dd Mon Sep 17 00:00:00 2001 From: faymanns Date: Thu, 19 Dec 2024 11:45:03 +0100 Subject: [PATCH 07/21] Fix lambda function for radius The hstack in the volume mesh failed when the radius was not callable. --- src/splinebox/spline_curves.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/splinebox/spline_curves.py b/src/splinebox/spline_curves.py index 3234f4c..5538265 100644 --- a/src/splinebox/spline_curves.py +++ b/src/splinebox/spline_curves.py @@ -926,7 +926,7 @@ def mesh( else: connections = connections[:-1] else: - _radius = (lambda t, phi: np.full((len(t),), radius)) if not callable(radius) else radius + _radius = (lambda t, phi: np.full(t.shape, radius)) if not callable(radius) else radius phi = np.arange(0, 360, angular_resolution) From 22b4ec9872e12e8838c9e3aa343c214f4615c3b9 Mon Sep 17 00:00:00 2001 From: faymanns Date: Thu, 19 Dec 2024 11:49:14 +0100 Subject: [PATCH 08/21] Rename vertices to points and connections to connectivity --- examples/plot_mesh.py | 58 ++++++++++++++-------------- src/splinebox/spline_curves.py | 70 +++++++++++++++++----------------- 2 files changed, 63 insertions(+), 65 deletions(-) diff --git a/examples/plot_mesh.py b/examples/plot_mesh.py index 8b12e19..8846bcc 100644 --- a/examples/plot_mesh.py +++ b/examples/plot_mesh.py @@ -24,16 +24,16 @@ # The resolution determines how fine the resulting mesh is an corresponds # to the step size in the spline parameter space t. # Setting the radius to :code:`None` or 0 results in points along the spline connected by lines. -vertices, connections = spline.mesh(resolution=0.1, radius=None) +points, connectivity = spline.mesh(resolution=0.1, radius=None) # %% -# To visulainze the "mesh" with pyvista we have to preprend the number of vertices -# that are connected by a specific connection. In this case the number of vertices +# To visulainze the "mesh" with pyvista we have to preprend the number of points +# that are connected by a specific connection. In this case the number of points # is two since we are simply connecing a two points with a line. -# Lastly, we have to turn the vertices and connections into a pyvista mesh using the +# Lastly, we have to turn the points and connectivity into a pyvista mesh using the # :code:`PolyData` class. -connections = np.hstack((np.full((connections.shape[0], 1), 2), connections)) -mesh = pyvista.PolyData(vertices, lines=connections) +connectivity = np.hstack((np.full((connectivity.shape[0], 1), 2), connectivity)) +mesh = pyvista.PolyData(points, lines=connectivity) mesh.plot() # %% @@ -41,12 +41,12 @@ # --------------- # Next, we will use a fixed radius to generate a surface mesh ("tube"). # We choose the Frenet-Serret frame in order to avoid having to choose an initial vector. -vertices, connections = spline.mesh(resolution=0.1, radius=0.2, frame="frenet") +points, connectivity = spline.mesh(resolution=0.1, radius=0.2, frame="frenet") # %% -# Once again, we have to prepend the number how vertices in each connection. +# Once again, we have to prepend the number how points in each connection. # Since, our mesh consists of triangles we choose 3. -connections = np.hstack((np.full((connections.shape[0], 1), 3), connections)) -mesh = pyvista.PolyData(vertices, faces=connections) +connectivity = np.hstack((np.full((connectivity.shape[0], 1), 3), connectivity)) +mesh = pyvista.PolyData(points, faces=connectivity) mesh.plot(show_edges=True) # %% @@ -66,9 +66,9 @@ def elliptical_radius(t, phi): return r -vertices, connections = spline.mesh(resolution=0.1, angular_resolution=36, radius=elliptical_radius, frame="frenet") -connections = np.hstack((np.full((connections.shape[0], 1), 3), connections)) -mesh = pyvista.PolyData(vertices, faces=connections) +points, connectivity = spline.mesh(resolution=0.1, angular_resolution=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) # %% @@ -79,9 +79,9 @@ def radius(t, phi): return 0.1 + 0.03 * np.sin(t / spline.M * 16 * np.pi) -vertices, connections = spline.mesh(resolution=0.1, angular_resolution=36, radius=radius, frame="frenet") -connections = np.hstack((np.full((connections.shape[0], 1), 3), connections)) -mesh = pyvista.PolyData(vertices, faces=connections) +points, connectivity = spline.mesh(resolution=0.1, angular_resolution=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) # %% @@ -103,9 +103,9 @@ def radius(t, phi): ] ) -vertices, connections = spline.mesh(resolution=0.1, angular_resolution=36, radius=elliptical_radius, frame="frenet") -connections = np.hstack((np.full((connections.shape[0], 1), 3), connections)) -mesh = pyvista.PolyData(vertices, faces=connections) +points, connectivity = spline.mesh(resolution=0.1, angular_resolution=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) # %% @@ -117,11 +117,11 @@ def radius(t, phi): initial_vector[1] = tangent[2] initial_vector[2] = -tangent[1] -vertices, connections = spline.mesh( +points, connectivity = spline.mesh( resolution=0.1, angular_resolution=36, radius=elliptical_radius, frame="bishop", initial_vector=initial_vector ) -connections = np.hstack((np.full((connections.shape[0], 1), 3), connections)) -mesh = pyvista.PolyData(vertices, faces=connections) +connectivity = np.hstack((np.full((connectivity.shape[0], 1), 3), connectivity)) +mesh = pyvista.PolyData(points, faces=connectivity) mesh.plot(show_edges=True) # %% @@ -129,11 +129,11 @@ def radius(t, phi): initial_vector = np.array([0.5, -0.5, 1]) -vertices, connections = spline.mesh( +points, connectivity = spline.mesh( resolution=0.1, angular_resolution=36, radius=elliptical_radius, frame="bishop", initial_vector=initial_vector ) -connections = np.hstack((np.full((connections.shape[0], 1), 3), connections)) -mesh = pyvista.PolyData(vertices, faces=connections) +connectivity = np.hstack((np.full((connectivity.shape[0], 1), 3), connectivity)) +mesh = pyvista.PolyData(points, faces=connectivity) mesh.plot(show_edges=True) # %% @@ -141,13 +141,13 @@ def radius(t, phi): # -------------- # Besides a surface mesh, we can also turn the spline into a volume mesh. -vertices, connections = spline.mesh( +points, connectivity = spline.mesh( radius=radius, resolution=0.5, angular_resolution=72, initial_vector=initial_vector, mesh_type="volume" ) -connections = np.hstack((np.full((connections.shape[0], 1), 4), connections)) -cell_types = np.full(len(connections), fill_value=pyvista.CellType.TETRA, dtype=np.uint8) +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(connections, cell_types, vertices) +mesh = pyvista.UnstructuredGrid(connectivity, cell_types, points) mesh.plot(show_edges=True) # %% diff --git a/src/splinebox/spline_curves.py b/src/splinebox/spline_curves.py index 5538265..ce994ec 100644 --- a/src/splinebox/spline_curves.py +++ b/src/splinebox/spline_curves.py @@ -918,13 +918,13 @@ def mesh( ): t = np.arange(0, self.M if self.closed else self.M - 1 + resolution, resolution) if radius is None or radius == 0: - vertices = self.eval(t) - connections = np.stack((np.arange(len(vertices)), np.arange(len(vertices)) + 1), axis=-1) + 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 - connections[-1, -1] = 0 + connectivity[-1, -1] = 0 else: - connections = connections[:-1] + connectivity = connectivity[:-1] else: _radius = (lambda t, phi: np.full(t.shape, radius)) if not callable(radius) else radius @@ -945,39 +945,37 @@ def mesh( + np.repeat(normals[:, 1], n_angles, axis=0) * np.cos(np.deg2rad(phiphi))[:, np.newaxis] ) - vertices = centers + rr[:, np.newaxis] * normals + points = centers + rr[:, np.newaxis] * normals if self.closed: - connections = np.zeros((2 * n_angles * n_t, 3), dtype=int) + connectivity = np.zeros((2 * n_angles * n_t, 3), dtype=int) else: - connections = np.zeros((2 * n_angles * (n_t - 1), 3), dtype=int) + connectivity = np.zeros((2 * n_angles * (n_t - 1), 3), dtype=int) face = 0 - n_vertices = len(vertices) + n_points = len(points) for i in range(n_t if self.closed else n_t - 1): for j in range(n_angles): - connections[face] = [ + connectivity[face] = [ i * n_angles + j, - ((i + 1) * n_angles + j) % n_vertices, - ((i + 1) * n_angles + (j + 1) % n_angles) % n_vertices, + ((i + 1) * n_angles + j) % n_points, + ((i + 1) * n_angles + (j + 1) % n_angles) % n_points, ] face += 1 - connections[face] = [ + connectivity[face] = [ i * n_angles + j, - ((i + 1) * n_angles + (j + 1) % n_angles) % n_vertices, - (i * n_angles + (j + 1) % n_angles) % n_vertices, + ((i + 1) * n_angles + (j + 1) % n_angles) % n_points, + (i * n_angles + (j + 1) % n_angles) % n_points, ] face += 1 if cap_ends and not self.closed: - vertices = np.concatenate( - (centers[0].reshape((1, -1)), vertices, centers[-1].reshape((1, -1))), axis=0 - ) - start_connections = np.zeros((n_angles, 3), dtype=int) - start_connections[:, 1] = np.arange(1, n_angles + 1) - start_connections[:, 2] = np.roll(start_connections[:, 1], -1) - end_connections = np.zeros((n_angles, 3), dtype=int) - end_connections[:, 0] = n_vertices + 1 - end_connections[:, 1] = np.arange(n_vertices, n_vertices - n_angles, -1) - end_connections[:, 2] = np.roll(end_connections[:, 1], -1) - connections = np.concatenate((start_connections, connections + 1, end_connections)) + 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) @@ -999,31 +997,31 @@ def mesh( 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] ) - vertices = centers + rr[:, np.newaxis] * normals + points = centers + rr[:, np.newaxis] * normals if self.closed: - connections = np.zeros((2 * n_angles * n_t, 4), dtype=int) + connectivity = np.zeros((2 * n_angles * n_t, 4), dtype=int) else: - connections = np.zeros((2 * n_angles * (n_t - 1), 4), dtype=int) + connectivity = np.zeros((2 * n_angles * (n_t - 1), 4), dtype=int) vol = 0 - n_vertices = len(vertices) + n_points = len(points) for i in range(n_t if self.closed else n_t - 1): for j in range(1, n_angles + 1): - connections[vol] = [ + connectivity[vol] = [ i * (n_angles + 1) + j, - ((i + 1) * (n_angles + 1) + j) % n_vertices, - ((i + 1) * (n_angles + 1) + 1 + j % n_angles) % n_vertices, + ((i + 1) * (n_angles + 1) + j) % n_points, + ((i + 1) * (n_angles + 1) + 1 + j % n_angles) % n_points, (i + 1) * (n_angles + 1), ] vol += 1 - connections[vol] = [ + connectivity[vol] = [ i * (n_angles + 1) + j, - ((i + 1) * (n_angles + 1) + 1 + j % n_angles) % n_vertices, - (i * (n_angles + 1) + 1 + j % n_angles) % n_vertices, + ((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 - return vertices, connections + return points, connectivity class HermiteSpline(Spline): From 8b7105fce9b93f357b23cb2b04f6f62e2a049238 Mon Sep 17 00:00:00 2001 From: faymanns Date: Thu, 19 Dec 2024 13:14:26 +0100 Subject: [PATCH 09/21] Add missing tetrahedron --- src/splinebox/spline_curves.py | 11 +++++++++-- 1 file changed, 9 insertions(+), 2 deletions(-) diff --git a/src/splinebox/spline_curves.py b/src/splinebox/spline_curves.py index ce994ec..f5a11e2 100644 --- a/src/splinebox/spline_curves.py +++ b/src/splinebox/spline_curves.py @@ -999,9 +999,9 @@ def mesh( ) points = centers + rr[:, np.newaxis] * normals if self.closed: - connectivity = np.zeros((2 * n_angles * n_t, 4), dtype=int) + connectivity = np.zeros((3 * n_angles * n_t, 4), dtype=int) else: - connectivity = np.zeros((2 * n_angles * (n_t - 1), 4), dtype=int) + connectivity = np.zeros((3 * n_angles * (n_t - 1), 4), dtype=int) vol = 0 n_points = len(points) for i in range(n_t if self.closed else n_t - 1): @@ -1020,6 +1020,13 @@ def mesh( 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), + i * (n_angles + 1), + ] + vol += 1 return points, connectivity From dac4e4b4215a119c5de5f7916d2320b71cd988ec Mon Sep 17 00:00:00 2001 From: faymanns Date: Thu, 19 Dec 2024 14:14:23 +0100 Subject: [PATCH 10/21] Add docstring to mesh --- src/splinebox/spline_curves.py | 46 ++++++++++++++++++++++++++++++---- 1 file changed, 41 insertions(+), 5 deletions(-) diff --git a/src/splinebox/spline_curves.py b/src/splinebox/spline_curves.py index f5a11e2..d1f0129 100644 --- a/src/splinebox/spline_curves.py +++ b/src/splinebox/spline_curves.py @@ -909,14 +909,50 @@ def _distance(t): def mesh( self, radius=None, - resolution=0.1, + step_t=0.1, + step_angle="auto", mesh_type="surface", - angular_resolution=10, cap_ends=False, frame="bishop", - initial_vector=None, + initial_vector="auto", ): - t = np.arange(0, self.M if self.closed else self.M - 1 + resolution, resolution) + """ + Creats a mesh around the spline curve in 3D. The distance of the mesh + from the spline can be changes using the :code:`radius` paramters. + The mesh can either be a surface with triangular cells or a volume mesh + with tetrahedral cells. + + Parameters + ---------- + radius : float or callable + A float value results in a constant radius. Alternatively, it can be a function that takes + the spline parameter t and the polar angle in the normal plane as arguments and returns + a float. + step_t : float or "auto" + The step size of the spline parameter t. If "auto", step_angle has to be a float. + step_angle : float or "auto" + The step size of the polar angle. If "auto", step_t has to be a float. + mesh_type : str + Can be "surface" or "volume". + cap_ends : boolean + If True, the ends of the surface mesh are closed with orthogonal planes. + frame : str + Can be "frenet" or "bishop". See :meth:`splinebox.spline_curves.moving_frame`. + initial_vector : numpy array or "auto" + The initial vector determining the orientation of the Bishop frame. + See :meth:`splinebox.spline_curves.moving_frame`. + + Returns + ------- + points : numpy array + A 2D numpy array. The second dimension encodes the coordinates of each point in 3D space. + connectivity : numpy array + A 2D numpy array. Each row corresponds to an element in the mesh. + The values are the indices of the points in :code:`points` that form the element. + """ + 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) @@ -928,7 +964,7 @@ def mesh( else: _radius = (lambda t, phi: np.full(t.shape, radius)) if not callable(radius) else radius - phi = np.arange(0, 360, angular_resolution) + phi = np.arange(0, 360, step_angle) if mesh_type == "surface": phiphi, tt = np.meshgrid(phi, t) From c298d62930a26aa8e6316b44a7e212ffab5a2b76 Mon Sep 17 00:00:00 2001 From: faymanns Date: Thu, 19 Dec 2024 14:29:11 +0100 Subject: [PATCH 11/21] Calculate initial vector for bishop frame automatically --- src/splinebox/spline_curves.py | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/src/splinebox/spline_curves.py b/src/splinebox/spline_curves.py index d1f0129..101e26b 100644 --- a/src/splinebox/spline_curves.py +++ b/src/splinebox/spline_curves.py @@ -668,7 +668,7 @@ def moving_frame(self, t, kind="frenet", initial_vector=None): 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. + the curve. If None, an initial vector is computed automatically. Returns ------- @@ -703,7 +703,11 @@ def moving_frame(self, t, kind="frenet", initial_vector=None): frame[:, 1] = np.cross(frame[:, 2], frame[:, 0]) elif kind == "bishop": if initial_vector is None: - raise ValueError("The bishop frame requires the keyword argument initial_vector.") + initial_vector = np.zeros(3) + tangent = frame[0, 0] + max_axis = np.argmax(tangent) + initial_vector[max_axis] = tangent[(max_axis + 1) % 1] + initial_vector[(max_axis + 1) % 1] = -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].") @@ -914,7 +918,7 @@ def mesh( mesh_type="surface", cap_ends=False, frame="bishop", - initial_vector="auto", + initial_vector=None, ): """ Creats a mesh around the spline curve in 3D. The distance of the mesh @@ -938,7 +942,7 @@ def mesh( If True, the ends of the surface mesh are closed with orthogonal planes. frame : str Can be "frenet" or "bishop". See :meth:`splinebox.spline_curves.moving_frame`. - initial_vector : numpy array or "auto" + initial_vector : numpy array or None The initial vector determining the orientation of the Bishop frame. See :meth:`splinebox.spline_curves.moving_frame`. From c96feaf750d55c8dfad43691aafdb5220cae6469 Mon Sep 17 00:00:00 2001 From: faymanns Date: Thu, 19 Dec 2024 14:49:12 +0100 Subject: [PATCH 12/21] Make default initial vector of Bishop consistent with Frenet whenever possible --- examples/plot_frame.py | 1 + src/splinebox/spline_curves.py | 12 ++++++++---- 2 files changed, 9 insertions(+), 4 deletions(-) diff --git a/examples/plot_frame.py b/examples/plot_frame.py index 892f278..38305e2 100644 --- a/examples/plot_frame.py +++ b/examples/plot_frame.py @@ -49,6 +49,7 @@ initial_vector = np.zeros(3) initial_vector[1] = -tangent[2] initial_vector[2] = tangent[1] +initial_vector = None bishop_frame = spline.moving_frame(t, kind="bishop", initial_vector=initial_vector) diff --git a/src/splinebox/spline_curves.py b/src/splinebox/spline_curves.py index 101e26b..08055df 100644 --- a/src/splinebox/spline_curves.py +++ b/src/splinebox/spline_curves.py @@ -703,11 +703,15 @@ def moving_frame(self, t, kind="frenet", initial_vector=None): frame[:, 1] = np.cross(frame[:, 2], frame[:, 0]) elif kind == "bishop": if initial_vector is None: - initial_vector = np.zeros(3) tangent = frame[0, 0] - max_axis = np.argmax(tangent) - initial_vector[max_axis] = tangent[(max_axis + 1) % 1] - initial_vector[(max_axis + 1) % 1] = -tangent[max_axis] + # 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].") From bc0c6d73973a30d4153b6ce0dd08b379646b4cbe Mon Sep 17 00:00:00 2001 From: faymanns Date: Thu, 19 Dec 2024 15:15:09 +0100 Subject: [PATCH 13/21] Remove auto for step sizes --- src/splinebox/spline_curves.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/splinebox/spline_curves.py b/src/splinebox/spline_curves.py index 08055df..e0f47e1 100644 --- a/src/splinebox/spline_curves.py +++ b/src/splinebox/spline_curves.py @@ -918,7 +918,7 @@ def mesh( self, radius=None, step_t=0.1, - step_angle="auto", + step_angle=5, mesh_type="surface", cap_ends=False, frame="bishop", @@ -936,10 +936,10 @@ def mesh( A float value results in a constant radius. Alternatively, it can be a function that takes the spline parameter t and the polar angle in the normal plane as arguments and returns a float. - step_t : float or "auto" - The step size of the spline parameter t. If "auto", step_angle has to be a float. - step_angle : float or "auto" - The step size of the polar angle. If "auto", step_t has to be a float. + step_t : float + The step size of the spline parameter t. + step_angle : float + The step size of the polar angle. mesh_type : str Can be "surface" or "volume". cap_ends : boolean From 7b87bb7101d5e35c4d8f79e407c573487a9004a1 Mon Sep 17 00:00:00 2001 From: faymanns Date: Thu, 19 Dec 2024 18:19:10 +0100 Subject: [PATCH 14/21] Test mesh --- tests/conftest.py | 38 +++++++++++++++++++++++++ tests/test_spline_curves.py | 55 +++++++++++++++++++++++++++++++++---- 2 files changed, 88 insertions(+), 5 deletions(-) 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..8173da6 100644 --- a/tests/test_spline_curves.py +++ b/tests/test_spline_curves.py @@ -823,8 +823,53 @@ 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_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 From cac663c7d342ba7bb610c32d7993f6811c1b8fea Mon Sep 17 00:00:00 2001 From: faymanns Date: Thu, 19 Dec 2024 18:20:12 +0100 Subject: [PATCH 15/21] Fix overflow of indices for closed volume meshes --- src/splinebox/spline_curves.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/splinebox/spline_curves.py b/src/splinebox/spline_curves.py index e0f47e1..acf865d 100644 --- a/src/splinebox/spline_curves.py +++ b/src/splinebox/spline_curves.py @@ -1054,7 +1054,7 @@ def mesh( 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), + (i + 1) * (n_angles + 1) % n_points, ] vol += 1 connectivity[vol] = [ @@ -1067,7 +1067,7 @@ def mesh( connectivity[vol] = [ i * (n_angles + 1) + j, ((i + 1) * (n_angles + 1) + 1 + j % n_angles) % n_points, - (i + 1) * (n_angles + 1), + (i + 1) * (n_angles + 1) % n_points, i * (n_angles + 1), ] vol += 1 From 9527c19ab098e59eb3d79b71514558aa8dcc2c37 Mon Sep 17 00:00:00 2001 From: faymanns Date: Fri, 20 Dec 2024 13:05:45 +0100 Subject: [PATCH 16/21] Remove test for missing initial vector since I implemented an automatic initial vector calculation --- tests/test_spline_curves.py | 4 ---- 1 file changed, 4 deletions(-) diff --git a/tests/test_spline_curves.py b/tests/test_spline_curves.py index 8173da6..1f576a9 100644 --- a/tests/test_spline_curves.py +++ b/tests/test_spline_curves.py @@ -798,10 +798,6 @@ 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)) From 6659b579846e1ece813e10b782d1a01ba85f4694 Mon Sep 17 00:00:00 2001 From: faymanns Date: Fri, 20 Dec 2024 13:06:30 +0100 Subject: [PATCH 17/21] Use numba.jit for mesh connectivity --- src/splinebox/spline_curves.py | 108 +++++++++++++++++++-------------- 1 file changed, 61 insertions(+), 47 deletions(-) diff --git a/src/splinebox/spline_curves.py b/src/splinebox/spline_curves.py index acf865d..a5b9595 100644 --- a/src/splinebox/spline_curves.py +++ b/src/splinebox/spline_curves.py @@ -990,26 +990,8 @@ def mesh( ) points = centers + rr[:, np.newaxis] * normals - if self.closed: - connectivity = np.zeros((2 * n_angles * n_t, 3), dtype=int) - else: - connectivity = np.zeros((2 * n_angles * (n_t - 1), 3), dtype=int) - face = 0 n_points = len(points) - for i in range(n_t if self.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 + 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) @@ -1020,6 +1002,7 @@ def mesh( 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) @@ -1042,38 +1025,69 @@ def mesh( + np.repeat(normals[:, 1], n_angles + 1, axis=0) * np.cos(np.deg2rad(phiphi))[:, np.newaxis] ) points = centers + rr[:, np.newaxis] * normals - if self.closed: - connectivity = np.zeros((3 * n_angles * n_t, 4), dtype=int) - else: - connectivity = np.zeros((3 * n_angles * (n_t - 1), 4), dtype=int) - vol = 0 n_points = len(points) - for i in range(n_t if self.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 + + 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): + 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): + 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): """ From e0aaf95ccdee1d911fdb019553bb541ea77c6b75 Mon Sep 17 00:00:00 2001 From: faymanns Date: Fri, 20 Dec 2024 16:57:30 +0100 Subject: [PATCH 18/21] Clean up examples --- examples/plot_frame.py | 57 +++++++++++-------- examples/plot_mesh.py | 123 ++++++++++++++++++++++------------------- 2 files changed, 99 insertions(+), 81 deletions(-) diff --git a/examples/plot_frame.py b/examples/plot_frame.py index 38305e2..0cbef22 100644 --- a/examples/plot_frame.py +++ b/examples/plot_frame.py @@ -2,15 +2,16 @@ Moving Frames ============= -`Moving frames`_ can be used to to implment a local coordinate system on a curve. -This coordinate system can them be used to interogate the space around a point on the spline, -for instance orthogonal image can be extracted from a volume as in the :ref:`Dendrite-Centric Coordinate System`. -Splinebox implements two different frames, the `Frenet-Serret`_ and the `Bishop`_ [Bishop1975]_ frame. -The Bishop frame has the advantage that it does not twist around the curve. +`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: https://en.wikipedia.org/wiki/Frenet-Serret_formulas -.. _Bishop: https://www.jstor.org/stable/2319846 +.. _Frenet-Serret Frame: https://en.wikipedia.org/wiki/Frenet-Serret_formulas +.. _Bishop Frame: https://www.jstor.org/stable/2319846 """ import numpy as np @@ -18,7 +19,8 @@ import splinebox # %% -# Let's construct a spline and visualize the two frames. +# 1. Create a Spline +# ^^^^^^^^^^^^^^^^^^ spline = splinebox.Spline(M=4, basis_function=splinebox.B3(), closed=False) spline.control_points = np.array( @@ -35,31 +37,28 @@ t = np.linspace(0, spline.M - 1, spline.M * 3) # %% -# First, we will compute the Frenet-Serret frame. +# 2. Compute the Frenet-Serret Frame +# ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + frenet_frame = spline.moving_frame(t, kind="frenet") -print(frenet_frame.shape) # %% -# The returned array contains three vectors that form an orthonormal basis -# for each t. -# To compute the Bishop frame we need an initial vector that fixes the roation -# around the curve. This vector has to be orthogonal to the tangent at :code:`t[0]`. +# 3. Compute the Bishop Frame +# ^^^^^^^^^^^^^^^^^^^^^^^^^^^ -tangent = spline.eval(t[0], derivative=1) -initial_vector = np.zeros(3) -initial_vector[1] = -tangent[2] -initial_vector[2] = tangent[1] -initial_vector = None - -bishop_frame = spline.moving_frame(t, kind="bishop", initial_vector=initial_vector) +bishop_frame = spline.moving_frame(t, kind="bishop") # %% -# To visualize the frames we can use a pyvista mesh. -# In the plot below it is clearly visible how the Frenet-Serret frame on the left -# twists around the curve. +# 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 @@ -67,7 +66,10 @@ 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") @@ -76,6 +78,8 @@ 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") @@ -84,5 +88,10 @@ 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 index 8846bcc..c19c115 100644 --- a/examples/plot_mesh.py +++ b/examples/plot_mesh.py @@ -1,7 +1,7 @@ """ Spline to mesh ============== -This example demonstrates how you can turn a spline curve in 3D into a mesh. +This guide demonstrates how to convert a 3D spline curve into various types of meshes using SplineBox. """ # sphinx_gallery_thumbnail_number = 4 @@ -11,51 +11,51 @@ import splinebox # %% -# 1. Construct a spline -# --------------------- -# We begin by constructing a circular spline. +# 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. Generate a mesh with no radius -# --------------------------------- -# The resolution determines how fine the resulting mesh is an corresponds -# to the step size in the spline parameter space t. -# Setting the radius to :code:`None` or 0 results in points along the spline connected by lines. -points, connectivity = spline.mesh(resolution=0.1, radius=None) +# 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. -# %% -# To visulainze the "mesh" with pyvista we have to preprend the number of points -# that are connected by a specific connection. In this case the number of points -# is two since we are simply connecing a two points with a line. -# Lastly, we have to turn the points and connectivity into a pyvista mesh using the -# :code:`PolyData` class. +# 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. Fixed radius -# --------------- -# Next, we will use a fixed radius to generate a surface mesh ("tube"). -# We choose the Frenet-Serret frame in order to avoid having to choose an initial vector. -points, connectivity = spline.mesh(resolution=0.1, radius=0.2, frame="frenet") -# %% -# Once again, we have to prepend the number how points in each connection. -# Since, our mesh consists of triangles we choose 3. +# 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. Eliptical cross section -# -------------------------- -# In order to create, a spline with a custom cross section we can specify -# the radius as a function. The function was to take the spline parameter :code:`t` -# and the angle polar angle :code:`phi` as arguments in that order. -# For this example we create an elliptical cross section. +# 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): @@ -66,84 +66,93 @@ def elliptical_radius(t, phi): return r -points, connectivity = spline.mesh(resolution=0.1, angular_resolution=36, radius=elliptical_radius, frame="frenet") +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) # %% -# Or you can change the radius along the spline. +# 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(resolution=0.1, angular_resolution=36, radius=radius, frame="frenet") +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) # %% -# 4. Bishop frame -# --------------- +# 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 +# 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. -# Let's create a spline with some unwated twist and see how the Bishop frame fixes it. +# +# 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], - [1.0, 0.0, 0.0], - [1.0, 1.0, 0.0], - [1.0, 1.0, 1.0], - [0.0, 1.0, 1.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(resolution=0.1, angular_resolution=36, radius=elliptical_radius, frame="frenet") +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 but requires an initial orientation for the ellipse, given by an initial vector. +# The Bishop frame eliminates this twist. -initial_vector = np.zeros(3) -tangent = spline.eval(0, derivative=1) -initial_vector[1] = tangent[2] -initial_vector[2] = -tangent[1] +points, connectivity = spline.mesh(step_t=0.1, step_angle=36, radius=elliptical_radius, frame="bishop") -points, connectivity = spline.mesh( - resolution=0.1, angular_resolution=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) # %% -# Changing the initial vector rotates the ellipse, because the initial vector is the reference for phi=0. +# Changing the initial vector rotates the frame and therefore the ellipse. initial_vector = np.array([0.5, -0.5, 1]) points, connectivity = spline.mesh( - resolution=0.1, angular_resolution=36, radius=elliptical_radius, frame="bishop", initial_vector=initial_vector + 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) # %% -# 5. Volume mesh +# 6. Volume Mesh # -------------- -# Besides a surface mesh, we can also turn the spline into a volume mesh. +# Finally, you can generate a volumetric mesh by setting the :code:`mesh_type` to :code:`"volume"`. points, connectivity = spline.mesh( - radius=radius, resolution=0.5, angular_resolution=72, initial_vector=initial_vector, mesh_type="volume" + 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) @@ -158,5 +167,5 @@ def radius(t, phi): # %% # Tips # ---- -# * To save the meshes for visulalisation in ParaView, you can use pyvista: :code:`mesh.save("mesh.vtk")` -# * By default the surface meshes are open at the end. You can close them using the :code:`cap_ends` keyword argument of :meth:`splinebox.spline_curves.Spline.mesh`. +# * 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. From c07a0a8e67404261a43bb3e8a4a7ba0e3685d07d Mon Sep 17 00:00:00 2001 From: faymanns Date: Mon, 6 Jan 2025 11:12:03 +0100 Subject: [PATCH 19/21] Change key word argument kind to method --- src/splinebox/spline_curves.py | 12 ++++++------ tests/test_spline_curves.py | 8 ++++---- 2 files changed, 10 insertions(+), 10 deletions(-) diff --git a/src/splinebox/spline_curves.py b/src/splinebox/spline_curves.py index a5b9595..39aa6f5 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,7 +649,7 @@ 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 @@ -661,7 +661,7 @@ def moving_frame(self, t, kind="frenet", initial_vector=None): t : np.array The parameters values of the spline for which the frame should be evaluated. - kind : string + method : 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 @@ -687,7 +687,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,7 +701,7 @@ 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: tangent = frame[0, 0] # Try to do the same as for the Frenet frame @@ -737,7 +737,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): diff --git a/tests/test_spline_curves.py b/tests/test_spline_curves.py index 1f576a9..5fe7290 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) @@ -800,13 +800,13 @@ def test_moving_frame(initialized_spline_curve, not_differentiable_twice): 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) From 637f559bad808dcc807aecc45c9a42adf862aa49 Mon Sep 17 00:00:00 2001 From: faymanns Date: Mon, 6 Jan 2025 12:46:18 +0100 Subject: [PATCH 20/21] Improve doc strings of mesh and moving_frame and fix kind in example --- examples/plot_frame.py | 6 +- src/splinebox/spline_curves.py | 169 ++++++++++++++++++++++++--------- 2 files changed, 129 insertions(+), 46 deletions(-) diff --git a/examples/plot_frame.py b/examples/plot_frame.py index 0cbef22..0847291 100644 --- a/examples/plot_frame.py +++ b/examples/plot_frame.py @@ -7,7 +7,7 @@ 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]_. +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 @@ -40,13 +40,13 @@ # 2. Compute the Frenet-Serret Frame # ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ -frenet_frame = spline.moving_frame(t, kind="frenet") +frenet_frame = spline.moving_frame(t, method="frenet") # %% # 3. Compute the Bishop Frame # ^^^^^^^^^^^^^^^^^^^^^^^^^^^ -bishop_frame = spline.moving_frame(t, kind="bishop") +bishop_frame = spline.moving_frame(t, method="bishop") # %% # 4. Visualize the Frames diff --git a/src/splinebox/spline_curves.py b/src/splinebox/spline_curves.py index 39aa6f5..1d906de 100644 --- a/src/splinebox/spline_curves.py +++ b/src/splinebox/spline_curves.py @@ -651,33 +651,65 @@ def normal(self, t, frame="bishop", 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. - method : 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. If None, an initial vector is computed automatically. + 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: @@ -925,38 +957,89 @@ def mesh( initial_vector=None, ): """ - Creats a mesh around the spline curve in 3D. The distance of the mesh - from the spline can be changes using the :code:`radius` paramters. - The mesh can either be a surface with triangular cells or a volume mesh - with tetrahedral cells. + 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 - A float value results in a constant radius. Alternatively, it can be a function that takes - the spline parameter t and the polar angle in the normal plane as arguments and returns - a float. - step_t : float - The step size of the spline parameter t. - step_angle : float - The step size of the polar angle. - mesh_type : str - Can be "surface" or "volume". - cap_ends : boolean - If True, the ends of the surface mesh are closed with orthogonal planes. - frame : str - Can be "frenet" or "bishop". See :meth:`splinebox.spline_curves.moving_frame`. - initial_vector : numpy array or None - The initial vector determining the orientation of the Bishop frame. - See :meth:`splinebox.spline_curves.moving_frame`. + 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 numpy array. The second dimension encodes the coordinates of each point in 3D space. + 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 numpy array. Each row corresponds to an element in the mesh. - The values are the indices of the points in :code:`points` that form the element. + 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.") From d91f6e21b206c41386a1749eb0fd5a89862ee9c7 Mon Sep 17 00:00:00 2001 From: faymanns Date: Mon, 6 Jan 2025 13:34:20 +0100 Subject: [PATCH 21/21] Improve test coverage --- src/splinebox/spline_curves.py | 4 ++-- tests/test_spline_curves.py | 21 +++++++++++++++++++++ 2 files changed, 23 insertions(+), 2 deletions(-) diff --git a/src/splinebox/spline_curves.py b/src/splinebox/spline_curves.py index 1d906de..d200d2f 100644 --- a/src/splinebox/spline_curves.py +++ b/src/splinebox/spline_curves.py @@ -1116,7 +1116,7 @@ def mesh( @staticmethod @numba.jit(nopython=True, nogil=True, cache=True) - def _surface_mesh_connectivity(closed, n_angles, n_t, n_points): + 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: @@ -1140,7 +1140,7 @@ def _surface_mesh_connectivity(closed, n_angles, n_t, n_points): @staticmethod @numba.jit(nopython=True, nogil=True, cache=True) - def _volume_mesh_connectivity(closed, n_angles, n_t, n_points): + 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: diff --git a/tests/test_spline_curves.py b/tests/test_spline_curves.py index 5fe7290..4c1dcff 100644 --- a/tests/test_spline_curves.py +++ b/tests/test_spline_curves.py @@ -820,6 +820,27 @@ def test_moving_frame(initialized_spline_curve, not_differentiable_twice): assert np.allclose(np.cross(frame[:, 0], spline.eval(t, derivative=1)), 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 ):