python - 如何评估 splprep 的样条导数？

Question

我希望使用 B 样条曲线计算 numpy 中闭合曲线的曲率。我想评估样条表示的导数而不是数据以获得平滑的结果。但是下面的代码返回错误:

    Traceback (most recent call last):
      File "<stdin>", line 16, in <module>
      File "/Users/jfl/anaconda/lib/python2.7/site-packages/scipy/interpolate/fitpack.py", line 1212, in splder
        c = (c[1:-1-k] - c[:-2-k]) * k / dt
    TypeError: unsupported operand type(s) for -: 'list' and 'list'

shell returned 1

所以我想知道我是否正确使用了 splder 函数...

import numpy as np
import scipy.interpolate as intplt
import matplotlib.pyplot as plt


t = np.linspace(0,2*np.pi,100)
x = np.cos(t)
y = np.sin(t)


pts = np.vstack((x,y))
tck, u = intplt.splprep(pts, u=None,k=3, s=0.0, per=1)
u_new = np.linspace(u.min(), u.max(), 1000)
x_new, y_new = intplt.splev(u_new, tck, der=0)

tck_der1 = intplt.splder(tck)
tck_der2 = intplt.splder(tck_der1)


xp, yp = intplt.splev(u_new, tck_der1, der=0)
xpp, ypp = intplt.splev(u_new, tck_der2, der=0)



plt.figure()
plt.plot(x,y,".")
plt.plot(x_new,y_new)

plt.figure()
curvature = np.abs(xp* ypp - yp* xpp) / np.power(xp** 2 + yp** 2, 3 / 2)
plt.plot(u_new,curvature)

plt.show()

Answer 1

splder 仅支持标量样条线(使用 splrep 插值函数 y = f(x))，不支持使用 获得的矢量样条线>splprep。无论如何，您都不需要它:只需使用 splev 的 der 参数即可:

xp, yp = intplt.splev(u_new, tck, der=1)
xpp, ypp = intplt.splev(u_new, tck, der=2)

这不是导数的通用有限差分评估； splev 实际上会使用样条结构来计算它，这正是您想要做的。

我猜您尝试了上述操作并对输出不满意:

但这只是 Matplotlib 的搞笑而已。绘图窗口从+9.998e-1 到+9.998e-1 + 0.0006。换句话说，曲率几乎恒定 (1)，因此 matplotlib 放大了来自插值的不可避免的噪声。只要设置合理的绘图窗口，问题就消失了。

(或者，使用诸如 x = 2*np.cos(t) 作为示例，以获得漂亮的非常数曲率。)