使用 numpy/scipy 的快速 b 样条算法 [英] Fast b-spline algorithm with numpy/scipy

问题描述

我需要在 python 中计算 bspline 曲线.我查看了 scipy.interpolate.splprep 和其他一些 scipy 模块，但找不到任何可以轻松满足我需要的东西.所以我在下面写了我自己的模块.代码工作正常，但速度很慢(测试函数在 0.03 秒内运行，考虑到我只要求 100 个具有 6 个控制顶点的样本，这似乎很多).

I need to compute bspline curves in python. I looked into scipy.interpolate.splprep and a few other scipy modules but couldn't find anything that readily gave me what I needed. So i wrote my own module below. The code works fine, but it is slow (test function runs in 0.03s, which seems like a lot considering i'm only asking for 100 samples with 6 control vertices).

有没有办法通过一些 scipy 模块调用来简化下面的代码，这可能会加快速度?如果没有，我可以对我的代码做些什么来提高其性能?

Is there a way to simplify the code below with a few scipy module calls, which presumably would speed it up? And if not, what could i do to my code to improve its performance?

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:

推荐答案

因此，在对我的问题进行了大量研究并进行了大量研究之后，我终于有了答案.一切都在 scipy 中可用，我将我的代码放在这里，希望其他人能发现它有用.

So after obsessing a lot about my question, and much research, i finally have my answer. Everything is available in scipy , and i'm putting my code here so hopefully someone else can find this useful.

该函数接受一组 N-d 个点、一个曲线度数、一个周期状态(打开或关闭)，并将沿该曲线返回 n 个样本.有多种方法可以确保曲线样本等距，但目前我将重点关注这个问题，因为这完全与速度有关.

The function takes in an array of N-d points, a curve degree, a periodic state (opened or closed) and will return n samples along that curve. There are ways to make sure the curve samples are equidistant but for the time being i'll focus on this question, as it is all about speed.

值得注意的是:我似乎无法超越 20 度的曲线.当然，这已经有点矫枉过正了，但我认为值得一提.

Worthy of note: I can't seem to be able to go beyond a curve of 20th degree. Granted, that's overkill already, but i figured it's worth mentioning.

另外值得注意的是:在我的机器上，下面的代码可以在 0.017 秒内计算 100,000 个样本

Also worthy of note: on my machine the code below can calculate 100,000 samples 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

测试:

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

开曲线或周期曲线的结果:

Results for both opened or periodic curves:

从 scipy-0.19.0 开始，有一个新的 scipy.interpolate.BSpline 可以使用的函数.

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

等效性测试:

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

这篇关于使用 numpy/scipy 的快速 b 样条算法的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！