在 scipy文档 中没有提到如何手动生成样条曲线的系数。然而,可以从现有的B样条文献中找到相关内容。以下函数bspleval
展示了如何构建B样条基函数(代码中的矩阵B
),通过将系数与最高阶基函数相乘并求和,可以轻松地生成样条曲线:
def bspleval(x, knots, coeffs, order, debug=False):
'''
Evaluate a B-spline at a set of points.
Parameters
----------
x : list or ndarray
The set of points at which to evaluate the spline.
knots : list or ndarray
The set of knots used to define the spline.
coeffs : list of ndarray
The set of spline coefficients.
order : int
The order of the spline.
Returns
-------
y : ndarray
The value of the spline at each point in x.
'''
k = order
t = knots
m = alen(t)
npts = alen(x)
B = zeros((m-1,k+1,npts))
if debug:
print('k=%i, m=%i, npts=%i' % (k, m, npts))
print('t=', t)
print('coeffs=', coeffs)
for i in range(m-1):
B[i,0,:] = float64(logical_and(x >= t[i], x < t[i+1]))
if (k == 0):
B[m-2,0,-1] = 1.0
for j in range(1,k+1):
for i in range(m-j-1):
if (t[i+j] - t[i] == 0.0):
first_term = 0.0
else:
first_term = ((x - t[i]) / (t[i+j] - t[i])) * B[i,j-1,:]
if (t[i+j+1] - t[i+1] == 0.0):
second_term = 0.0
else:
second_term = ((t[i+j+1] - x) / (t[i+j+1] - t[i+1])) * B[i+1,j-1,:]
B[i,j,:] = first_term + second_term
B[m-j-2,j,-1] = 1.0
if debug:
plt.figure()
for i in range(m-1):
plt.plot(x, B[i,k,:])
plt.title('B-spline basis functions')
y = zeros(npts)
for i in range(m-k-1):
y += coeffs[i] * B[i,k,:]
if debug:
plt.figure()
plt.plot(x, y)
plt.title('spline curve')
plt.show()
return(y)
为了举例说明如何将此方法与Scipy现有的单变量样条函数结合使用,以下是一个示例脚本。它采用输入数据并使用Scipy的函数和面向对象的方法进行样条拟合。从这两个方法中获取系数和结点,并将其作为输入传递给我们手动计算的
bspleval
程序,我们可以复制它们生成的相同曲线。请注意,手动评估曲线和Scipy评估方法之间的差异非常小,几乎肯定是浮点噪声。
x = array([-273.0, -176.4, -79.8, 16.9, 113.5, 210.1, 306.8, 403.4, 500.0])
y = array([2.25927498e-53, 2.56028619e-03, 8.64512988e-01, 6.27456769e+00, 1.73894734e+01,
3.29052124e+01, 5.14612316e+01, 7.20531200e+01, 9.40718450e+01])
x_nodes = array([-273.0, -263.5, -234.8, -187.1, -120.3, -34.4, 70.6, 194.6, 337.8, 500.0])
y_nodes = array([2.25927498e-53, 3.83520726e-46, 8.46685318e-11, 6.10568083e-04, 1.82380809e-01,
2.66344008e+00, 1.18164677e+01, 3.01811501e+01, 5.78812583e+01, 9.40718450e+01])
k = 3
tck = splrep(x_nodes, y_nodes, k=k, s=0)
knots = tck[0]
coeffs = tck[1]
print('knot points=', knots)
print('coefficients=', coeffs)
uspline = UnivariateSpline(x_nodes, y_nodes, s=0)
uspline_knots = uspline.get_knots()
uspline_coeffs = uspline.get_coeffs()
ytck = splev(x, tck)
y_uspline = uspline(x)
y_knots = uspline(knots)
y_eval = bspleval(x, knots, coeffs, k, debug=False)
plt.plot(x, ytck, label='tck')
plt.plot(x, y_uspline, label='uspline')
plt.plot(x, y_eval, label='manual')
plt.plot(x_nodes, y_nodes, 'ko', markersize=7, label='input nodes')
plt.plot(knots, y_knots, 'mo', markersize=5, label='tck knots')
plt.xlim((-300.0,530.0))
plt.legend(loc='best', prop={'size':14})
plt.figure()
plt.title('difference')
plt.plot(x, ytck-y_uspline, label='tck-uspl')
plt.plot(x, ytck-y_eval, label='tck-manual')
plt.legend(loc='best', prop={'size':14})
plt.show()
![enter image description here](https://istack.dev59.com/LuRmE.webp)
x = np.linspace(0, 100, 1000)
,d = np.sin(x * 0.5) + 2 + np.cos(x * 0.1)
,x_sample = x[::50]
,d_sample = d[::50]
,s = UnivariateSpline(x_sample, d_sample, k=3, s=0.005)
,spl = bspleval(x, s.get_knots(), s.get_coeffs(), 3, debug=False)
。生成的图表显示d
和spl
看起来相当不同。有什么想法吗? - Clebx_sample
和d_sample
作为样条拟合的输入,但是你将其与x
和d
进行比较。如果你让x_sample=x
和d_sample=d
,那么你会得到更好的拟合效果。但仍然存在误差。接下来,如果你将平滑因子从s=0.005
降低到非常小的值(比如s=1.0e-15
),那么你会得到一个很好的拟合效果。 - nzh