使用 Python 进行 B 样条插值
声明：本页面是StackOverFlow热门问题的中英对照翻译，遵循CC BY-SA 4.0协议，如果您需要使用它，必须同样遵循CC BY-SA许可，注明原文地址和作者信息，同时你必须将它归于原作者（不是我）：StackOverFlow
原文地址: http://stackoverflow.com/questions/24612626/
Warning: these are provided under cc-by-sa 4.0 license. You are free to use/share it, But you must attribute it to the original authors (not me):
StackOverFlow
B-spline interpolation with Python
提问by zeus300
I am trying to reproduce a Mathematica example for a B-spline with Python.
我正在尝试使用 Python 为 B 样条重现 Mathematica 示例。
The code of the mathematica example reads
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]}]
and produces what I expect. Now I try to do the same with Python/scipy:
并产生我所期望的。现在我尝试用 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()
This results in somethng that is also interpolated but doesn't look quite right. I parameterize and spline the x- and y-components separately, using the same knots as mathematica. However, I get over- and undershoots, which make my interpolated curve bow outside of the convex hull of the control points. What's the right way to do this/how does mathematica do it?
这会导致一些也被插值但看起来不太正确的东西。我使用与 mathematica 相同的节点分别对 x 和 y 分量进行参数化和样条化。然而，我得到了过冲和下冲，这使得我的插值曲线在控制点的凸包之外弯曲。这样做的正确方法是什么/mathematica 是如何做到的？
回答by zeus300
I was able to recreate the Mathematica example I asked about in the previous post using Python/scipy. Here's the result:
我能够使用 Python/scipy 重新创建我在上一篇文章中询问的 Mathematica 示例。结果如下：
B-Spline, Aperiodic
B样条，非周期性
The trick was to either intercept the coefficients, i.e. element 1 of the tuple returned by scipy.interpolate.splrep, and to replace them with the control point values before handing them to scipy.interpolate.splev, or, if you are fine with creating the knots yourself, you can also do without splrepand create the entire tuple yourself.
诀窍是要么截取系数，即由 返回的元组的元素 1 scipy.interpolate.splrep，并在将它们交给 之前用控制点值替换它们scipy.interpolate.splev，或者，如果您自己创建结很好，您也可以不做splrep和自己创建整个元组。
What is strange about this all, though, is that, according to the manual, splrepreturns (and splevexpects) a tuple containing, among others, a spline coefficients vector with one coefficient per knot. However, according to all sources I found, a spline is defined as the weighted sum of the N_control_pointsbasis splines, so I would expect the coefficients vector to have as many elements as control points, not knot positions.
然而，这一切的奇怪之处在于，根据手册，splrep返回（并splev期望）一个元组，其中包含一个每个结一个系数的样条系数向量。然而，根据我发现的所有来源，样条被定义为N_control_points基础样条的加权和，所以我希望系数向量具有与控制点一样多的元素，而不是节点位置。
In fact, when supplying splrep's result tuple with the coefficients vector modified as described above to scipy.interpolate.splev, it turns out that the first N_control_pointsof that vector actually are the expected coefficients for the N_control_pointsbasis splines. The last degree + 1elements of that vector seem to have no effect. I'm stumped as to why it's done this way. If anyone can clarify that, that would be great. Here's the source that generates the above plots:
事实上，当提供splrep具有如上所述修改为 的系数向量的结果元组时scipy.interpolate.splev，结果证明该向量的前N_control_points实际上是N_control_points基样条的预期系数。该向量的最后一个度 + 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-Spline, Periodic
B样条，周期性
Now in order to create a closed curve like the following, which is another Mathematica example that can be found on the web,
现在为了创建如下所示的闭合曲线，这是另一个可以在网上找到的 Mathematica 示例，
it is necessary to set the perparameter in the splrepcall, if you use that. After padding the list of control points with degree+1values at the end, this seems to work well enough, as the images show.
如果您使用该per参数，则有必要在splrep调用中设置该参数。在最后用degree+1值填充控制点列表后，这似乎工作得很好，如图所示。
The next peculiarity here, however, is that the first and the last degreeelements in the coefficients vector have no effect, meaning that the control points must be put in the vector starting at the second position, i.e. position 1. Only then are the results ok. For degrees k=4 and k=5, that position even changes to position 2.
然而，这里的下一个特点是系数向量中的第一个和最后一个度元素没有影响，这意味着控制点必须放在从第二个位置开始的向量中，即位置 1。只有这样才能得到结果好的。对于度数 k=4 和 k=5，该位置甚至会更改为位置 2。
Here's the source for generating the closed curve:
这是生成闭合曲线的来源：
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-Spline, Periodic, Higher Degree
B样条，周期性，更高阶
Lastly, there is an effect that I can not explain either, and this is when going to degree 5, there is a small discontinuity that appears in the splined curve, see the upper right panel, which is a close-up of that 'half-moon-with-nose-shape'. The source code that produces this is listed below.
最后，还有一个我也无法解释的效果，这是在达到 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()
Given that b-splines are ubiquitous in the scientific community, and that scipy is such a comprehensive toolbox, and that I have not been able to find much about what I'm asking here on the web, leads me to believe I'm on the wrong track or overlooking something. Any help would be appreciated.
鉴于 b 样条在科学界无处不在，而且 scipy 是一个如此全面的工具箱，而且我无法在网上找到很多关于我在这里询问的内容，这让我相信我在错误的轨道或忽视的东西。任何帮助，将不胜感激。
回答by spamduck
I believe scipy's fitpack Library is doing something more complicated than what Mathematica is doing. I was confused as to what was going on as well.
我相信 scipy 的 fitpack 库正在做比 Mathematica 正在做的更复杂的事情。我也很困惑到底发生了什么。
There is the smoothing parameter in these functions, and the default interpolation behavior is to try to make points go through lines. That's what this fitpack software does, so I guess scipy just inherited it? (http://www.netlib.org/fitpack/all-- I'm not sure this is the right fitpack)
这些函数中有平滑参数，默认的插值行为是尝试让点穿过线。这就是 fitpack 软件所做的，所以我猜 scipy 只是继承了它？（http://www.netlib.org/fitpack/all——我不确定这是不是合适的 fitpack）
I took some ideas from http://research.microsoft.com/en-us/um/people/ablake/contours/and coded up your example with the B-splines in there.
我从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()
Anyway I think other folks have the same problems, have a look at: Behavior of scipy's splrep
无论如何，我认为其他人也有同样的问题，请看： scipy's splrep 的行为
回答by Fnord
Use this function i wrote for another question i asked here.
使用我为这里提出的另一个问题编写的这个函数。
In my question i was looking for ways to calculate bsplines with scipy (this is how i actually stumbled upon your question).
在我的问题中，我正在寻找使用 scipy 计算 bsplines 的方法（这就是我实际上偶然发现您的问题的方式）。
After much obsession, i came up with the function below. It'll evaluate any curve up to the 20th degree (way more than we need). And speed wise i tested it for 100,000 samples and it took 0.017s
经过一番痴迷，我想出了下面的功能。它将评估任何高达 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
Results for both open and periodic curves:
开放曲线和周期曲线的结果：
cv = np.array([[ 50., 25.],
[ 59., 12.],
[ 50., 10.],
[ 57., 2.],
[ 40., 4.],
[ 40., 14.]])
