B-Spline derivative using de Boor's algorithm

Question

Wikipedia gives us a Python implementation for the de Boor's algorithm:

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]

Is there a similar algorithm calculating the derivative of the B-Spline interpolated curve (or even n-th derivative)?

I know that mathematically it is reduced to using a spline of the lower order but can't apply it to the de Boor's algorithm.

Answer 1

I think I found the right way to re-use the de Boor's algorithm for curve derivatives.

First, we consider the definition of the B-Spline curve. It is a linear combination of control points:       (1)

Hence, the derivative is a linear combination of the basis-function derivatives

      (2)

The derivative of the basis function is defined as follows:

      (3)

We plug-in (3) into (2) and after some algebra kung-fu, described here http://public.vrac.iastate.edu/~oliver/courses/me625/week5b.pdf, we obtain:

      (4), where

The derivative of the B-Spline curve is nothing else but a new B-Spline curve of (p-1) degree built on top of the new control points Q. Now, to employ the de Boor's algorithm we compute the new control point set and lower the spline degree p by 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]

Test:

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)

Output:

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

Answer 2

If you want to calculate the value and its derivative at the same time, then it's reasonable to just replace all the values v that depend on x with tuples (v, dv/dx)

Then you can just apply the sum and product rules when you add or multiply them: https://en.wikipedia.org/wiki/Product_rule

The function you provided would turn into this, for example:

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]