如何正确计算scipy.CubicSpline的局部曲率？

问题背景 #

使用scipy.interpolate.CubicSpline拟合二维参数化曲线，尝试通过一阶、二阶导数计算曲率并绘制曲率圆，但结果存在明显偏差，怀疑是样条参数化（点密度/速度）影响计算结果。

错误原因 #

你误用了一元函数的曲率公式，但这里的曲线是二维参数化曲线（以t为参数，r(t)=(x(t), y(t))），两者的曲率计算公式完全不同：

• 一元函数y=f(x)的曲率公式：k = |f''(x)| / (1 + f'(x)²)^(3/2)
• 二维参数化曲线的正确曲率公式：

  k(t) = |x’(t)y''(t) - x''(t)y’(t)| / (x’(t)² + y’(t)²)^(3/2)
  

  分子部分是一阶导数和二阶导数的二维叉积（行列式）的绝对值，反映曲线的转向速率；分母是一阶导数模长的三次方，用于消除参数化速度的干扰——这正是你怀疑的“点密度/速度”问题：参数化曲线的曲率是几何属性，与参数t的选取无关，正确公式会自动消除参数速度的影响。

修正后的代码 #

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)
tangent = cs(t, nu=1)
tangent_normed = tangent / np.linalg.norm(tangent)

# 计算参数化曲线的曲率
x1, y1 = tangent
x2, y2 = cs(t, nu=2)
cross = x1 * y2 - x2 * y1  # 二维叉积
denominator = (x1**2 + y1**2) ** 1.5
k = abs(cross) / denominator

r = 1/k
# 曲率中心方向：垂直于切线，指向曲线凹侧（方向相反则替换为[tangent_normed[1], -tangent_normed[0]]）
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, alpha=0.3)
plt.gca().add_artist(circle)

plt.show()

关键修正点 #

1. 替换曲率计算逻辑：用一阶导数和二阶导数的叉积（x1*y2 - x2*y1）替代原代码中二阶导数的模长，这是参数化曲线曲率计算的核心。
2. 分母使用一阶导数模长的三次方，自动消除参数t的速度影响，确保曲率是曲线的几何属性，与采样密度无关。
3. 曲率中心的方向需根据曲线凹侧调整，若生成的圆在曲线凸侧，只需将方向向量改为[tangent_normed[1], -tangent_normed[0]]即可。

内容的提问来源于stack exchange，提问作者FooTheBar