Algoritmo b-spline veloce con numpy/scipy

Question

Ho bisogno di calcolare curve bspline 3D in python. Ho esaminato scipy.interpolate.splprep e alcuni altri moduli scipy ma non ho trovato nulla che mi fornisse prontamente ciò di cui avevo bisogno. Quindi ho scritto il mio modulo qui sotto. Il codice funziona bene, ma è lento (la funzione di test gira in 0.03s, il che sembra molto considerando che sto solo chiedendo 100 campioni con 6 vertici di controllo).

C'è un modo per semplificare il codice sottostante con alcune chiamate scipy di modulo, che presumibilmente lo renderebbero più veloce? E se no, cosa potrei fare al mio codice per migliorare le sue prestazioni?

import numpy as np 

# cv = np.array of 3d control vertices 
# n = number of samples (default: 100) 
# d = curve degree (default: cubic) 
# closed = is the curve closed (periodic) or open? (default: open) 
def bspline(cv, n=100, d=3, closed=False): 

    # Create a range of u values 
    count = len(cv) 
    knots = None 
    u = None 
    if not closed: 
     u = np.arange(0,n,dtype='float')/(n-1) * (count-d) 
     knots = np.array([0]*d + range(count-d+1) + [count-d]*d,dtype='int') 
    else: 
     u = ((np.arange(0,n,dtype='float')/(n-1) * count) - (0.5 * (d-1))) % count # keep u=0 relative to 1st cv 
     knots = np.arange(0-d,count+d+d-1,dtype='int') 


    # Simple Cox - DeBoor recursion 
    def coxDeBoor(u, k, d): 

     # Test for end conditions 
     if (d == 0): 
      if (knots[k] <= u and u < knots[k+1]): 
       return 1 
      return 0 

     Den1 = knots[k+d] - knots[k] 
     Den2 = knots[k+d+1] - knots[k+1] 
     Eq1 = 0; 
     Eq2 = 0; 

     if Den1 > 0: 
      Eq1 = ((u-knots[k])/Den1) * coxDeBoor(u,k,(d-1)) 
     if Den2 > 0: 
      Eq2 = ((knots[k+d+1]-u)/Den2) * coxDeBoor(u,(k+1),(d-1)) 

     return Eq1 + Eq2 


    # Sample the curve at each u value 
    samples = np.zeros((n,3)) 
    for i in xrange(n): 
     if not closed: 
      if u[i] == count-d: 
       samples[i] = np.array(cv[-1]) 
      else: 
       for k in xrange(count): 
        samples[i] += coxDeBoor(u[i],k,d) * cv[k] 

     else: 
      for k in xrange(count+d): 
       samples[i] += coxDeBoor(u[i],k,d) * cv[k%count] 


    return samples 




if __name__ == "__main__": 
    import matplotlib.pyplot as plt 
    def test(closed): 
     cv = np.array([[ 50., 25., -0.], 
       [ 59., 12., -0.], 
       [ 50., 10., 0.], 
       [ 57., 2., 0.], 
       [ 40., 4., 0.], 
       [ 40., 14., -0.]]) 

     p = bspline(cv,closed=closed) 
     x,y,z = p.T 
     cv = cv.T 
     plt.plot(cv[0],cv[1], 'o-', label='Control Points') 
     plt.plot(x,y,'k-',label='Curve') 
     plt.minorticks_on() 
     plt.legend() 
     plt.xlabel('x') 
     plt.ylabel('y') 
     plt.xlim(35, 70) 
     plt.ylim(0, 30) 
     plt.gca().set_aspect('equal', adjustable='box') 
     plt.show() 

    test(False) 

Le due immagini qui sotto mostra ciò che il mio codice restituisce con entrambe le condizioni chiusi:

Answer 1

Quindi dopo aver ossessionato molto sulla mia domanda e molte ricerche, ho finalmente la mia risposta. Tutto è disponibile in Scipy, e sto mettendo il mio codice qui, quindi spero che qualcun altro possa trovarlo utile.

La funzione comprende una serie di punti N-d, un grado di curva, uno stato periodico (aperto o chiuso) e restituirà n campioni lungo tale curva. Ci sono modi per assicurarci che i campioni delle curve siano equidistanti, ma per il momento mi concentrerò su questa domanda, poiché si tratta di velocità.

Degno di nota: non riesco a superare la curva del 20 ° grado. Certo, è già eccessivo, ma ho pensato che sia degno di nota.

