用Python实现样条曲线（使用控制节点和端点）

Question

我正在尝试做类似以下的事情（图片来自维基百科）

#!/usr/bin/env python
from scipy import interpolate
import numpy as np
import matplotlib.pyplot as plt

# sampling
x = np.linspace(0, 10, 10)
y = np.sin(x)

# spline trough all the sampled points
tck = interpolate.splrep(x, y)
x2 = np.linspace(0, 10, 200)
y2 = interpolate.splev(x2, tck)

# spline with all the middle points as knots (not working yet)
# knots = x[1:-1]  # it should be something like this
knots = np.array([x[1]])  # not working with above line and just seeing what this line does
weights = np.concatenate(([1],np.ones(x.shape[0]-2)*.01,[1]))
tck = interpolate.splrep(x, y, t=knots, w=weights)
x3 = np.linspace(0, 10, 200)
y3 = interpolate.splev(x2, tck)

# plot
plt.plot(x, y, 'go', x2, y2, 'b', x3, y3,'r')
plt.show()

这段代码的第一部分是从主要参考文献中提取的代码，但没有解释如何使用点作为控制节点。

该代码的结果是以下图像。

这些点是样本，蓝线是考虑到所有点的样条线。而红线是对我无效的那一条。我正在尝试将所有中间点作为控制节点考虑在内，但我做不到。如果我尝试使用knots=x[1:-1]，它就不起作用。我会感激任何帮助。

简短问题：如何在样条函数中使用所有中间点作为控制节点？

注意：这张最后的图片正是我需要的，它展示了我所拥有的（通过所有点的样条线）和我需要的（带有控制节点的样条线）之间的差异。有什么想法吗？

- silgon

---

哇！你是我在这个社区里见到的第一个使用NX（或Unigraphics）的人。 :-) - fang

你想要做的是通用的B样条实现。如果scipy可以进行B样条插值，那么他们一定已经有了通用的B样条实现。但是是否向客户公开具有适当API的实现尚不确定。 - fang

时间有点紧迫 =)，本来应该完成了，但是没有文档 :/ - silgon

我花了一些时间研究这个问题，但并没有真正取得任何进展。在每种情况下，knots = x[2:-2]都有效。底层Fortran函数的错误代码表明可能违反了Schoenberg-Whitney条件，但我无法理解它。您可以在此处找到实现方式（http://www.netlib.org/dierckx/）。这是curfit.f函数，错误代码为`ier = 10`。至少Fortran函数有很好的文档说明。 - spfrnd

Matplotlib 知道贝塞尔曲线：http://matplotlib.org/users/path_tutorial.html#bezier-example - Dietrich

Answer 1

如果您想要评估 B样条，您需要确定适当的节点向量，然后手动重建tck以满足您的需求。 tck代表节点t+系数c+曲线度数k。 splrep为通过给定控制点的三次曲线计算tck。因此，您不能将其用于您想要的内容。

下面的函数将展示我对我之前提出的一个类似问题的解决方案。进行了调整以适应您的需求。

有趣的事实：该代码适用于任何维度的曲线（1D、2D、3D、...、nD）。

import numpy as np
import scipy.interpolate as si


def bspline(cv, n=100, degree=3):
    """ Calculate n samples on a bspline

        cv :      Array ov control vertices
        n  :      Number of samples to return
        degree:   Curve degree
    """
    cv = np.asarray(cv)
    count = cv.shape[0]

    # Prevent degree from exceeding count-1, otherwise splev will crash
    degree = np.clip(degree,1,count-1)

    # Calculate knot vector
    kv = np.array([0]*degree + list(range(count-degree+1)) + [count-degree]*degree,dtype='int')

    # Calculate query range
    u = np.linspace(0,(count-degree),n)

    # Calculate result
    return np.array(si.splev(u, (kv,cv.T,degree))).T

测试一下：

import matplotlib.pyplot as plt
colors = ('b', 'g', 'r', 'c', 'm', 'y', 'k')

cv = np.array([[ 50.,  25.],
   [ 59.,  12.],
   [ 50.,  10.],
   [ 57.,   2.],
   [ 40.,   4.],
   [ 40.,   14.]])

plt.plot(cv[:,0],cv[:,1], 'o-', label='Control Points')

for d in range(1,5):
    p = bspline(cv,n=100,degree=d)
    x,y = p.T
    plt.plot(x,y,'k-',label='Degree %s'%d,color=colors[d%len(colors)])

plt.minorticks_on()
plt.legend()
plt.xlabel('x')
plt.ylabel('y')
plt.xlim(35, 70)
plt.ylim(0, 30)
plt.gca().set_aspect('equal', adjustable='box')
plt.show()

结果：

- Fnord

