0
votes

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:

Planes following spline

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

1 Answers

0
votes

I think I ended up producing something that seems to work in the end:

import numpy as np
import vg
from pytransform3d.rotations import matrix_from_axis_angle

def _rotmat(self, vector, points):
    """
    Rotates a 3xn array of 3D coordinates from the +z normal to an
    arbitrary new normal vector.
    """
    
    vector = vg.normalize(vector)
    axis = vg.perpendicular(vg.basis.z, vector)
    angle = vg.angle(vg.basis.z, vector, units='rad')
    
    a = np.hstack((axis, (angle,)))
    R = matrix_from_axis_angle(a)
    
    r = Rot.from_matrix(R)
    rotmat = r.apply(points)
    
    return rotmat

Not too insanely complicated, just start with a plane of points aligned with the x-y plane (assuming you're using x-y as your horizontal like me here apparently, please don't hate me), then it'll rotate it along the vector and not really care about rotation about the axis. Seems to work ok.