[Python] Calculating the curvature of a discrete function in 3D space

I have a set of points representing a curve in 3D space. The goal is to detect the point with the maximum curvature.

When looking on the curvature page on Wikipedia, I find the curvature can be found as the magnitude of the acceleration of the parametric function.

My idea of solving the solution is to interpolate a B-spline over the 3D points. Next, I discretize the function with equal spaced points. At last, I calculate the acceleration (second derivation) over the discrete points. The curvature should be found by the magnitude of the acceleration. The Radius curvature can be found by inverting the curvature.

This is my understanding of the problem. If I'm incorrect, please correct me.

I have written a python function for this problem.

def curvature(points: np.ndarray) -> np.ndarray:
tck, u = splprep(points.T)
t = np.linspace(0, 1, len(points))
x, y, z = splev(t, tck)
parametric_points = np.stack([x, y, z]).T

tangent = np.diff(parametric_points, axis=0)

acceleration = np.diff(tangent, axis=0)

magnitude = np.array([np.linalg.norm(a) for a in acceleration])
radius_curvature = 1 / magnitude
return radius_curvature


To test my function, I generate a circle in 3D space with the code below and test to see the radius curvature for each point:

def generate_circle_by_angles(t, C, r, theta, phi):
# Source: https://meshlogic.github.io/posts/jupyter/curve-fitting/fitting-a-circle-to-cluster-of-3d-points/
# Orthonormal vectors n, u, <n,u>=0
n = np.array([np.cos(phi) * np.sin(theta), np.sin(phi) * np.sin(theta), np.cos(theta)])
u = np.array([-np.sin(phi), np.cos(phi), 0])

# P(t) = r*cos(t)*u + r*sin(t)*(n x u) + C
p_circle = r * np.cos(t)[:, np.newaxis] * u + r * np.sin(t)[:, np.newaxis] * np.cross(n, u) + C
return p_circle

r = 2.5 # Radius
c = np.array([3, 3, 4]) # Center
theta = 0 / 180 * np.pi # Azimuth
phi = 0 / 180 * np.pi # Zenith

t = np.linspace(0, np.pi, 100)
p = generate_circle_by_angles(t, c, r, theta, phi)

Rs = curvature(p)


However, the result is not correct.

[252.36949094 256.16299957 ... 260.04828741 256.16299957 252.36949094]


Has anyone a solution or remarks on my solution?

Continue reading...