基于numpy/scipy的快速b样条算法

Question

我需要用python计算bspline曲线。我看了看scipy.interpolate.splprep公司和其他一些scipy模块，但找不到任何能给我所需要的东西。所以我在下面写了我自己的模块。
代码运行得很好，但速度很慢（test函数运行时间为0.03s，看起来
就像很多考虑到我只要求100个样品和6个对照品
顶点）。有没有一种方法可以通过几个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)

The two images below shows what my code returns with both closed conditions:

Answer 1

经过研究，我终于有了很多问题
我的回答。在scipy中一切都是可用的，我把我的代码放在这里
希望其他人能发现这个有用。
该函数采用N-d点数组、曲线阶数和周期状态
（打开或关闭）并将沿该曲线返回n个样本。有很多方法
以确保曲线样本是等距的，但目前我将
把注意力集中在这个问题上，因为它完全是关于速度的。
值得注意的是：我似乎无法超越20度的曲线。
当然，这已经是过分了，但我觉得值得一提。
同样值得注意的是：在我的机器上下面的代码可以计算100000
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

To test it:

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

As of scipy-0.19.0 there is a new
scipy.interpolate.BSpline
function that can be used.

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

Testing for equivalency:

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