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?
Fundamentally, each successive three points must lie on their own circle and the local radius of curvature (reciprocal of the “curvature”) is just the radius of that circle.
Take points i-1, i, i+1 and form successive tangent vectors tA and tB.
Three points must lie in a unique plane, unless they are collinear.
Within that plane the three points must lie on a circle, whose centre is the intersection of the perpendicular bisectors; (again, unless they are collinear, in which case Rc is infinite.) Find this centre and its distance from any one point: this is the radius of curvature.
I found the centre by noting the vector equations of the two bisectors:
centre = mA + p nA = mB + q nB
where mA is the midpoint of the line joining the first two points and nA is in the direction of the perpendicular bisector to those two points; similarly mB and nB for the latter pair of points. Isolate p (say) by dot-producting with tB, the tangent vector for the later segment. Then you get the centre of the circle.
import numpy as np
def radius_of_curvature( points ): # Points are [ [x0,y0,z0], [x1,y1,z1], ... ]
N = len( points )
R = np.zeros( N )
for i in range( 1, N - 1 ):
tA, tB = points[i] - points[i-1], points[i+1] - points[i] # tangent vectors
z = np.cross( tA, tB ) # normal to plane
if np.linalg.norm( z ) < 1e-30: # HELP! points are co-linear
R[i] = 1e30 # my favourite approximation to infinity
continue
nA, nB = np.cross( z, tA ), np.cross( z, tB ) # normals to successive segments
p = np.dot( 0.5 * ( tA + tB ), tB ) / np.dot( nA, tB ) # parameter along normal A
C = 0.5 * ( points[i-1] + points[i] ) + p * nA # centre of circle
R[i] = np.linalg.norm( C - points[i] ) # radius of circle
R[0], R[-1] = R[1], R[-2] # arbitrary choice at end points
return R
def generate_circle_by_angles(t, C, r, theta, phi):
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_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]) # centre
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 = radius_of_curvature(p)
print( Rs )
Output:
[2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5
2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5
2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5
2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5
2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5
2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5]