How can I get the local curvature of a scipy.CubicSpline?

Question

I'm using the scipy.interpolate.CubicSpline to compute a 2d function and want to analyse the curvature of the graph. The CubicSpline can compute a first and second derivative and my approach was to use curvature as k(t) = |f''(t)| / (1+f'(t)**2)**1.5 (https://math.stackexchange.com/questions/1155398/difference-between-second-order-derivative-and-curvature)

I then visualize the curvature by plotting a circle with r=1/r, which looks reasonable at first sight (larger circle at flatter regions), but is off noticeably.

The plot looks like this:

The point density on the spline also depends on the curvature and I suspect that this 'velocity' influences the curvature computation as well.

The plot was generated with this small script:

import numpy as np
from scipy.interpolate import CubicSpline
import matplotlib.pyplot as plt

start = [0, 0]
end = [0.6, 0.2]

d_start = np.array([3.0, 0.0])
d_end = np.array([3.0, 0.0])

cs = CubicSpline([0, 1], [start, end], bc_type=((1, d_start), (1, d_end)))

samples = np.linspace(0, 1, 100)
positions = cs(samples)

plt.plot(positions[:, 0], positions[:, 1], 'bx-')
plt.axis('equal')

t = samples[28]

touch_point = cs(t) # circle should be tangent to the curve at this point
tangent = cs(t, nu=1)
tangent_normed = tangent / np.linalg.norm(tangent)

cs1 = np.linalg.norm(cs(t, nu=1))
cs2 = np.linalg.norm(cs(t, nu=2))

k = abs(cs2) / ((1 + cs1**2)**1.5)

r = 1/k
center = touch_point + r * np.array([-tangent_normed[1], tangent_normed[0]])
plt.plot(touch_point[0], touch_point[1], 'go')
circle = plt.Circle(center, r, color='r', fill=True)
plt.gca().add_artist(circle)

plt.show()

Answer 1

The formula used is only a special case for a function in the form of x=t, y=f(t). In the 2d version f(t)=(x(t), y(t)) the general formula is noted below.

See: https://en.wikipedia.org/wiki/Curvature#In_terms_of_arc-length_parametrization

in my case, the curvature can be computed as:

x_, y_ = cs(t, nu=1)
x__, y__ = cs(t, nu=2)

k = (x_ * y__ - y_ * x__) / (x_**2 + y_**2)**1.5
r = 1/k 

and now the circle fits nicely to the plot