如何获取scipy.CubicSpline的局部曲率？

Question

我正在使用scipy.interpolate.CubicSpline来计算一个二维函数，并想分析这个图形的曲率。CubicSpline可以计算一阶和二阶导数，我的想法是用曲率来表示：

k(t) = |f''(t)| / (1+f'(t)**2)**1.5

(https://math.stackexchange.com/questions/1155398/difference-between-second-order-derivative-and-curvature)

然后我通过绘制一个半径为1/r的圆来可视化曲率，乍一看这个方法似乎合理（在比较平坦的区域，圆会更大），但实际上有明显的偏差。

这个图看起来是这样的：

样条曲线上的点密度也与曲率有关，我怀疑这种“速度”也会影响曲率的计算。

这个图是用下面这个小脚本生成的：

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

这个公式其实是一个特定情况，适用于一种形式的函数，具体是x=t，y=f(t)。在二维情况下，f(t)=(x(t), y(t))，下面是更一般的公式。

可以参考这里: https://en.wikipedia.org/wiki/Curvature#In_terms_of_arc-length_parametrization

在我的例子中，曲率可以这样计算：

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 

现在这个圆形很好地适应了图表。