如何获得 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）

然后，我通过绘制一个 r=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) 形式的函数的特例。在 2d 版本中，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 

现在圆圈非常适合情节