De Boor's algorithm

Short description: Method of evaluating spline curves

In the mathematical subfield of numerical analysis, de Boor's algorithm is a polynomial-time and numerically stable algorithm for evaluating spline curves in B-spline form. It is a generalization of de Casteljau's algorithm for Bézier curves. The algorithm was devised by German-American mathematician Carl R. de Boor. Simplified, potentially faster variants of the de Boor algorithm have been created but they suffer from comparatively lower stability.

Introduction

A general introduction to B-splines is given in the main article. Here we discuss de Boor's algorithm, an efficient and numerically stable scheme to evaluate a spline curve [math]\displaystyle{ \mathbf{S}(x) }[/math] at position [math]\displaystyle{ x }[/math]. The curve is built from a sum of B-spline functions [math]\displaystyle{ B_{i,p}(x) }[/math] multiplied with potentially vector-valued constants [math]\displaystyle{ \mathbf{c}_i }[/math], called control points,

[math]\displaystyle{ \mathbf{S}(x) = \sum_i \mathbf{c}_i B_{i,p}(x). }[/math]

B-splines of order [math]\displaystyle{ p + 1 }[/math] are connected piece-wise polynomial functions of degree [math]\displaystyle{ p }[/math] defined over a grid of knots [math]\displaystyle{ {t_0, \dots, t_i, \dots, t_m} }[/math] (we always use zero-based indices in the following). De Boor's algorithm uses O(p²) + O(p) operations to evaluate the spline curve. Note: the main article about B-splines and the classic publications use a different notation: the B-spline is indexed as [math]\displaystyle{ B_{i,n}(x) }[/math] with [math]\displaystyle{ n = p + 1 }[/math].

Local support

B-splines have local support, meaning that the polynomials are positive only in a finite domain and zero elsewhere. The Cox-de Boor recursion formula shows this:

[math]\displaystyle{ B_{i,0}(x) := \begin{cases} 1 & \text{if } \quad t_i \leq x \lt t_{i+1} \\ 0 & \text{otherwise} \end{cases} }[/math]

[math]\displaystyle{ B_{i,p}(x) := \frac{x - t_i}{t_{i+p} - t_i} B_{i,p-1}(x) + \frac{t_{i+p+1} - x}{t_{i+p+1} - t_{i+1}} B_{i+1,p-1}(x). }[/math]

Let the index [math]\displaystyle{ k }[/math] define the knot interval that contains the position, [math]\displaystyle{ x \in [t_{k},t_{k+1}) }[/math]. We can see in the recursion formula that only B-splines with [math]\displaystyle{ i = k-p, \dots, k }[/math] are non-zero for this knot interval. Thus, the sum is reduced to:

[math]\displaystyle{ \mathbf{S}(x) = \sum_{i=k-p}^{k} \mathbf{c}_i B_{i,p}(x). }[/math]

It follows from [math]\displaystyle{ i \geq 0 }[/math] that [math]\displaystyle{ k \geq p }[/math]. Similarly, we see in the recursion that the highest queried knot location is at index [math]\displaystyle{ k + 1 + p }[/math]. This means that any knot interval [math]\displaystyle{ [t_k,t_{k+1}) }[/math] which is actually used must have at least [math]\displaystyle{ p }[/math] additional knots before and after. In a computer program, this is typically achieved by repeating the first and last used knot location [math]\displaystyle{ p }[/math] times. For example, for [math]\displaystyle{ p = 3 }[/math] and real knot locations [math]\displaystyle{ (0, 1, 2) }[/math], one would pad the knot vector to [math]\displaystyle{ (0,0,0,0,1,2,2,2,2) }[/math].

The algorithm

With these definitions, we can now describe de Boor's algorithm. The algorithm does not compute the B-spline functions [math]\displaystyle{ B_{i,p}(x) }[/math] directly. Instead it evaluates [math]\displaystyle{ \mathbf{S}(x) }[/math] through an equivalent recursion formula.

