numpy/scipyを使った高速bスプラインアルゴリズム

Question

私はPythonで3D bspline曲線を計算する必要があります。私はscipy.interpolate.splprepといくつかの他のscipyモジュールを調べましたが、私が必要とするものを簡単に私に与えたものは見つかりませんでした。だから私は自分の下に自分のモジュールを書きました。コードは正常に動作しますが、遅いです（テスト関数は0.03秒で実行されますが、これは6つのコントロール頂点を持つ100個のサンプルしか求めていないと思われます）。

いくつかのscipyモジュール呼び出しを使って以下のコードを単純化する方法はありますか？そうでない場合は、パフォーマンスを向上させるためにコードにどうすればよいですか？

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) 

以下の二つの画像は、私のコードは、閉じた状態の両方で返すものを示しています データをプロファイリングすることなく、最適化のヒントを与える

Answer 1

私の質問と多くの研究について多くのことに執着した後、ついに私は答えを得ました。すべてがscipyで利用可能です、私はここに自分のコードを入れていますので、他の誰かがこれを見つけることができればうれしいです。

この関数は、N-d点、曲線の次数、周期的な状態（開または閉）の配列をとり、その曲線に沿ってn個のサンプルを返します。曲線のサンプルが等距離であることを確認する方法はありますが、当面はこの問題に焦点を当てます。スピードについてはすべてです。

値する価値：私は20度のカーブを越えることができないようです。確かに、それはすでに過労ですが、私はそれが言及する価値があると考えました。

ノートのも値する

：両方のオープンまたは定期的なカーブの

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

結果：：

私のマシン上で以下のコードは、それをテストするには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 

10万個のサンプルを計算することができます

追加

scipy-0.19.0からは、新しい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)) 

テスト：

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

はビット暗闇での撮影のようなものです。しかし、関数coxDeBoorは非常に頻繁に呼び出されるようです。ここで私は最適化を開始します。

Pythonの関数呼び出しare expensive。過剰な関数呼び出しを避けるために、反復をcoxDeBoorの繰り返しに置き換えるようにしてください。これを行う方法の一般的な情報は、this questionへの回答に記載されています。スタック/キューとしてはcollections.dequeを使用できます。