Мне нужно оценить b-сплайны в python. Для этого я написал код, который работает очень хорошо.Как получить сплайн-основу, используемую scipy.interpolate.splev
import numpy as np
import scipy.interpolate as si
def scipy_bspline(cv,n,degree):
""" bspline basis function
c = list of control points.
n = number of points on the curve.
degree = curve degree
"""
# Create a range of u values
c = cv.shape[0]
kv = np.clip(np.arange(c+degree+1)-degree,0,c-degree)
u = np.linspace(0,c-degree,n)
# Calculate result
return np.array(si.splev(u, (kv,cv.T,degree))).T
Придав ему 6 контрольных точек и попросить его оценить 100k точек на кривой ветер:
# Control points
cv = np.array([[ 50., 25., 0.],
[ 59., 12., 0.],
[ 50., 10., 0.],
[ 57., 2., 0.],
[ 40., 4., 0.],
[ 40., 14., 0.]])
n = 100000 # 100k Points
degree = 3 # Curve degree
points_scipy = scipy_bspline(cv,n,degree) #cProfile clocks this at 0.012 seconds
Вот сюжет "points_scipy": 
Теперь, чтобы сделать это быстрее я мог бы вычислить основу для всех 100k точек на кривой, сохранить это в памяти, и когда мне нужно нарисовать кривую, все, что мне нужно будет сделать, - это умножить новые позиции контрольной точки с сохраненной основой, чтобы получить новую кривую , Для того, чтобы доказать свою точку зрения я написал функцию, которая использует DeBoor's algorithm вычислить свою базу:
def basis(c, n, degree):
""" bspline basis function
c = number of control points.
n = number of points on the curve.
degree = curve degree
"""
# Create knot vector and a range of samples on the curve
kv = np.array([0]*degree + range(c-degree+1) + [c-degree]*degree,dtype='int') # knot vector
u = np.linspace(0,c-degree,n) # samples range
# Cox - DeBoor recursive function to calculate basis
def coxDeBoor(u, k, d):
# Test for end conditions
if (d == 0):
if (kv[k] <= u and u < kv[k+1]):
return 1
return 0
Den1 = kv[k+d] - kv[k]
Den2 = 0
Eq1 = 0
Eq2 = 0
if Den1 > 0:
Eq1 = ((u-kv[k])/Den1) * coxDeBoor(u,k,(d-1))
try:
Den2 = kv[k+d+1] - kv[k+1]
if Den2 > 0:
Eq2 = ((kv[k+d+1]-u)/Den2) * coxDeBoor(u,(k+1),(d-1))
except:
pass
return Eq1 + Eq2
# Compute basis for each point
b = np.zeros((n,c))
for i in xrange(n):
for k in xrange(c):
b[i][k%c] += coxDeBoor(u[i],k,degree)
b[n-1][-1] = 1
return b
Теперь позволяет использовать это, чтобы вычислить новую основу, умножаем на контрольных точках и подтверждают, что мы получаем те же результаты, splev:
b = basis(len(cv),n,degree) #5600011 function calls (600011 primitive calls) in 10.975 seconds
points_basis = np.dot(b,cv) #3 function calls in 0.002 seconds
print np.allclose(points_basis,points_scipy) # Returns True
Моя очень медленно функция вернула 100k базисные значения в 11 секунд, но так как эти значения должны только быть вычислен один раз, вычисления точки на кривой заканчивал тем, что в 6 раз быстрее, таким образом, чем делать это через splev ,
Тот факт, что я смог получить точные результаты как от моего метода, так и от splev, заставляет меня поверить, что внутренне splev, вероятно, вычисляет базу, как я, за исключением гораздо более быстрого.
Итак, моя цель - выяснить, как быстро вычислить мой базис, сохранить его в памяти и просто использовать np.dot() для вычисления новых точек на кривой, и мой вопрос: можно ли использовать пряный. интерполировать, чтобы получить базовые значения, которые (я полагаю) splev использует для вычисления его результата? И если да, то как?
[ДОБАВЛЕНИЕ]
После unutbu годов и очень полезное понимание ЭВ-БР о том, как SciPy вычисляет основы сплайна, я посмотрел Фортран код и написал эквивалент в меру моих способностей:
def fitpack_basis(c, n=100, d=3, rMinOffset=0, rMaxOffset=0):
""" fitpack's spline basis function
c = number of control points.
n = number of points on the curve.
d = curve degree
"""
# Create knot vector
kv = np.array([0]*d + range(c-d+1) + [c-d]*d, dtype='int')
# Create sample range
u = np.linspace(rMinOffset, rMaxOffset + c - d, n) # samples range
# Create buffers
b = np.zeros((n,c)) # basis
bb = np.zeros((n,c)) # basis buffer
left = np.clip(np.floor(u),0,c-d-1).astype(int) # left knot vector indices
right = left+d+1 # right knot vector indices
# Go!
nrange = np.arange(n)
b[nrange,left] = 1.0
for j in xrange(1, d+1):
crange = np.arange(j)[:,None]
bb[nrange,left+crange] = b[nrange,left+crange]
b[nrange,left] = 0.0
for i in xrange(j):
f = bb[nrange,left+i]/(kv[right+i] - kv[right+i-j])
b[nrange,left+i] = b[nrange,left+i] + f * (kv[right+i] - u)
b[nrange,left+i+1] = f * (u - kv[right+i-j])
return b
Тестирование против версии unutbu в моей исходной базисной функции:
fb = fitpack_basis(c,n,d) #22 function calls in 0.044 seconds
b = basis(c,n,d) #81 function calls (45 primitive calls) in 0.013 seconds ~5 times faster
print np.allclose(b,fb) # Returns True
Моя функция в 5 раз медленнее, но все еще относительно быстро. Что мне нравится в этом, так это то, что он позволяет использовать диапазоны образцов, выходящие за границы, что является чем-то полезным в моем приложении. Например:
print fitpack_basis(c,5,d,rMinOffset=-0.1,rMaxOffset=.2)
[[ 1.331 -0.3468 0.0159 -0.0002 0. 0. ]
[ 0.0208 0.4766 0.4391 0.0635 0. 0. ]
[ 0. 0.0228 0.4398 0.4959 0.0416 0. ]
[ 0. 0. 0.0407 0.3621 0.5444 0.0527]
[ 0. 0. -0.0013 0.0673 -0.794 1.728 ]]
Поэтому по этой причине я, вероятно, буду использовать fitpack_basis, так как это относительно быстро. Но мне бы хотелось, чтобы предложения по улучшению его производительности и, надеюсь, приблизились к версии оригинальной исходной функции unutbu, которую я написал.
Пожалуйста, разместите код, который вы использовали для расчета основы для каждого образца. – unutbu
Что именно вы подразумеваете под «на выборку, используемой splev при оценке сплайна»? Вы поняли, как дать splev узлы и коэффициенты, что именно вы ищете? –
@unutbu я полностью переписал свой вопрос + добавленный код (медленно) получить базу – Fnord