# [python] Wie bekomme ich die Spline-Basis von scipy.interpolate.splev?

Hier ist eine Version von `coxDeBoor` die (auf meinem Rechner) 840x schneller als das Original ist.

``````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)``````

Der Großteil der Beschleunigung wird dadurch erreicht, dass coxDeBoor einen Vektor von Ergebnissen anstelle von jeweils einem einzelnen Wert berechnet. Beachten Sie, dass Sie aus den an `coxDeBoor` Argumenten `coxDeBoor` . Stattdessen berechnet das neue `coxDeBoor(k, d)` , was vor `np.array([coxDeBoor(u[i], k, d) for i in xrange(n)])` .

Da NumPy-Arithmetik die gleiche Syntax wie Skalararithmetik hat, ist nur sehr wenig Code zum Ändern erforderlich. Die einzige syntaktische Änderung war in der Endbedingung:

``````if (d == 0):
return ((u - kv[k] >= 0) & (u - kv[k + 1] < 0)).astype(int)``````

`(u - kv[k] >= 0)` und `(u - kv[k + 1] < 0)` sind boolesche Arrays. `astype` ändert die Array-Werte auf 0 und 1. Während also zuvor eine einzige 0 oder 1 zurückgegeben wurde, wird jetzt ein ganzes Array von 0s und 1s zurückgegeben - eins für jeden Wert in `u` .

Memoisierung kann auch die Leistung verbessern, da die Wiederholungsrelationen bewirken, dass `coxDeBoor(k, d)` für dieselben Werte von `k` und `d` mehr als einmal aufgerufen wird. Die Dekoratorsyntax

``````@memo
def coxDeBoor(k, d):
...``````

ist äquivalent zu

``````def coxDeBoor(k, d):
...
coxDeBoor = memo(coxDeBoor)``````

Der `memo` Decorator bewirkt, dass `coxDeBoor` in einem `cache` eine Abbildung von `(k, d)` Paaren von Argumenten auf zurückgegebene Werte aufzeichnet. Wenn `coxDeBoor(k, d)` erneut aufgerufen wird, wird der Wert aus dem `cache` zurückgegeben, anstatt den Wert neu zu berechnen.

`scipy_bspline` ist zwar immer noch schneller, aber mindestens `bspline_basis` plus `np.dot` ist im `np.dot` und kann nützlich sein, wenn du `b` mit vielen Kontrollpunkten, `cv` , wiederverwenden willst.

``````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``````
Question

Ich muss B-Splines in Python bewerten. Dazu habe ich den folgenden Code geschrieben, der sehr gut funktioniert.

``````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``````

Wenn Sie ihm 6 Kontrollpunkte geben und ihn bitten, 100k Punkte auf der Kurve zu bewerten, ist das ein Kinderspiel:

``````# 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``````

Hier ist eine Handlung von "points_scipy":

Jetzt, um das schneller zu machen, könnte ich eine Basis für alle 100k Punkte auf der Kurve berechnen, diese im Speicher speichern, und wenn ich eine Kurve zeichnen muss, muss ich nur die neuen Kontrollpunktpositionen mit der gespeicherten Basis multiplizieren hol dir die neue Kurve. Um meinen Standpunkt zu beweisen, habe ich eine Funktion geschrieben, die den Algorithmus von DeBoor benutzt, um meine Basis zu berechnen:

``````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``````

Nun können wir damit eine neue Basis berechnen, diese mit den Kontrollpunkten multiplizieren und bestätigen, dass wir die gleichen Ergebnisse wie splev erhalten:

``````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``````

Meine extrem langsame Funktion lieferte 100k Basiswerte in 11 Sekunden, aber da diese Werte nur einmal berechnet werden müssen, war die Berechnung der Punkte auf der Kurve auf diese Weise 6 mal schneller als über splev.

Die Tatsache, dass ich in der Lage war, die gleichen Ergebnisse sowohl von meiner Methode und Splev zu bekommen, führt mich zu der Annahme, dass intern Splev wahrscheinlich eine Basis wie ich, außer viel schneller berechnet.

Mein Ziel ist es, herauszufinden, wie ich meine Basis schnell berechnen, im Speicher speichern und np.dot () verwenden kann, um die neuen Punkte auf der Kurve zu berechnen, und meine Frage lautet: Ist es möglich, pixy.interpolate zu verwenden, um die Basiswerte, die (wie ich vermute) splev verwendet, um das Ergebnis zu berechnen? Und wenn ja, wie?

# [NACHTRAG]

Nach unutbus und ev-brs sehr nützlicher Einsicht, wie scipy die Basis eines Splines berechnet, habe ich den Fortran-Code nachgeschlagen und ein Äquivalent zu den besten meiner Fähigkeiten geschrieben:

``````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``````

Testen gegen die Unutbu-Version meiner ursprünglichen Basisfunktion:

``````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``````

Meine Funktion ist 5 mal langsamer, aber immer noch relativ schnell. Was ich daran mag, ist, dass es mir erlaubt, Musterbereiche zu verwenden, die über die Grenzen hinausgehen, was in meiner Anwendung nützlich ist. Beispielsweise:

``````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 ]]``````

Aus diesem Grund werde ich wahrscheinlich fitpack_basis verwenden, da es relativ schnell ist. Aber ich würde gerne Vorschläge zur Verbesserung seiner Leistung und hoffe, näher an die unutbu-Version der ursprünglichen Basisfunktion, die ich schrieb, zu kommen.