问题描述
在 Scipy 中是否有更有效的方法来生成类似于 b-spline 函数示例中给出的递归代码的 b-spline 基函数:
def B(x,k,i,t):
if k == 0:
return 1.0 if t[i] <= x < t[i+1] else 0.0
if t[i+k] == t[i]:
c1 = 0.0
else:
c1 = (x - t[i])/(t[i+k] - t[i]) * B(x,k-1,t)
if t[i+k+1] == t[i+1]:
c2 = 0.0
else:
c2 = (t[i+k+1] - x)/(t[i+k+1] - t[i+1]) * B(x,i+1,t)
return c1 + c2
我尝试使用 scipy.interpolate.BSpline.basis_element,但无法生成与函数“B”相同的结果。
解决方法
我以前这样做的方法如下:
import scipy.interpolate as si
import numpy as np
import matplotlib.pyplot as plt
nt = 3 # the number of knots
x = np.linspace(0,24,10)
y = np.ones((x.shape)) # the function we wish to interpolate (very trivial in this case)
t = np.linspace(x[0],x[-1],nt)
t = t[1:-1]
tt = np.linspace(0.0,100)
y_rep = si.splrep(x,y,t=t,k=2)
y_i = si.splev(tt,y_rep)
spl = np.zeros((100,))
plt.figure()
for i in range(nt+1):
vec = np.zeros(nt+1)
vec[i] = 1.0
y_list = list(y_rep)
y_list[1] = vec.tolist()
y_i = si.splev(tt,y_list) # your basis spline function
spl = spl + y_i*y_rep[1][i] # the interpolated function
plt.plot(tt,y_i)
您可以查看两个函数 scipy.interpolate.splrep
和 scipy.interpolate.splev
。他们分别进行 B 样条表示并评估 B 样条。上面的代码生成(在这种情况下)四个基本样条函数:
我希望这就是你想要的。