【问题标题】:B-spline interpolation with Python使用 Python 进行 B 样条插值
【发布时间】:2014-08-28 01:20:00
【问题描述】:

我正在尝试使用 Python 重现 B 样条的 Mathematica 示例。

mathematica 示例代码如下

pts = {{0, 0}, {0, 2}, {2, 3}, {4, 0}, {6, 3}, {8, 2}, {8, 0}};
Graphics[{BSplineCurve[pts, SplineKnots -> {0, 0, 0, 0, 2, 3, 4, 6, 6, 6, 6}], Green, Line[pts], Red, Point[pts]}]

并产生我所期望的。现在我尝试用 Python/scipy 做同样的事情:

import numpy as np
import matplotlib.pyplot as plt
import scipy.interpolate as si

points = np.array([[0, 0], [0, 2], [2, 3], [4, 0], [6, 3], [8, 2], [8, 0]])
x = points[:,0]
y = points[:,1]

t = range(len(x))
knots = [2, 3, 4]
ipl_t = np.linspace(0.0, len(points) - 1, 100)

x_tup = si.splrep(t, x, k=3, t=knots)
y_tup = si.splrep(t, y, k=3, t=knots)
x_i = si.splev(ipl_t, x_tup)
y_i = si.splev(ipl_t, y_tup)

print 'knots:', x_tup

fig = plt.figure()
ax = fig.add_subplot(111)
plt.plot(x, y, label='original')
plt.plot(x_i, y_i, label='spline')
plt.xlim([min(x) - 1.0, max(x) + 1.0])
plt.ylim([min(y) - 1.0, max(y) + 1.0])
plt.legend()
plt.show()

这会导致一些也被插值但看起来不太正确的东西。我使用与mathematica 相同的结分别对x 和y 分量进行参数化和样条化。但是,我得到了过冲和下冲,这使我的插值曲线弯曲到控制点的凸包之外。这样做的正确方法是什么/mathematica 是如何做到的?

