使用numpy/scipy的快速B样条算法

Question

我需要在Python中计算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)

下面的两张图片显示了我的代码在两种关闭条件下返回的内容：

- Fnord

---

我写了一个更好的实现你的工作 https://gist.github.com/soulslicer/1224bfc6a81f25835054cadf18325251 - raaj

Answer 1

因为一直困扰着我这个问题，所以我经过了大量研究，终于得到了答案。scipy 库中提供了一切需要的内容，我将在此放置我的代码，希望能对其他人有所帮助。

该函数接收一个 N 维点的数组、曲线度数、周期状态（开放式或封闭式），并返回沿曲线的 n 个样本。有方法可以确保曲线采样均匀分布，但目前我将专注于这个问题，因为它涉及速度问题。

值得注意的是：我似乎无法超过 20 次方的曲线。当然，那已经相当夸张了，但我认为这是值得一提的。

另值得注意的是：在我的电脑上，下面的代码可以在 0.017 秒内计算出 100,000 个样本。

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

测试方法：

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

针对开放或周期性曲线的结果：

附录

从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

- Fnord

---

太棒了，非常感谢！在我的应用程序中，我有两个端点的三维坐标和一堆已知的三维控制点。它能够非常出色地绘制样条曲线！太棒了！我正在处理3D图像数据的ndarrays。 - kabammi

太好了！你促使我编辑了我的答案并删除了最后的for循环，这是不必要的。我还在最后添加了一个附录，提到了在scipy 0.19.0中添加的官方BSpline函数。 - Fnord

嗯...我在使用你的scipy_bspline函数时遇到了错误。我将一个列表作为CV传递，所以在你原来的函数中使用cv = np.asarray(cv)很有帮助。然后我使用degree=5，新函数抛出一个错误并告诉我我需要至少12个节点...旧代码不关心这个问题，而且可以正常工作。所以对我来说，旧代码胜出。 :) - kabammi

@kabammi 哦，是的，我忘记了将曲线的度数夹紧以保持其稳定。已修复 :) - Fnord

@Abhijeet，我的用例中控制顶点“cv”是三维的。这就是我所说的3D。但是代码可以处理任何维度（例如2D，就像我的示例所示）。 - Fnord

Answer 2

没有性能数据的优化建议就像是盲目射击一样。然而，函数coxDeBoor似乎被频繁调用。这是我开始优化的地方。

在Python中，函数调用很耗费资源。您应该尝试使用迭代替换coxDeBoor递归以避免过多的函数调用。有关如何进行此操作的一些常规信息可以在这个问题的答案中找到。您可以使用collections.deque作为堆栈/队列。

- MB-F