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

Question

维基百科为我们提供了 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]

有没有类似的算法计算B样条插值曲线的导数(甚至n阶导数)？

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

Answer 1

我认为我找到了将 de Boor 算法重新用于曲线导数的正确方法。

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

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

(2)

基函数的导数定义如下:

(3)

我们将 (3) 代入 (2) 并在进行了一些代数功夫之后，此处描述 http://public.vrac.iastate.edu/~oliver/courses/me625/week5b.pdf ，我们得到:

(4), 其中

B 样条曲线的导数不过是在新控制点 Q 之上构建的 (p-1) 次的新 B 样条曲线。 现在，为了使用 de Boor 算法，我们计算新的控制点集并将样条次数 p 降低 1:

def deBoorDerivative(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
    """
    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)]

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

    return q[p-1]

测试:

import numpy as np
import math as m

points = np.array([[i, m.sin(i / 3.0), m.cos(i / 2)] for i in range(0, 11)])
knots = 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])


def finiteDifferenceDerivative(k, x, t, c, p):
    """ Third order finite difference derivative """

    f = lambda xx : deBoor(k, xx, t, c, p)

    dx = 1e-7

    return (- f(x + 2 * dx) \
            + 8 * f(x + dx) \
            - 8 * f(x - dx) \
            + f(x - 2 * dx)) / ( 12 * dx )


print "Derivatives: "·
print "De Boor:\t", deBoorDerivative(7, 0.44, knots, points, 3)
print "Finite Difference:\t", finiteDifferenceDerivative(7, 0.44, knots, points, 3)

输出:

Derivatives: 
De Boor:              [10. 0.36134438  2.63969004]
Finite Difference:    [9.99999999 0.36134438 2.63969004]