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

Question

维基百科为我们提供了德布尔算法的 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-Spline插值曲线的导数（甚至n阶导数）是否有类似的算法计算？

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

Answer 1

如果你想同时计算值和它的导数，那么只需要替换所有依赖于x[的值v是合理的=21=] 与元组 (v, dv/dx)

然后您可以在加法或乘法时应用加法和乘法规则：https://en.wikipedia.org/wiki/Product_rule

你提供的函数会变成这样，例如：

def deBoorWithDerivative(k, x, t, c, p): """ Evaluates (S(x), dS(x)/dx). 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],0) for j in range(0, p+1)] for r in range(1, p+1): for j in range(p, r-1, -1): dalpha = 1.0/(t[j+1+k-r] - t[j+k-p]) alpha = (x - t[j+k-p]) * dalpha d[j] = ( (1.0 - alpha) * d[j-1][0] + alpha * d[j][0], -dalpha * d[j-1][0] + (1.0 - alpha) * d[j-1][1] +dalpha * d[j][0] + alpha*d[j][1] ) return d[p]

Answer 2

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

首先，我们考虑B样条曲线的定义。它是控制点的线性组合： $C(u) = \sum_{i=0}^{n}N_{i,p}(u)P_i$ (1)

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

$C'(u) = \sum_{i=0}^{n}N'_{i,p}(u)P_i$ (2)

基函数的导数定义如下：

$N'(u) = \frac{p}{u_{i+p} - u_i} N_{i,p-1} - \frac{p}{u_{i+p+1}-u_{i+1}}N_{i+1,p-1}$ (3)

我们将 (3) 插入到 (2) 中，经过这里 http://public.vrac.iastate.edu/~oliver/courses/me625/week5b.pdf 描述的一些代数功夫后，我们得到：

$C'$u$ = \sum_{i=0}^{n-1}N_{i,p-1}$u$Q_i$ (4)，其中 $Q_i = p\frac{P_{i+1}-P{i}}{u_{i+p+1} - u_{i+1}}$

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]

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

B-Spline derivative using de Boor's algorithm

python

algorithm

interpolation

bspline