如何获得 scipy.interpolate.splev 使用的样条基
How to get the spline basis used by scipy.interpolate.splev
我需要评估 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 个点是一件轻而易举的事：
# 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
我的极慢函数在 11 秒内返回了 100k 基值，但由于这些值只需要计算一次，因此计算曲线上的点最终比通过 splev 快 6 倍。
我能够从我的方法和 splev 中得到完全相同的结果这一事实让我相信内部 splev 可能会像我一样计算一个基础，只是速度要快得多。
所以我的目标是找出如何快速计算我的基础，将其存储在内存中并仅使用 np.dot() 来计算曲线上的新点，我的问题是：是否有可能使用 spicy.interpolate 获取（我认为）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 版本进行测试：
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 版本。
这是 coxDeBoor 的一个版本，（在我的机器上）比原始版本快 840 倍。
In [102]: %timeit basis(len(cv), 10**5, degree)
1 loops, best of 3: 24.5 s per loop
In [104]: %timeit bspline_basis(len(cv), 10**5, degree)
10 loops, best of 3: 29.1 ms per loop
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 = len(cv)
kv = np.array(
[0] * degree + range(c - degree + 1) + [c - degree] * degree, dtype='int')
u = np.linspace(0, c - degree, n)
# Calculate result
arange = np.arange(n)
points = np.zeros((n, cv.shape[1]))
for i in xrange(cv.shape[1]):
points[arange, i] = si.splev(u, (kv, cv[:, i], degree))
return points
def memo(f):
# Peter Norvig's
"""Memoize the return value for each call to f(args).
Then when called again with same args, we can just look it up."""
cache = {}
def _f(*args):
try:
return cache[args]
except KeyError:
cache[args] = result = f(*args)
return result
except TypeError:
# some element of args can't be a dict key
return f(*args)
_f.cache = cache
return _f
def bspline_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
@memo
def coxDeBoor(k, d):
# Test for end conditions
if (d == 0):
return ((u - kv[k] >= 0) & (u - kv[k + 1] < 0)).astype(int)
denom1 = kv[k + d] - kv[k]
term1 = 0
if denom1 > 0:
term1 = ((u - kv[k]) / denom1) * coxDeBoor(k, d - 1)
denom2 = kv[k + d + 1] - kv[k + 1]
term2 = 0
if denom2 > 0:
term2 = ((-(u - kv[k + d + 1]) / denom2) * coxDeBoor(k + 1, d - 1))
return term1 + term2
# Compute basis for each point
b = np.column_stack([coxDeBoor(k, degree) for k in xrange(c)])
b[n - 1][-1] = 1
return b
# Control points
cv = np.array([[50., 25., 0.],
[59., 12., 0.],
[50., 10., 0.],
[57., 2., 0.],
[40., 4., 0.],
[40., 14., 0.]])
n = 10 ** 6
degree = 3 # Curve degree
points_scipy = scipy_bspline(cv, n, degree)
b = bspline_basis(len(cv), n, degree)
points_basis = np.dot(b, cv)
print np.allclose(points_basis, points_scipy)
大部分加速是通过让 coxDeBoor 计算一个向量来实现的
结果而不是一次单个值。请注意，u 已从
传递给 coxDeBoor 的参数。相反，新的 coxDeBoor(k, d) 计算
np.array([coxDeBoor(u[i], k, d) for i in xrange(n)]).
之前是什么
由于 NumPy 数组运算与标量运算具有相同的语法，因此非常
需要更改的代码很少。唯一的句法变化是在最后
条件：
if (d == 0):
return ((u - kv[k] >= 0) & (u - kv[k + 1] < 0)).astype(int)
(u - kv[k] >= 0) 和 (u - kv[k + 1] < 0) 是布尔数组。 astype
将数组值更改为 0 和 1。因此，在单个 0 或 1 之前
返回，现在返回一整组 0 和 1——每个值对应一个
u.
Memoization也可以提高性能，因为递归关系
导致 coxDeBoor(k, d) 被调用以获得 k 和 d 的相同值
不止一次。装饰器语法
@memo
def coxDeBoor(k, d):
...
等同于
def coxDeBoor(k, d):
...
coxDeBoor = memo(coxDeBoor)
和 memo 装饰器导致 coxDeBoor 在 cache 映射中记录
从 (k, d) 对参数到返回值。如果 coxDeBoor(k, d) 是
再次调用，则返回 cache 中的值而不是
重新计算值。
scipy_bspline 仍然更快，但至少 bspline_basis 加上 np.dot 在大概范围内，
如果你想重复使用 b 和许多控制点，cv.
可能会很有用
In [109]: %timeit scipy_bspline(cv, 10**5, degree)
100 loops, best of 3: 17.2 ms per loop
In [104]: %timeit b = bspline_basis(len(cv), 10**5, degree)
10 loops, best of 3: 29.1 ms per loop
In [111]: %timeit points_basis = np.dot(b, cv)
100 loops, best of 3: 2.46 ms per loop
fitpack_basis 使用双循环迭代修改 bb 中的元素
和 b。我看不到使用 NumPy 来矢量化这些循环的方法，因为值
bb 和 b 在每个迭代阶段的值取决于来自
以前的迭代。在这种情况下，有时可以使用 Cython
提高循环的性能。
这是 fitpack_basis 的 Cython 化版本，运行 的速度与
。主要思想
用于提高使用 Cython 的速度的是声明每个变量的类型，并且
使用普通整数索引将 NumPy 花式索引的所有使用重写为循环。
参见这个
第页为
关于如何构建代码和 运行 它的说明来自 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 代码对其性能进行基准测试，
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 代码 运行 与 一样快（也许快一点）：
NB.bspline_basis: 0.3322
CB.cython_fitpack_basis: 0.2939
NB.fitpack_basis: 0.9182
我需要评估 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 个点是一件轻而易举的事：
# 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
我的极慢函数在 11 秒内返回了 100k 基值，但由于这些值只需要计算一次，因此计算曲线上的点最终比通过 splev 快 6 倍。
我能够从我的方法和 splev 中得到完全相同的结果这一事实让我相信内部 splev 可能会像我一样计算一个基础，只是速度要快得多。
所以我的目标是找出如何快速计算我的基础，将其存储在内存中并仅使用 np.dot() 来计算曲线上的新点，我的问题是：是否有可能使用 spicy.interpolate 获取（我认为）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 版本进行测试：
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 版本。
这是 coxDeBoor 的一个版本，（在我的机器上）比原始版本快 840 倍。
In [102]: %timeit basis(len(cv), 10**5, degree)
1 loops, best of 3: 24.5 s per loop
In [104]: %timeit bspline_basis(len(cv), 10**5, degree)
10 loops, best of 3: 29.1 ms per loop
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 = len(cv)
kv = np.array(
[0] * degree + range(c - degree + 1) + [c - degree] * degree, dtype='int')
u = np.linspace(0, c - degree, n)
# Calculate result
arange = np.arange(n)
points = np.zeros((n, cv.shape[1]))
for i in xrange(cv.shape[1]):
points[arange, i] = si.splev(u, (kv, cv[:, i], degree))
return points
def memo(f):
# Peter Norvig's
"""Memoize the return value for each call to f(args).
Then when called again with same args, we can just look it up."""
cache = {}
def _f(*args):
try:
return cache[args]
except KeyError:
cache[args] = result = f(*args)
return result
except TypeError:
# some element of args can't be a dict key
return f(*args)
_f.cache = cache
return _f
def bspline_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
@memo
def coxDeBoor(k, d):
# Test for end conditions
if (d == 0):
return ((u - kv[k] >= 0) & (u - kv[k + 1] < 0)).astype(int)
denom1 = kv[k + d] - kv[k]
term1 = 0
if denom1 > 0:
term1 = ((u - kv[k]) / denom1) * coxDeBoor(k, d - 1)
denom2 = kv[k + d + 1] - kv[k + 1]
term2 = 0
if denom2 > 0:
term2 = ((-(u - kv[k + d + 1]) / denom2) * coxDeBoor(k + 1, d - 1))
return term1 + term2
# Compute basis for each point
b = np.column_stack([coxDeBoor(k, degree) for k in xrange(c)])
b[n - 1][-1] = 1
return b
# Control points
cv = np.array([[50., 25., 0.],
[59., 12., 0.],
[50., 10., 0.],
[57., 2., 0.],
[40., 4., 0.],
[40., 14., 0.]])
n = 10 ** 6
degree = 3 # Curve degree
points_scipy = scipy_bspline(cv, n, degree)
b = bspline_basis(len(cv), n, degree)
points_basis = np.dot(b, cv)
print np.allclose(points_basis, points_scipy)
大部分加速是通过让 coxDeBoor 计算一个向量来实现的
结果而不是一次单个值。请注意，u 已从
传递给 coxDeBoor 的参数。相反，新的 coxDeBoor(k, d) 计算
np.array([coxDeBoor(u[i], k, d) for i in xrange(n)]).
由于 NumPy 数组运算与标量运算具有相同的语法，因此非常 需要更改的代码很少。唯一的句法变化是在最后 条件：
if (d == 0):
return ((u - kv[k] >= 0) & (u - kv[k + 1] < 0)).astype(int)
(u - kv[k] >= 0) 和 (u - kv[k + 1] < 0) 是布尔数组。 astype
将数组值更改为 0 和 1。因此，在单个 0 或 1 之前
返回，现在返回一整组 0 和 1——每个值对应一个
u.
Memoization也可以提高性能，因为递归关系
导致 coxDeBoor(k, d) 被调用以获得 k 和 d 的相同值
不止一次。装饰器语法
@memo
def coxDeBoor(k, d):
...
等同于
def coxDeBoor(k, d):
...
coxDeBoor = memo(coxDeBoor)
和 memo 装饰器导致 coxDeBoor 在 cache 映射中记录
从 (k, d) 对参数到返回值。如果 coxDeBoor(k, d) 是
再次调用，则返回 cache 中的值而不是
重新计算值。
scipy_bspline 仍然更快，但至少 bspline_basis 加上 np.dot 在大概范围内，
如果你想重复使用 b 和许多控制点，cv.
In [109]: %timeit scipy_bspline(cv, 10**5, degree)
100 loops, best of 3: 17.2 ms per loop
In [104]: %timeit b = bspline_basis(len(cv), 10**5, degree)
10 loops, best of 3: 29.1 ms per loop
In [111]: %timeit points_basis = np.dot(b, cv)
100 loops, best of 3: 2.46 ms per loop
fitpack_basis 使用双循环迭代修改 bb 中的元素
和 b。我看不到使用 NumPy 来矢量化这些循环的方法，因为值
bb 和 b 在每个迭代阶段的值取决于来自
以前的迭代。在这种情况下，有时可以使用 Cython
提高循环的性能。
这是 fitpack_basis 的 Cython 化版本，运行 的速度与
参见这个 第页为 关于如何构建代码和 运行 它的说明来自 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 代码对其性能进行基准测试，
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 代码 运行 与
NB.bspline_basis: 0.3322
CB.cython_fitpack_basis: 0.2939
NB.fitpack_basis: 0.9182