Let [math]\displaystyle{ \mathbf{d}_{i,r} }[/math] be new control points with [math]\displaystyle{ \mathbf{d}_{i,0} := \mathbf{c}_{i} }[/math] for [math]\displaystyle{ i = k-p, \dots, k }[/math]. For [math]\displaystyle{ r = 1, \dots, p }[/math] the following recursion is applied:

[math]\displaystyle{ \mathbf{d}_{i,r} = (1-\alpha_{i,r}) \mathbf{d}_{i-1,r-1} + \alpha_{i,r} \mathbf{d}_{i,r-1}; \quad i=k-p+r,\dots,k }[/math]

[math]\displaystyle{ \alpha_{i,r} = \frac{x-t_i}{t_{i+1+p-r}-t_i}. }[/math]

Once the iterations are complete, we have [math]\displaystyle{ \mathbf{S}(x) = \mathbf{d}_{k,p} }[/math], meaning that [math]\displaystyle{ \mathbf{d}_{k,p} }[/math] is the desired result.

De Boor's algorithm is more efficient than an explicit calculation of B-splines [math]\displaystyle{ B_{i,p}(x) }[/math] with the Cox-de Boor recursion formula, because it does not compute terms which are guaranteed to be multiplied by zero.

Optimizations

The algorithm above is not optimized for the implementation in a computer. It requires memory for [math]\displaystyle{ (p + 1) + p + \dots + 1 = (p + 1)(p + 2)/2 }[/math] temporary control points [math]\displaystyle{ \mathbf{d}_{i,r} }[/math]. Each temporary control point is written exactly once and read twice. By reversing the iteration over [math]\displaystyle{ i }[/math] (counting down instead of up), we can run the algorithm with memory for only [math]\displaystyle{ p + 1 }[/math] temporary control points, by letting [math]\displaystyle{ \mathbf{d}_{i,r} }[/math] reuse the memory for [math]\displaystyle{ \mathbf{d}_{i,r-1} }[/math]. Similarly, there is only one value of [math]\displaystyle{ \alpha }[/math] used in each step, so we can reuse the memory as well.

Furthermore, it is more convenient to use a zero-based index [math]\displaystyle{ j = 0, \dots, p }[/math] for the temporary control points. The relation to the previous index is [math]\displaystyle{ i = j + k - p }[/math]. Thus we obtain the improved algorithm:

Let [math]\displaystyle{ \mathbf{d}_{j} := \mathbf{c}_{j+k - p} }[/math] for [math]\displaystyle{ j = 0, \dots, p }[/math]. Iterate for [math]\displaystyle{ r = 1, \dots, p }[/math]:

[math]\displaystyle{ \mathbf{d}_{j} := (1-\alpha_j) \mathbf{d}_{j-1} + \alpha_j \mathbf{d}_{j}; \quad j=p, \dots, r \quad }[/math] ([math]\displaystyle{ j }[/math] must be counted down)

[math]\displaystyle{ \alpha_j := \frac{x-t_{j + k - p}}{t_{j+1+k-r}-t_{j+k-p}}. }[/math]

After the iterations are complete, the result is [math]\displaystyle{ \mathbf{S}(x) = \mathbf{d}_{p} }[/math].

Example implementation

The following code in the Python programming language is a naive implementation of the optimized algorithm.

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

    Arguments
    ---------
    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]

See also

• De Casteljau's algorithm
• Bézier curve
• Non-uniform rational B-spline

External links

• De Boor's Algorithm
• The DeBoor-Cox Calculation

Computer code

• PPPACK: contains many spline algorithms in Fortran
• GNU Scientific Library: C-library, contains a sub-library for splines ported from PPPACK
• SciPy: Python-library, contains a sub-library scipy.interpolate with spline functions based on FITPACK
• TinySpline: C-library for splines with a C++ wrapper and bindings for C#, Java, Lua, PHP, Python, and Ruby
• Einspline: C-library for splines in 1, 2, and 3 dimensions with Fortran wrappers

References

Works cited

• Carl de Boor (2003). A Practical Guide to Splines, Revised Edition. Springer-Verlag. ISBN 0-387-95366-3.