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



Answers

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.




Links