Python中的B样条插值
我正在尝试用Python重现一个Mathematica的B样条示例。
Mathematica示例的代码是：
pts = {{0, 0}, {0, 2}, {2, 3}, {4, 0}, {6, 3}, {8, 2}, {8, 0}};
Graphics[{BSplineCurve[pts, SplineKnots -> {0, 0, 0, 0, 2, 3, 4, 6, 6, 6, 6}], Green, Line[pts], Red, Point[pts]}]
运行后得到了我预期的结果。现在我想用Python/scipy来做同样的事情：
import numpy as np
import matplotlib.pyplot as plt
import scipy.interpolate as si
points = np.array([[0, 0], [0, 2], [2, 3], [4, 0], [6, 3], [8, 2], [8, 0]])
x = points[:,0]
y = points[:,1]
t = range(len(x))
knots = [2, 3, 4]
ipl_t = np.linspace(0.0, len(points) - 1, 100)
x_tup = si.splrep(t, x, k=3, t=knots)
y_tup = si.splrep(t, y, k=3, t=knots)
x_i = si.splev(ipl_t, x_tup)
y_i = si.splev(ipl_t, y_tup)
print 'knots:', x_tup
fig = plt.figure()
ax = fig.add_subplot(111)
plt.plot(x, y, label='original')
plt.plot(x_i, y_i, label='spline')
plt.xlim([min(x) - 1.0, max(x) + 1.0])
plt.ylim([min(y) - 1.0, max(y) + 1.0])
plt.legend()
plt.show()
结果是得到了一些插值，但看起来不太对。我分别对x和y的部分进行参数化和样条处理，使用的节点和Mathematica是一样的。然而，我得到了过高和过低的值，这导致我的插值曲线在控制点的凸包外面弯曲。正确的做法是什么？Mathematica是怎么做到的？
3 个回答
1
我觉得scipy的fitpack库做的事情比Mathematica要复杂一些。我自己也搞不太明白到底发生了什么。
这些函数里有一个平滑参数，默认的插值行为是试图让点通过直线。这就是fitpack软件的工作原理，所以我猜scipy只是继承了这个功能？（http://www.netlib.org/fitpack/all -- 我不确定这是不是正确的fitpack）
我从http://research.microsoft.com/en-us/um/people/ablake/contours/那里获取了一些想法，并用里面的B样条代码实现了你的例子。
import numpy
import matplotlib.pyplot as plt
# This is the basis function described in eq 3.6 in http://research.microsoft.com/en-us/um/people/ablake/contours/
def func(x, offset):
out = numpy.ndarray((len(x)))
for i, v in enumerate(x):
s = v - offset
if s >= 0 and s < 1:
out[i] = s * s / 2.0
elif s >= 1 and s < 2:
out[i] = 3.0 / 4.0 - (s - 3.0 / 2.0) * (s - 3.0 / 2.0)
elif s >= 2 and s < 3:
out[i] = (s - 3.0) * (s - 3.0) / 2.0
else:
out[i] = 0.0
return out
# We have 7 things to fit, so let's do 7 basis functions?
y = numpy.array([0, 2, 3, 0, 3, 2, 0])
# We need enough x points for all the basis functions... That's why the weird linspace max here
x = numpy.linspace(0, len(y) + 2, 100)
B = numpy.ndarray((len(x), len(y)))
for k in range(len(y)):
B[:, k] = func(x, k)
plt.plot(x, B.dot(y))
# The x values in the next statement are the maximums of each basis function. I'm not sure at all this is right
plt.plot(numpy.array(range(len(y))) + 1.5, y, '-o')
plt.legend('B-spline', 'Control points')
plt.show()
for k in range(len(y)):
plt.plot(x, B[:, k])
plt.title('Basis functions')
plt.show()
总之，我觉得其他人也有同样的问题，可以看看这个链接：scipy的splrep行为
13
我写了这个函数，是为了回答我在这里提的另一个问题。
在我的问题中，我想找一些方法来用scipy计算B样条（这也是我偶然发现你们问题的原因）。
经过一番研究，我想出了下面这个函数。它可以评估任何曲线，最高支持到20次方（其实远远超过我们的需求）。在速度方面，我测试了100,000个样本，结果只用了0.017秒。
import numpy as np
import scipy.interpolate as si
def bspline(cv, n=100, degree=3, periodic=False):
""" Calculate n samples on a bspline
cv : Array ov control vertices
n : Number of samples to return
degree: Curve degree
periodic: True - Curve is closed
False - Curve is open
"""
# If periodic, extend the point array by count+degree+1
cv = np.asarray(cv)
count = len(cv)
if periodic:
factor, fraction = divmod(count+degree+1, count)
cv = np.concatenate((cv,) * factor + (cv[:fraction],))
count = len(cv)
degree = np.clip(degree,1,degree)
# If opened, prevent degree from exceeding count-1
else:
degree = np.clip(degree,1,count-1)
# Calculate knot vector
kv = None
if periodic:
kv = np.arange(0-degree,count+degree+degree-1,dtype='int')
else:
kv = np.concatenate(([0]*degree, np.arange(count-degree+1), [count-degree]*degree))
# Calculate query range
u = np.linspace(periodic,(count-degree),n)
# Calculate result
return np.array(si.splev(u, (kv,cv.T,degree))).T
这是对开放曲线和周期曲线的结果：
cv = np.array([[ 50., 25.],
[ 59., 12.],
[ 50., 10.],
[ 57., 2.],
[ 40., 4.],
[ 40., 14.]])
27
我用Python和scipy重新做了之前提到的Mathematica例子，结果如下：
B样条，非周期性
关键在于要么拦截系数，也就是从scipy.interpolate.splrep返回的元组中的第一个元素，然后把它们替换成控制点的值，再传给scipy.interpolate.splev；要么如果你愿意自己创建节点，也可以不使用splrep，自己创建整个元组。
不过奇怪的是，根据手册，splrep返回的（而splev期望的）元组中包含一个样条系数向量，每个节点对应一个系数。然而，根据我找到的所有资料，样条是定义为N_control_points个基样条的加权和，所以我本以为系数向量的元素数量应该和控制点一样，而不是节点位置。
实际上，当把splrep的结果元组中修改过的系数向量传给scipy.interpolate.splev时，发现这个向量的前N_control_points个元素确实是N_control_points个基样条的预期系数。最后degree + 1个元素似乎没有任何效果。我对为什么会这样感到困惑。如果有人能解释一下，那就太好了。下面是生成上述图形的源代码：
import numpy as np
import matplotlib.pyplot as plt
import scipy.interpolate as si
points = [[0, 0], [0, 2], [2, 3], [4, 0], [6, 3], [8, 2], [8, 0]];
points = np.array(points)
x = points[:,0]
y = points[:,1]
t = range(len(points))
ipl_t = np.linspace(0.0, len(points) - 1, 100)
x_tup = si.splrep(t, x, k=3)
y_tup = si.splrep(t, y, k=3)
x_list = list(x_tup)
xl = x.tolist()
x_list[1] = xl + [0.0, 0.0, 0.0, 0.0]
y_list = list(y_tup)
yl = y.tolist()
y_list[1] = yl + [0.0, 0.0, 0.0, 0.0]
x_i = si.splev(ipl_t, x_list)
y_i = si.splev(ipl_t, y_list)
#==============================================================================
# Plot
#==============================================================================
fig = plt.figure()
ax = fig.add_subplot(231)
plt.plot(t, x, '-og')
plt.plot(ipl_t, x_i, 'r')
plt.xlim([0.0, max(t)])
plt.title('Splined x(t)')
ax = fig.add_subplot(232)
plt.plot(t, y, '-og')
plt.plot(ipl_t, y_i, 'r')
plt.xlim([0.0, max(t)])
plt.title('Splined y(t)')
ax = fig.add_subplot(233)
plt.plot(x, y, '-og')
plt.plot(x_i, y_i, 'r')
plt.xlim([min(x) - 0.3, max(x) + 0.3])
plt.ylim([min(y) - 0.3, max(y) + 0.3])
plt.title('Splined f(x(t), y(t))')
ax = fig.add_subplot(234)
for i in range(7):
vec = np.zeros(11)
vec[i] = 1.0
x_list = list(x_tup)
x_list[1] = vec.tolist()
x_i = si.splev(ipl_t, x_list)
plt.plot(ipl_t, x_i)
plt.xlim([0.0, max(t)])
plt.title('Basis splines')
plt.show()
B样条，周期性
现在，为了创建一个像下面这样的闭合曲线，这是另一个可以在网上找到的Mathematica例子，
如果你使用splrep，需要在调用时设置per参数。在控制点列表的末尾添加degree+1个值后，这似乎效果不错，如图所示。
不过这里的另一个特殊情况是，系数向量中的第一个和最后degree个元素没有效果，这意味着控制点必须从向量的第二个位置开始放入，也就是位置1。只有这样结果才会正常。对于k=4和k=5的情况，这个位置甚至变成了位置2。
下面是生成闭合曲线的源代码：
import numpy as np
import matplotlib.pyplot as plt
import scipy.interpolate as si
points = [[-2, 2], [0, 1], [-2, 0], [0, -1], [-2, -2], [-4, -4], [2, -4], [4, 0], [2, 4], [-4, 4]]
degree = 3
points = points + points[0:degree + 1]
points = np.array(points)
n_points = len(points)
x = points[:,0]
y = points[:,1]
t = range(len(x))
ipl_t = np.linspace(1.0, len(points) - degree, 1000)
x_tup = si.splrep(t, x, k=degree, per=1)
y_tup = si.splrep(t, y, k=degree, per=1)
x_list = list(x_tup)
xl = x.tolist()
x_list[1] = [0.0] + xl + [0.0, 0.0, 0.0, 0.0]
y_list = list(y_tup)
yl = y.tolist()
y_list[1] = [0.0] + yl + [0.0, 0.0, 0.0, 0.0]
x_i = si.splev(ipl_t, x_list)
y_i = si.splev(ipl_t, y_list)
#==============================================================================
# Plot
#==============================================================================
fig = plt.figure()
ax = fig.add_subplot(231)
plt.plot(t, x, '-og')
plt.plot(ipl_t, x_i, 'r')
plt.xlim([0.0, max(t)])
plt.title('Splined x(t)')
ax = fig.add_subplot(232)
plt.plot(t, y, '-og')
plt.plot(ipl_t, y_i, 'r')
plt.xlim([0.0, max(t)])
plt.title('Splined y(t)')
ax = fig.add_subplot(233)
plt.plot(x, y, '-og')
plt.plot(x_i, y_i, 'r')
plt.xlim([min(x) - 0.3, max(x) + 0.3])
plt.ylim([min(y) - 0.3, max(y) + 0.3])
plt.title('Splined f(x(t), y(t))')
ax = fig.add_subplot(234)
for i in range(n_points - degree - 1):
vec = np.zeros(11)
vec[i] = 1.0
x_list = list(x_tup)
x_list[1] = vec.tolist()
x_i = si.splev(ipl_t, x_list)
plt.plot(ipl_t, x_i)
plt.xlim([0.0, 9.0])
plt.title('Periodic basis splines')
plt.show()
B样条，周期性，高阶
最后，还有一个我无法解释的现象，就是当阶数达到5时，样条曲线出现了一个小的间断，见右上角的面板，这是那个“半月形带鼻子”的特写。产生这个效果的源代码如下。
import numpy as np
import matplotlib.pyplot as plt
import scipy.interpolate as si
points = [[-2, 2], [0, 1], [-2, 0], [0, -1], [-2, -2], [-4, -4], [2, -4], [4, 0], [2, 4], [-4, 4]]
degree = 5
points = points + points[0:degree + 1]
points = np.array(points)
n_points = len(points)
x = points[:,0]
y = points[:,1]
t = range(len(x))
ipl_t = np.linspace(1.0, len(points) - degree, 1000)
knots = np.linspace(-degree, len(points), len(points) + degree + 1).tolist()
xl = x.tolist()
coeffs_x = [0.0, 0.0] + xl + [0.0, 0.0, 0.0]
yl = y.tolist()
coeffs_y = [0.0, 0.0] + yl + [0.0, 0.0, 0.0]
x_i = si.splev(ipl_t, (knots, coeffs_x, degree))
y_i = si.splev(ipl_t, (knots, coeffs_y, degree))
#==============================================================================
# Plot
#==============================================================================
fig = plt.figure()
ax = fig.add_subplot(231)
plt.plot(t, x, '-og')
plt.plot(ipl_t, x_i, 'r')
plt.xlim([0.0, max(t)])
plt.title('Splined x(t)')
ax = fig.add_subplot(232)
plt.plot(t, y, '-og')
plt.plot(ipl_t, y_i, 'r')
plt.xlim([0.0, max(t)])
plt.title('Splined y(t)')
ax = fig.add_subplot(233)
plt.plot(x, y, '-og')
plt.plot(x_i, y_i, 'r')
plt.xlim([min(x) - 0.3, max(x) + 0.3])
plt.ylim([min(y) - 0.3, max(y) + 0.3])
plt.title('Splined f(x(t), y(t))')
ax = fig.add_subplot(234)
for i in range(n_points - degree - 1):
vec = np.zeros(11)
vec[i] = 1.0
x_i = si.splev(ipl_t, (knots, vec, degree))
plt.plot(ipl_t, x_i)
plt.xlim([0.0, 9.0])
plt.title('Periodic basis splines')
plt.show()
考虑到b样条在科学界的广泛应用，以及scipy是如此全面的工具箱，而我在网上找不到太多关于我在问的问题的信息，这让我觉得我可能走错了方向或者忽略了什么。如果有人能提供帮助，我将非常感激。