如何获得scipy.interpolate.splev使用的样条线基础 [英] How to get the spline basis used by scipy.interpolate.splev
问题描述
我需要在python中评估b样条曲线.为此，我在下面编写了很好的代码.
I need to evaluate b-splines in python. To do so i wrote the code below which works very well.
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个控制点并要求它评估曲线上的10万个点很容易:
Giving it 6 control points and asking it to evaluate 100k points on curve is a breeze:
# 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"的图:
Here's a plot of "points_scipy":
现在，要使其更快，我可以计算曲线上所有100k点的基础，并将其存储在内存中，当我需要绘制曲线时，我要做的就是将新的控制点位置乘以存储基础以获取新曲线.为了证明我的观点，我编写了一个函数，该函数使用 DeBoor算法来计算我的依据:
Now, to make this faster i could compute a basis for all 100k points on the curve, store that in memory, and when i need to draw a curve all i would need to do is multiply the new control point positions with the stored basis to get the new curve. To prove my point i wrote a function that uses DeBoor's algorithm to compute my basis:
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相同:
Now lets use this to compute a new basis, multiply that by the control points and confirm that we're getting the same results as 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
我的极慢函数在11秒内返回了100k个基值，但是由于这些值只需要计算一次，因此，计算曲线上的点最终比通过splev进行计算快6倍.
My extremely slow function returned 100k basis values in 11 seconds, but since these values only need to be computed once, calculating the points on curve ended up being 6 times faster this way than doing it via splev.
我能够从我的方法和splev中获得完全相同的结果，这使我相信内部splev可能会像我一样计算基数，只是速度要快得多.
The fact that I was able to get the exact same results from both my method and splev leads me to believe that internally splev probably calculates a basis like i do, except much faster.
所以我的目标是找出如何快速计算我的基础，将其存储在内存中，然后仅使用np.dot()来计算曲线上的新点，我的问题是:是否可以使用辣.获取(我认为)splev用于计算结果的基值?如果可以，怎么办?
So my goal is to find out how to calculate my basis fast, store it in memory and just use np.dot() to compute the new points on curve, and my question is: Is it possible to use spicy.interpolate to get the basis values that (i presume) splev uses to compute it's result? And if so, how?
在unutbu和ev-br关于scipy如何计算样条曲线的基础的非常有用的见解之后，我查找了fortran代码并编写了与我所能达到的最佳功能等效的代码:
Following unutbu's and ev-br's very useful insight on how scipy calculates a spline's basis, I looked up the fortran code and wrote a equivalent to the best of my abilities:
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版本进行测试:
Testing against unutbu's version of my original basis function:
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倍，但仍然相对快.我喜欢它的原因是它使我可以使用超出范围的样本范围，这在我的应用程序中很有用.例如:
My function is 5 times slower, but still relatively fast. What I like about it is that it lets me use sample ranges that go beyond the boundaries, which is something of use in my application. For example:
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版本的原始基函数更接近.
So for that reason I will probably use fitpack_basis since it's relatively fast. But i would love suggestions for improving its performance and hopefully get closer to unutbu's version of the original basis function i wrote.
推荐答案
fitpack_basis使用double循环迭代地修改bb中的元素
和b.我看不到使用NumPy向量化这些循环的方法，因为该值
迭代的每个阶段的bb和b的值取决于
以前的迭代.在这种情况下，Cython有时可以用来
改善循环的性能.
fitpack_basis uses a double loop which iteratively modifies elements in bb
and b. I don't see a way to use NumPy to vectorize these loops since the value
of bb and b at each stage of the iteration depends on the values from
previous iterations. In situations like this sometimes Cython can be used to
improve the performance of the loops.
这里是fitpack_basis的Cython版本，运行速度与
bspline_basis .主要思想
使用Cython提升速度的方法是声明每个变量的类型，并且
使用纯整数索引将NumPy花式索引的所有用法重写为循环.
Here is a Cython-ized version of fitpack_basis, which runs as fast as
bspline_basis. The main ideas
used to boost speed using Cython is to declare the type of every variable, and
rewrite all uses of NumPy fancy indexing as loops using plain integer indexing.
请参见此 页面 有关如何构建代码并从python运行代码的说明.
See this page for instructions on how to build the code and run it from python.
import numpy as np
cimport numpy as np
cimport cython
ctypedef np.float64_t DTYPE_f
ctypedef np.int64_t DTYPE_i
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.nonecheck(False)
def cython_fitpack_basis(int c, int n=100, int d=3,
double rMinOffset=0, double rMaxOffset=0):
""" fitpack's spline basis function
c = number of control points.
n = number of points on the curve.
d = curve degree
"""
cdef Py_ssize_t i, j, k, l
cdef double f
# Create knot vector
cdef np.ndarray[DTYPE_i, ndim=1] kv = np.array(
[0]*d + range(c-d+1) + [c-d]*d, dtype=np.int64)
# Create sample range
cdef np.ndarray[DTYPE_f, ndim=1] u = np.linspace(
rMinOffset, rMaxOffset + c - d, n)
# basis
cdef np.ndarray[DTYPE_f, ndim=2] b = np.zeros((n,c))
# basis buffer
cdef np.ndarray[DTYPE_f, ndim=2] bb = np.zeros((n,c))
# left knot vector indices
cdef np.ndarray[DTYPE_i, ndim=1] left = np.clip(np.floor(u), 0, c-d-1).astype(np.int64)
# right knot vector indices
cdef np.ndarray[DTYPE_i, ndim=1] right = left+d+1
for k in range(n):
b[k, left[k]] = 1.0
for j in range(1, d+1):
for l in range(j):
for k in range(n):
bb[k, left[k] + l] = b[k, left[k] + l]
b[k, left[k]] = 0.0
for i in range(j):
for k in range(n):
f = bb[k, left[k]+i] / (kv[right[k]+i] - kv[right[k]+i-j])
b[k, left[k]+i] = b[k, left[k]+i] + f * (kv[right[k]+i] - u[k])
b[k, left[k]+i+1] = f * (u[k] - kv[right[k]+i-j])
return b
使用此timeit代码对性能进行基准测试，
Using this timeit code to benchmark it's performance,
import timeit
import numpy as np
import cython_bspline as CB
import numpy_bspline as NB
c = 6
n = 10**5
d = 3
fb = NB.fitpack_basis(c, n, d)
bb = NB.bspline_basis(c, n, d)
cfb = CB.cython_fitpack_basis(c,n,d)
assert np.allclose(bb, fb)
assert np.allclose(cfb, fb)
# print(NB.fitpack_basis(c,5,d,rMinOffset=-0.1,rMaxOffset=.2))
timing = dict()
timing['NB.fitpack_basis'] = timeit.timeit(
stmt='NB.fitpack_basis(c, n, d)',
setup='from __main__ import NB, c, n, d',
number=10)
timing['NB.bspline_basis'] = timeit.timeit(
stmt='NB.bspline_basis(c, n, d)',
setup='from __main__ import NB, c, n, d',
number=10)
timing['CB.cython_fitpack_basis'] = timeit.timeit(
stmt='CB.cython_fitpack_basis(c, n, d)',
setup='from __main__ import CB, c, n, d',
number=10)
for func_name, t in timing.items():
print "{:>25}: {:.4f}".format(func_name, t)
看来，Cython可以使fitpack_basis代码的运行速度与 bspline_basis 一样快(可能更快). :
it appears Cython can make the fitpack_basis code run as fast as (and perhaps a bit faster) than bspline_basis:
NB.bspline_basis: 0.3322
CB.cython_fitpack_basis: 0.2939
NB.fitpack_basis: 0.9182
这篇关于如何获得scipy.interpolate.splev使用的样条线基础的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