Degno di nota: sulla mia macchina il codice qui sotto può calcolare 100.000 campioni in 0.017s

import numpy as np 
import scipy.interpolate as si 


def bspline(cv, n=100, degree=3, periodic=False): 
    """ Calculate n samples on a bspline 

     cv :  Array ov control vertices 
     n :  Number of samples to return 
     degree: Curve degree 
     periodic: True - Curve is closed 
        False - Curve is open 
    """ 

    # If periodic, extend the point array by count+degree+1 
    cv = np.asarray(cv) 
    count = len(cv) 

    if periodic: 
     factor, fraction = divmod(count+degree+1, count) 
     cv = np.concatenate((cv,) * factor + (cv[:fraction],)) 
     count = len(cv) 
     degree = np.clip(degree,1,degree) 

    # If opened, prevent degree from exceeding count-1 
    else: 
     degree = np.clip(degree,1,count-1) 


    # Calculate knot vector 
    kv = None 
    if periodic: 
     kv = np.arange(0-degree,count+degree+degree-1) 
    else: 
     kv = np.clip(np.arange(count+degree+1)-degree,0,count-degree) 

    # Calculate query range 
    u = np.linspace(periodic,(count-degree),n) 


    # Calculate result 
    return np.array(si.splev(u, (kv,cv.T,degree))).T 

Per provarlo:

import matplotlib.pyplot as plt 
colors = ('b', 'g', 'r', 'c', 'm', 'y', 'k') 

cv = np.array([[ 50., 25.], 
    [ 59., 12.], 
    [ 50., 10.], 
    [ 57., 2.], 
    [ 40., 4.], 
    [ 40., 14.]]) 

plt.plot(cv[:,0],cv[:,1], 'o-', label='Control Points') 

for d in range(1,21): 
    p = bspline(cv,n=100,degree=d,periodic=True) 
    x,y = p.T 
    plt.plot(x,y,'k-',label='Degree %s'%d,color=colors[d%len(colors)]) 

plt.minorticks_on() 
plt.legend() 
plt.xlabel('x') 
plt.ylabel('y') 
plt.xlim(35, 70) 
plt.ylim(0, 30) 
plt.gca().set_aspect('equal', adjustable='box') 
plt.show() 

Risultati per entrambe le curve aperte o periodiche:

ADDENDUM

A partire da scipy-0.19.0 è disponibile una nuova funzione scipy.interpolate.BSpline.

import numpy as np 
import scipy.interpolate as si 
def scipy_bspline(cv, n=100, degree=3, periodic=False): 
    """ Calculate n samples on a bspline 

     cv :  Array ov control vertices 
     n :  Number of samples to return 
     degree: Curve degree 
     periodic: True - Curve is closed 
    """ 
    cv = np.asarray(cv) 
    count = cv.shape[0] 

    # Closed curve 
    if periodic: 
     kv = np.arange(-degree,count+degree+1) 
     factor, fraction = divmod(count+degree+1, count) 
     cv = np.roll(np.concatenate((cv,) * factor + (cv[:fraction],)),-1,axis=0) 
     degree = np.clip(degree,1,degree) 

    # Opened curve 
    else: 
     degree = np.clip(degree,1,count-1) 
     kv = np.clip(np.arange(count+degree+1)-degree,0,count-degree) 

    # Return samples 
    max_param = count - (degree * (1-periodic)) 
    spl = si.BSpline(kv, cv, degree) 
    return spl(np.linspace(0,max_param,n)) 

Test per equivalenza:

p1 = bspline(cv,n=10**6,degree=3,periodic=True) # 1 million samples: 0.0882 sec 
p2 = scipy_bspline(cv,n=10**6,degree=3,periodic=True) # 1 million samples: 0.0789 sec 
print np.allclose(p1,p2) # returns True 

Answer 2

Dare suggerimenti per l'ottimizzazione senza profilazione dei dati è un po 'come sparare al buio. Tuttavia, la funzione coxDeBoor sembra essere chiamata molto spesso. Questo è dove vorrei iniziare l'ottimizzazione.

Chiamate di funzione in Python are expensive. Si dovrebbe provare a sostituire la ricorsione coxDeBoor con iterazione per evitare chiamate di funzione eccessive. Alcune informazioni generali su come eseguire questa operazione sono disponibili nelle risposte a this question. Come stack/coda puoi usare collections.deque.