So I'm currently trying to take slices on a plane orthogonal to a spline. Direction doesn't really matter too much since I'm using the points to interpolate 3D scans
I'm mainly unsure about the rotmat method (this is a stripped down version of my class, technically a NURBS-Python surface derived class), where I'm rotating the plane mesh from a flat x/y plane (all z=0) to match the new normal vector (tangent of the spline, stored in the der variable).
Anyone have an idea how to rotate a set of points to go from one normal vector to another? The angle around the axis of the new vector doesn't matter than much to me.
(sorry for vg, kind of an obscure library but somewhat convenient actually):
from scipy.interpolate import splprep, splev
import numpy as np
import vg
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from scipy.spatial.transform import Rotation as R
class SplineTube():
_points = np.array(
[[0, 0, 0],
[0, 1, 0],
[1, 1, 0],
[1, 0, 0]],
) - np.array([0.5, 0.5, 0])
_normal = np.array([0, 0, 1])
def __init__(self, x, y, z, n = 3, degree=2, **kwargs):
assert n >= 3
tck, u = splprep([x, y, z], s=0, k=2)
evalpts = np.linspace(0, 1, n)
pts = np.array(splev(evalpts, tck))
der = np.array(splev(evalpts, tck, der=1))
points = []
for i in range(n):
points_slice = self.rotmat(der[:, i], self._points)
points_slice = points_slice + pts[:, i]
points.append(points_slice)
points = np.stack(points)
return points
def rotmat(self, vector, points):
perpen = vg.perpendicular(self._normal, vector)
r = R.from_rotvec(perpen)
rotmat = r.apply(points)
return rotmat
Here's an example where I used a meshgrid instead of the _points, but is very similar:
x = [0, 1, 2, 3, 6]
y = [0, 2, 5, 6, 2]
z = [0, 3, 5, 7, 10]
tck, u = splprep([x, y, z], s=0, k=2)
evalpts = np.linspace(0, 1, 10)
pts = splev(evalpts, tck)
der = splev(evalpts, tck, der=1)
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot(pts[0], pts[1], pts[2])
ax.quiver(*pts, *der, length=0.05)
ax.scatter(x, y, z)
planes = SplineTube(x, y, z, n=10)
ax.scatter(planes[:, :, 0], planes[:, :, 1], planes[:, :, 2])