Python 的样条曲线（使用控制节点和端点）
Splines with Python (using control knots and endpoints)
我正在尝试做类似下面的事情（图片从维基百科中提取）
#!/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()
代码的第一部分是从the main reference中提取的代码，但没有解释如何使用这些点作为控制结。
这段代码的结果是下图。
点是样本，蓝线是考虑了所有点的样条。红线是对我不起作用的那条。我正在尝试将 所有中间点都作为控制节点 考虑在内，但我做不到。如果我尝试使用 knots=x[1:-1] 它就是行不通的。如果有任何帮助，我将不胜感激。
简而言之问题：如何在样条函数中使用所有中间点作为控制节点？
注意：这最后一张图片正是我所需要的，它是我所拥有的（样条曲线通过所有点）和我所需要的（带控制结的样条曲线）之间的区别。有任何想法吗？
您的示例函数是周期性的，您需要将 per=True 选项添加到 interpolate.splrep 方法。
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]，它是完整范围的 non-periodic 子集。
我认为问题出在你的节点向量上。如果你 select 太多结似乎会导致问题，它需要在结之间有一些数据点。本题解决问题Bug (?) on selecting knots on scipy.insterpolate's splrep function
#!/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 节给出奇怪的结果，再多则给出错误。
我刚刚在 this link 中发现了一些非常有趣的答案，我需要一个贝塞尔曲线。然后我用代码自己试了一下。它显然工作正常。这是我的实现：
#! /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()
还有一张图片作为示例。
红色的点代表控制点。
就是这样 =)
在此 IPython Notebook http://nbviewer.ipython.org/github/empet/geom_modeling/blob/master/FP-Bezier-Bspline.ipynb 中，您可以找到生成 B 样条曲线所涉及数据的详细描述，以及 Python de Boor 算法的实现.
如果您想要计算 bspline，您需要为您的样条找出合适的节点向量，然后手动重建 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 + 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,periodic=True)
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()
结果：
我正在尝试做类似下面的事情（图片从维基百科中提取）
#!/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()
代码的第一部分是从the main reference中提取的代码，但没有解释如何使用这些点作为控制结。
这段代码的结果是下图。
点是样本，蓝线是考虑了所有点的样条。红线是对我不起作用的那条。我正在尝试将 所有中间点都作为控制节点 考虑在内，但我做不到。如果我尝试使用 knots=x[1:-1] 它就是行不通的。如果有任何帮助，我将不胜感激。
简而言之问题：如何在样条函数中使用所有中间点作为控制节点？
注意：这最后一张图片正是我所需要的，它是我所拥有的（样条曲线通过所有点）和我所需要的（带控制结的样条曲线）之间的区别。有任何想法吗？
您的示例函数是周期性的，您需要将 per=True 选项添加到 interpolate.splrep 方法。
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]，它是完整范围的 non-periodic 子集。
我认为问题出在你的节点向量上。如果你 select 太多结似乎会导致问题，它需要在结之间有一些数据点。本题解决问题Bug (?) on selecting knots on scipy.insterpolate's splrep function
#!/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 节给出奇怪的结果，再多则给出错误。
我刚刚在 this link 中发现了一些非常有趣的答案，我需要一个贝塞尔曲线。然后我用代码自己试了一下。它显然工作正常。这是我的实现：
#! /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()
还有一张图片作为示例。
红色的点代表控制点。 就是这样 =)
在此 IPython Notebook http://nbviewer.ipython.org/github/empet/geom_modeling/blob/master/FP-Bezier-Bspline.ipynb 中，您可以找到生成 B 样条曲线所涉及数据的详细描述，以及 Python de Boor 算法的实现.
如果您想要计算 bspline，您需要为您的样条找出合适的节点向量，然后手动重建 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 + 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,periodic=True)
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()
结果：