---

TypeError：bspline() 函数传入了一个意外的关键字参数 'periodic'， TypeError：只能将列表（而不是“range”）与列表相加。 - kimstik

@kimstik 感谢您指出这一点。 - Fnord

Answer 2

在这个IPython笔记本http://nbviewer.ipython.org/github/empet/geom_modeling/blob/master/FP-Bezier-Bspline.ipynb中，您可以找到生成B样条曲线所涉及的数据的详细描述，以及de Boor算法的Python实现。

- xecafe

---

+1，甚至可以加10分！看起来非常不错，我周末一定会去试试（在此之前无法尝试）。非常感谢。;) - silgon

Answer 3

我在这个链接中发现了一些关于贝塞尔曲线的有趣答案，它提供了我所需要的答案。然后我使用了这段代码来尝试自己做一个。看起来它运行得很好。这是我的实现：

#! /usr/bin/python
# -*- coding: utf-8 -*-
import numpy as np
import matplotlib.pyplot as plt
from scipy.special import binom

def Bernstein(n, k):
    """Bernstein polynomial.

    """
    coeff = binom(n, k)

    def _bpoly(x):
        return coeff * x ** k * (1 - x) ** (n - k)

    return _bpoly


def Bezier(points, num=200):
    """Build Bézier curve from points.

    """
    N = len(points)
    t = np.linspace(0, 1, num=num)
    curve = np.zeros((num, 2))
    for ii in range(N):
        curve += np.outer(Bernstein(N - 1, ii)(t), points[ii])
    return curve
xp = np.array([2,3,4,5])
yp = np.array([2,1,4,0])
x, y = Bezier(list(zip(xp, yp))).T

plt.plot(x,y)
plt.plot(xp,yp,"ro")
plt.plot(xp,yp,"b--")

plt.show()

以下为例子的图片。

红点代表控制点。 就是这样 =)

- silgon

---

你的链接无效。已更新为使用Python进行交互式贝塞尔曲线。 - Bijan

谢谢@Bijan，我已经在帖子中更新了它=) - silgon

@silgon我们能否打印拟合曲线的多项式函数？ - Natasha

Answer 4

我认为问题与你的节点向量有关。如果选择太多的节点，它似乎会引起问题，需要在节点之间有一些数据点。这个问题解决了 scipy.interpolate的splrep函数选择节点时是否存在Bug？

#!/usr/bin/env python
from scipy import interpolate
import numpy as np
import matplotlib.pyplot as plt

# sampling
x = np.linspace(0, 10, 10)
y = np.sin(x)

# spline trough all the sampled points
tck = interpolate.splrep(x, y)
print tck
x2 = np.linspace(0, 10, 200)
y2 = interpolate.splev(x2, tck)

# spline with all the middle points as knots (not working yet)
knots = np.asarray(x[1:-1])  # it should be something like this
#knots = np.array([x[1]])  # not working with above line and just seeing what this line does
nknots = 5
idx_knots = (np.arange(1,len(x)-1,(len(x)-2)/np.double(nknots))).astype('int')
knots = x[idx_knots]
print knots

weights = np.concatenate(([1],np.ones(x.shape[0]-2)*.01,[1]))
tck = interpolate.splrep(x, y,  t=knots, w=weights)
x3 = np.linspace(0, 10, 200)
y3 = interpolate.splev(x2, tck)

# plot
plt.plot(x, y, 'go', x2, y2, 'b', x3, y3,'r')
plt.show()

选择5节时似乎能够正常工作，6节会产生奇怪的结果，再多就会出现错误。

- Salix alba

---

抱歉之前没有留言，我会在明天一早看一下。 - silgon

我刚刚遇到了一个丑陋的错误http://pastebin.com/GNVGLsmF，顺便说一下，有一个括号前面多了一个1，它不应该在那里=) - silgon

我刚刚更新了示例。我已经剪切并粘贴了在我的机器上运行的代码。 - Salix alba

Answer 5

你的示例函数是周期性的，需要在 interpolate.splrep 方法中添加 per=True 选项。

knots = x[1:-1]
weights = np.concatenate(([1],np.ones(x.shape[0]-2)*.01,[1]))
tck = interpolate.splrep(x, y, t=knots, w=weights, per=True)

这给了我以下内容：

编辑：这也解释了为什么使用knots = x[-2:2]，它是完整范围的非周期子集，也能够正常工作。

- John

---

抱歉我还没有检查它，我明天早上会试一下。不过从我看到的来看，样条曲线的最后一个点并不是真正的结束点，这在应用程序中是非常重要的一件事情。 - silgon

刚在我的电脑上测试了一下，它没有将数据的端点视为样条的端点 :/，所以... - silgon