我需要在python中计算b样条函数。为了做到这一点,我写了下面的代码,这是非常好的工作。在
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点很简单:
^{pr2}$现在,为了更快,我可以计算出曲线上所有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
我的速度非常慢的函数在11秒内返回了100k个基本值,但是由于这些值只需要计算一次,所以用这种方法计算曲线上的点的速度比通过splev快6倍。在
事实上,我能够从我的方法和splev得到完全相同的结果,这让我相信splev内部可能会像我一样计算基础,只是速度快得多。在
所以我的目标是找出如何快速计算我的基,把它存储在内存中,然后直接使用美国运输部()来计算曲线上的新点,我的问题是:是否可以使用辛辣。内插为了得到splev用来计算结果的基础值(我猜是这样的)?如果是,怎么办?在
根据unutbu和ev-br关于scipy如何计算样条曲线基的非常有用的见解,我查阅了fortran代码,并尽我所能写了一个等效的代码:
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版本的原始basis函数进行测试:
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_基,因为它相对比较快。但我希望能有改进其性能的建议,并希望能更接近unutbu最初编写的basis函数的版本。在
fitpack_basis
使用一个双循环,它迭代地修改bb
中的元素 和b
。我看不到使用NumPy向量化这些循环的方法,因为 在迭代的每个阶段,bb
和{这是
fitpack_basis
的Cython版本,它的运行速度与 bspline_basis。主要思想 用Cython声明每一个变量的速度 使用普通整数索引将NumPy花式索引的所有用法重写为循环。在请参见this page 关于如何从python构建和运行代码的说明。在
使用这个timeit代码来测试它的性能
^{pr2}$看来Cython可以使
fitpack_basis
代码的运行速度与bspline_basis一样快(也许还快一点):这是
coxDeBoor
的一个版本,它(在我的机器上)比原始版本快840倍。在^{pr2}$
大部分加速是通过使coxDeBoor计算一个 结果而不是一次一个值。请注意,
u
已从 传递给coxDeBoor
的参数。相反,新的coxDeBoor(k, d)
计算 在np.array([coxDeBoor(u[i], k, d) for i in xrange(n)])
之前是什么。在因为NumPy数组算法与标量算法有相同的语法,所以 需要修改的代码很少。唯一的句法变化是在最后 条件:
(u - kv[k] >= 0)
和(u - kv[k + 1] < 0)
是布尔数组。astype
将数组值更改为0和1。所以在单个0或1之前 返回,现在返回一个由0和1组成的完整数组,其中每个值对应一个u
。在Memoization还可以提高性能,因为递归关系 使
coxDeBoor(k, d)
被调用为k
和d
的相同值 不止一次。修饰符语法相当于
并且
memo
修饰符使coxDeBoor
在cache
中记录一个映射 从(k, d)
对参数到返回值。如果coxDeBoor(k, d)
是 再次调用,则返回cache
中的值,而不是 重新计算价值。在scipy_bspline
仍然更快,但至少bspline_basis
加上np.dot
就可以了, 如果您想将b
与多个控制点cv
重复使用,则可能会很有用。在相关问题 更多 >
编程相关推荐