【问题讨论】:

    标签: python bspline


    【解决方案1】:

    我能够使用 Python/scipy 重新创建我在上一篇文章中询问的 Mathematica 示例。结果如下:

    B 样条,非周期

    诀窍是截取系数,即scipy.interpolate.splrep 返回的元组的元素 1,并在将它们交给scipy.interpolate.splev 之前将它们替换为控制点值,或者,如果您可以创建自己打结,你也可以不用splrep,自己创建整个元组。

    然而,这一切的奇怪之处在于,根据手册,splrep 返回(并且splev 期望)一个元组,其中包含一个样条系数向量,每个结一个系数。但是,根据我发现的所有资料,样条被定义为 N_control_points 个基本样条的加权和,因此我希望系数向量具有与控制点一样多的元素,而不是节点位置。

    事实上,当将splrep 的结果元组与如上所述修改的系数向量提供给scipy.interpolate.splev 时,结果表明该向量的第一个 N_control_points 实际上是预期的系数N_control_points 基本样条。该向量的最后一个 degree + 1 元素似乎没有效果。我很困惑为什么会这样。如果有人能澄清这一点,那就太好了。以下是生成上述图的来源:

    import numpy as np
    import matplotlib.pyplot as plt
    import scipy.interpolate as si
    
    points = [[0, 0], [0, 2], [2, 3], [4, 0], [6, 3], [8, 2], [8, 0]];
    points = np.array(points)
    x = points[:,0]
    y = points[:,1]
    
    t = range(len(points))
    ipl_t = np.linspace(0.0, len(points) - 1, 100)
    
    x_tup = si.splrep(t, x, k=3)
    y_tup = si.splrep(t, y, k=3)
    
    x_list = list(x_tup)
    xl = x.tolist()
    x_list[1] = xl + [0.0, 0.0, 0.0, 0.0]
    
    y_list = list(y_tup)
    yl = y.tolist()
    y_list[1] = yl + [0.0, 0.0, 0.0, 0.0]
    
    x_i = si.splev(ipl_t, x_list)
    y_i = si.splev(ipl_t, y_list)
    
    #==============================================================================
    # Plot
    #==============================================================================
    
    fig = plt.figure()
    
    ax = fig.add_subplot(231)
    plt.plot(t, x, '-og')
    plt.plot(ipl_t, x_i, 'r')
    plt.xlim([0.0, max(t)])
    plt.title('Splined x(t)')
    
    ax = fig.add_subplot(232)
    plt.plot(t, y, '-og')
    plt.plot(ipl_t, y_i, 'r')
    plt.xlim([0.0, max(t)])
    plt.title('Splined y(t)')
    
    ax = fig.add_subplot(233)
    plt.plot(x, y, '-og')
    plt.plot(x_i, y_i, 'r')
    plt.xlim([min(x) - 0.3, max(x) + 0.3])
    plt.ylim([min(y) - 0.3, max(y) + 0.3])
    plt.title('Splined f(x(t), y(t))')
    
    ax = fig.add_subplot(234)
    for i in range(7):
        vec = np.zeros(11)
        vec[i] = 1.0
        x_list = list(x_tup)
        x_list[1] = vec.tolist()
        x_i = si.splev(ipl_t, x_list)
        plt.plot(ipl_t, x_i)
    plt.xlim([0.0, max(t)])
    plt.title('Basis splines')
    plt.show()
    

    B 样条,周期性

    现在为了创建如下所示的闭合曲线,这是另一个可以在网络上找到的 Mathematica 示例,

    有必要在splrep调用中设置per参数,如果你使用它。在最后用 degree+1 值填充控制点列表后,这似乎工作得很好,如图所示。

    然而,这里的下一个特点是系数向量中的第一个和最后一个 degree 元素没有影响,这意味着控制点必须从第二个位置开始放入向量中,即位置1。只有这样结果才ok。对于度数 k=4 和 k=5,该位置甚至会更改为位置 2。

    这里是生成闭合曲线的来源:

    import numpy as np
    import matplotlib.pyplot as plt
    import scipy.interpolate as si
    
    points = [[-2, 2], [0, 1], [-2, 0], [0, -1], [-2, -2], [-4, -4], [2, -4], [4, 0], [2, 4], [-4, 4]]
    
    degree = 3
    
    points = points + points[0:degree + 1]
    points = np.array(points)
    n_points = len(points)
    x = points[:,0]
    y = points[:,1]
    
    t = range(len(x))
    ipl_t = np.linspace(1.0, len(points) - degree, 1000)
    
    x_tup = si.splrep(t, x, k=degree, per=1)
    y_tup = si.splrep(t, y, k=degree, per=1)
    x_list = list(x_tup)
    xl = x.tolist()
    x_list[1] = [0.0] + xl + [0.0, 0.0, 0.0, 0.0]
    
    y_list = list(y_tup)
    yl = y.tolist()
    y_list[1] = [0.0] + yl + [0.0, 0.0, 0.0, 0.0]
    
    x_i = si.splev(ipl_t, x_list)
    y_i = si.splev(ipl_t, y_list)
    
    #==============================================================================
    # Plot
    #==============================================================================
    
    fig = plt.figure()
    
    ax = fig.add_subplot(231)
    plt.plot(t, x, '-og')
    plt.plot(ipl_t, x_i, 'r')
    plt.xlim([0.0, max(t)])
    plt.title('Splined x(t)')
    
    ax = fig.add_subplot(232)
    plt.plot(t, y, '-og')
    plt.plot(ipl_t, y_i, 'r')
    plt.xlim([0.0, max(t)])
    plt.title('Splined y(t)')
    
    ax = fig.add_subplot(233)
    plt.plot(x, y, '-og')
    plt.plot(x_i, y_i, 'r')
    plt.xlim([min(x) - 0.3, max(x) + 0.3])
    plt.ylim([min(y) - 0.3, max(y) + 0.3])
    plt.title('Splined f(x(t), y(t))')
    
    ax = fig.add_subplot(234)
    for i in range(n_points - degree - 1):
        vec = np.zeros(11)
        vec[i] = 1.0
        x_list = list(x_tup)
        x_list[1] = vec.tolist()
        x_i = si.splev(ipl_t, x_list)
        plt.plot(ipl_t, x_i)
    plt.xlim([0.0, 9.0])
    plt.title('Periodic basis splines')
    
    plt.show()
    

    B 样条,周期性,更高次

    最后,还有一个我也无法解释的效果,就是在去5度的时候,样条曲线出现了一个小的不连续性,见右上图,这是一个特写那个“半月形鼻形”。下面列出了产​​生这个的源代码。

    import numpy as np
    import matplotlib.pyplot as plt
    import scipy.interpolate as si
    
    points = [[-2, 2], [0, 1], [-2, 0], [0, -1], [-2, -2], [-4, -4], [2, -4], [4, 0], [2, 4], [-4, 4]]
    
    degree = 5
    
    points = points + points[0:degree + 1]
    points = np.array(points)
    n_points = len(points)
    x = points[:,0]
    y = points[:,1]
    
    t = range(len(x))
    ipl_t = np.linspace(1.0, len(points) - degree, 1000)
    
    knots = np.linspace(-degree, len(points), len(points) + degree + 1).tolist()
    
    xl = x.tolist()
    coeffs_x = [0.0, 0.0] + xl + [0.0, 0.0, 0.0]
    
    yl = y.tolist()
    coeffs_y = [0.0, 0.0] + yl + [0.0, 0.0, 0.0]
    
    x_i = si.splev(ipl_t, (knots, coeffs_x, degree))
    y_i = si.splev(ipl_t, (knots, coeffs_y, degree))
    
    #==============================================================================
    # Plot
    #==============================================================================
    
    fig = plt.figure()
    
    ax = fig.add_subplot(231)
    plt.plot(t, x, '-og')
    plt.plot(ipl_t, x_i, 'r')
    plt.xlim([0.0, max(t)])
    plt.title('Splined x(t)')
    
    ax = fig.add_subplot(232)
    plt.plot(t, y, '-og')
    plt.plot(ipl_t, y_i, 'r')
    plt.xlim([0.0, max(t)])
    plt.title('Splined y(t)')
    
    ax = fig.add_subplot(233)
    plt.plot(x, y, '-og')
    plt.plot(x_i, y_i, 'r')
    plt.xlim([min(x) - 0.3, max(x) + 0.3])
    plt.ylim([min(y) - 0.3, max(y) + 0.3])
    plt.title('Splined f(x(t), y(t))')
    
    ax = fig.add_subplot(234)
    for i in range(n_points - degree - 1):
        vec = np.zeros(11)
        vec[i] = 1.0
        x_i = si.splev(ipl_t, (knots, vec, degree))
        plt.plot(ipl_t, x_i)
    plt.xlim([0.0, 9.0])
    plt.title('Periodic basis splines')
    
    plt.show()
    

    鉴于 b 样条曲线在科学界无处不在,并且 scipy 是一个如此全面的工具箱,而且我无法在网上找到很多关于我在此询问的内容,这让我相信我我在错误的轨道上或俯瞰某些东西。任何帮助将不胜感激。

    【讨论】:

    • 这里给出了它为什么会这样工作的解释:brnt.eu/phd/node11.html(您可能想从“样条线作为 B 样条线的线性组合”开始阅读)。为了稍微捍卫 scipy,这不是 scupy 自己的算法,而仅仅是著名的 fortran 库 FITPACK 的包装。
    【解决方案2】:

    使用我为another question i asked here.写的这个函数

    在我的问题中,我正在寻找使用 scipy 计算 bsplines 的方法(这就是我实际上偶然发现您的问题的方式)。

    经过深思熟虑,我想出了下面的功能。它会评估任何高达 20 度的曲线(比我们需要的多)。在速度方面,我测试了 100,000 个样本,耗时 0.017 秒

    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,dtype='int')
        else:
            kv = np.concatenate(([0]*degree, np.arange(count-degree+1), [count-degree]*degree))
    
    
        # Calculate query range
        u = np.linspace(periodic,(count-degree),n)
    
    
        # Calculate result
        return np.array(si.splev(u, (kv,cv.T,degree))).T
    

    开放曲线和周期曲线的结果:

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

    【讨论】:

    • 这适用于 3 个维度(或任何维度)吗?
    • @kureta 任何维度
    • 谢谢。这使我免于数小时的麻烦。
    • @Fnord 适用于封闭的 B 样条曲线,非周期性的问题出现在 kv = np.array([0]*degree + range(count-degree+1) + [count-degree]*degree,dtype='int')
    • @nzou 在 Python 3 中特别是在该行中遇到错误。知道为什么吗? “TypeError: can only concatenate list (not "range") to list”尖叫着那些方括号试图成为数组。
    【解决方案3】:

    我相信 scipy 的 fitpack 库正在做的事情比 Mathematica 所做的更复杂。我也对发生的事情感到困惑。

    这些函数中都有smoothing参数,默认的插值行为是尽量让点穿过线。这就是 fitpack 软件所做的,所以我猜 scipy 只是继承了它? (http://www.netlib.org/fitpack/all -- 我不确定这是不是合适的背包)

    我从http://research.microsoft.com/en-us/um/people/ablake/contours/ 中获取了一些想法,并使用其中的 B 样条对您的示例进行了编码。

    import numpy
    
    import matplotlib.pyplot as plt
    
    # This is the basis function described in eq 3.6 in http://research.microsoft.com/en-us/um/people/ablake/contours/
    def func(x, offset):
        out = numpy.ndarray((len(x)))
    
        for i, v in enumerate(x):
            s = v - offset
    
            if s >= 0 and s < 1:
                out[i] = s * s / 2.0
            elif s >= 1 and s < 2:
                out[i] = 3.0 / 4.0 - (s - 3.0 / 2.0) * (s - 3.0 / 2.0)
            elif s >= 2 and s < 3:
                out[i] = (s - 3.0) * (s - 3.0) / 2.0
            else:
                out[i] = 0.0
    
        return out
    
    # We have 7 things to fit, so let's do 7 basis functions?
    y = numpy.array([0, 2, 3, 0, 3, 2, 0])
    
    # We need enough x points for all the basis functions... That's why the weird linspace max here
    x = numpy.linspace(0, len(y) + 2, 100)
    
    B = numpy.ndarray((len(x), len(y)))
    
    for k in range(len(y)):
        B[:, k] = func(x, k)
    
    plt.plot(x, B.dot(y))
    # The x values in the next statement are the maximums of each basis function. I'm not sure at all this is right
    plt.plot(numpy.array(range(len(y))) + 1.5, y, '-o')
    plt.legend('B-spline', 'Control points')
    plt.show()
    
    for k in range(len(y)):
        plt.plot(x, B[:, k])
    plt.title('Basis functions')
    plt.show()
    

    无论如何,我认为其他人也有同样的问题,看看: Behavior of scipy's splrep

    【讨论】:

      猜你喜欢
      • 1970-01-01
      • 1970-01-01
      • 1970-01-01
      • 2020-08-25
      • 2014-11-07
      • 2018-02-12
      • 2016-08-07
      • 1970-01-01
      • 2015-08-11
      相关资源
      最近更新 更多