使用 de Boor 算法的 B 样条导数

pac*_*der 5 python algorithm interpolation bspline

维基百科为我们提供了 de Boor 算法的 Python 实现:

def deBoor(k, x, t, c, p):
    """
    Evaluates S(x).

    Args
    ----
    k: index of knot interval that contains x
    x: position
    t: array of knot positions, needs to be padded as described above
    c: array of control points
    p: degree of B-spline
    """
    d = [c[j + k - p] for j in range(0, p+1)]

    for r in range(1, p+1):
        for j in range(p, r-1, -1):
            alpha = (x - t[j+k-p]) / (t[j+1+k-r] - t[j+k-p])
            d[j] = (1.0 - alpha) * d[j-1] + alpha * d[j]

    return d[p]
Run Code Online (Sandbox Code Playgroud)

是否有类似的算法计算 B-Spline 插值曲线的导数(甚至 n 阶导数)?

我知道在数学上它被简化为使用低阶样条,但不能将其应用于 de Boor 算法。

pac*_*der 5

我想我找到了重新使用 de Boor 曲线导数算法的正确方法。

\n\n

首先,我们考虑 B 样条曲线的定义。它是控制点的线性组合:\n曲线      (1)

\n\n

因此,导数是基函数导数的线性组合

\n\n

曲线      (2)

\n\n

基函数的导数定义如下:

\n\n

曲线      (3)

\n\n

我们将(3)插入(2),经过一些代数功夫(此处描述)http://public.vrac.iastate.edu/~oliver/courses/me625/week5b.pdf,我们得到:

\n\n

曲线      (4),\n其中曲线

\n\n

B 样条曲线的导数只不过是建立在新控制点 Q 之上的 (p-1) 次新 B 样条曲线。\n现在,为了采用 de Boor 算法,我们计算设置新的控制点并将样条曲线度数 p 降低 1:

\n\n
def deBoorDerivative(k, x, t, c, p):\n    """\n    Evaluates S(x).\n\n    Args\n    ----\n    k: index of knot interval that contains x\n    x: position\n    t: array of knot positions, needs to be padded as described above\n    c: array of control points\n    p: degree of B-spline\n    """\n    q = [p * (c[j+k-p+1] - c[j+k-p]) / (t[j+k+1] - t[j+k-p+1]) for j in range(0, p)]\n\n    for r in range(1, p):\n        for j in range(p-1, r-1, -1):\n            right = j+1+k-r\n            left = j+k-(p-1)\n            alpha = (x - t[left]) / (t[right] - t[left])\n            q[j] = (1.0 - alpha) * q[j-1] + alpha * q[j]\n\n    return q[p-1]\n
Run Code Online (Sandbox Code Playgroud)\n\n

测试:

\n\n
import numpy as np\nimport math as m\n\npoints = np.array([[i, m.sin(i / 3.0), m.cos(i / 2)] for i in range(0, 11)])\nknots = np.array([0, 0, 0, 0, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 1.0, 1.0, 1.0, 1.0])\n\n\ndef finiteDifferenceDerivative(k, x, t, c, p):\n    """ Third order finite difference derivative """\n\n    f = lambda xx : deBoor(k, xx, t, c, p)\n\n    dx = 1e-7\n\n    return (- f(x + 2 * dx) \\\n            + 8 * f(x + dx) \\\n            - 8 * f(x - dx) \\\n            + f(x - 2 * dx)) / ( 12 * dx )\n\n\nprint "Derivatives: "\xc2\xb7\nprint "De Boor:\\t", deBoorDerivative(7, 0.44, knots, points, 3)\nprint "Finite Difference:\\t", finiteDifferenceDerivative(7, 0.44, knots, points, 3)\n\n
Run Code Online (Sandbox Code Playgroud)\n\n

输出:

\n\n
Derivatives: \nDe Boor:              [10. 0.36134438  2.63969004]\nFinite Difference:    [9.99999999 0.36134438 2.63969004]\n
Run Code Online (Sandbox Code Playgroud)\n