均匀三次b样条曲线_B样条

B-Spline

Beizer Curve

在之前我写过一篇关于 Beizer Curve 的文章，即使现在看起来依旧觉得 Beizer Curve 按照interpolation 定义推导出来结果真是美啊：

解：

这个有好几种看到方式，比如我们把

看作基，然后我们就是用

作为 weights 来控制 Beizer Curve 的形状：

或者把它写成：

本质上也就是多项式，而且可以验证：

• t = 0,
• t = 1,

曲线并不过

. 几种常见的曲线形状：

图片来自《Mathematics for 3D Game Programming and Computer Graphics》

其实我的最爱还是写成矩阵形式：

B-Spline

样条的定义是： (n+1)个点， 包含

段，其中每一段都是由

来定义的，所以有：

样条为了美观我们需要满足两个条件：

• 连续
• 曲线应该跟坐标系选择无关

由以上两个条件我们可以解得样条每一段的函数表示：

依旧可以模仿 Beizer Curve 我们把它写成：

它的基可以看作：

同时我们可以算出开始点和结束点：

不同于 Bezier Curve，这里是完全不过控制点的:

图片来自《Mathematics for 3D Game Programming and Computer Graphics》

问题在于大部分情况我们并不是给出了控制点，给出的是曲线上的点（如上图的knots），给出knots我们想求B-Spline，假设我们有 n+1 个 knots:

, 那么根据上面的 1式和2式 我们可以列 n+1 个等式，控制点要多两个

，这里我们一般的办法是我们给出开始点和结束点的切线

然后就可以写出矩阵来求解:

图片来自《Curves and Surfaces for Computer Graphics》

或者：

图片来自《Curves and Surfaces for Computer Graphics》

Scipy 库

我们也先来看一下用样条的原因，假设我这里有6个点，其实当然也可以用6次多项式来拟合：

高次多项式的不好之处在于它波动太大了。

interp1d

使用 interp1d

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


points = [[0,1], [1, 1.1],[2, 1.2],[3,3],[4,2.9],[5,2.8],[6,2.7]]

x = [point[0] for point in points]
y = [point[1] for point in points]

f = np.polyfit(x, y, deg=6)
p = np.poly1d(f)

ls = interpolate.interp1d(x, y, kind = 'linear')
qs = interpolate.interp1d(x, y, kind = 'quadratic')
cs = interpolate.interp1d(x, y, kind = 'cubic')

x_new = np.linspace(x[0], x[-1], 100)

plt.scatter(x, y, label='data')
plt.plot(x_new, p(x_new), label = 'polynomial p6')
plt.plot(x_new, ls(x_new), label = 'linear spline')
plt.plot(x_new, qs(x_new), label = 'quadratic spline')
plt.plot(x_new, cs(x_new), label = 'cubic spline')

plt.legend(loc = 'upper left')

可以看到使用 interp1d 的 linear, quadratic, cubic 的不同效果，其实 quadratic 在这里也还表现不错。

发现还有一个scipy.interpolate.CubicSpline, 效果基本等同于 interpolate.interp1d(x, y, kind = 'cubic').

Akima1DInterpolator

用这个 Akima1DInterpolator 来看这个例子的话:

(⊙o⊙)…，不过大概是这个例子比较特殊，如果我们用最常用的正弦函数的话表现还可以，真是o(╯□╰)o：

splprep

然后发现 Scipy 还有一个 interpolate.splXXX, 按照官方的 doc 说法要用 spline interpolation需要两步：

1. a spline representation of the curve is computed
2. the spline is evaluated at the desired points

使用

spl = interpolate.splrep(x,y)
spline = interpolate.splev(x_new,spl)

和 cubic spline 的效果依旧类似，但是如果我们使用interpolate.splprep 看起来就不一样了，(๑˘ ˘๑)

tck,u = interpolate.splprep([x,y], k = 3, s = 0)
spline = interpolate.splev(np.linspace(0,1,100),tck)

plt.scatter(x, y, label='data')
plt.plot(x_new, cs(x_new), label = 'cubic spline')
plt.plot(x_new, f(x_new), label = 'Akima Interpolator')
plt.plot(spline[0], spline[1], label = "splprep")
plt.legend(loc = 'upper left')

Σ(っ °Д °;)っ 神奇！ 感觉这个绿色的 splprep 看起来更自然一点啊，文档上：

• splprep: Find the B-spline representation of an N-dimensional curve.
• splrep: Find the B-spline representation of 1-D curve.

N-dimensional vs 1-D，感觉这个地方我看的也不是很懂(´◔ω◔)

使用 splprep 还有的附加好处是如果我们的点并不是向右的，我们传入的x不是递增的，那么 CubicSpline 会报错 ‘x must be strictly increasing sequence’, 而 splprep 在这种情况下依旧能过输出合理的结果。

参考：

• pythonでxy座標上の離散点をスプライン補間
• <Curves and Surfaces for Computer Graphics>
• <Mathematics for 3D Game Programming and Computer Graphics, Third Edition>