De Boor's algorithm

(Redirected from De Boor algorithm)

In the mathematical subfield of numerical analysis, de Boor's algorithm[1] 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.[2][3]

Introduction

edit

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   at position  . The curve is built from a sum of B-spline functions   multiplied with potentially vector-valued constants  , called control points,   B-splines of order   are connected piece-wise polynomial functions of degree   defined over a grid of knots   (we always use zero-based indices in the following). De Boor's algorithm uses O(p2) + O(p) operations to evaluate the spline curve. Note: the main article about B-splines and the classic publications[1] use a different notation: the B-spline is indexed as   with  .

Local support

edit

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[4] shows this:    

Let the index   define the knot interval that contains the position,  . We can see in the recursion formula that only B-splines with   are non-zero for this knot interval. Thus, the sum is reduced to:  

It follows from   that  . Similarly, we see in the recursion that the highest queried knot location is at index  . This means that any knot interval   which is actually used must have at least   additional knots before and after. In a computer program, this is typically achieved by repeating the first and last used knot location   times. For example, for   and real knot locations  , one would pad the knot vector to  .

The algorithm

edit

With these definitions, we can now describe de Boor's algorithm. The algorithm does not compute the B-spline functions   directly. Instead it evaluates   through an equivalent recursion formula.

Let   be new control points with   for  . For   the following recursion is applied:    

Once the iterations are complete, we have  , meaning that   is the desired result.

De Boor's algorithm is more efficient than an explicit calculation of B-splines   with the Cox-de Boor recursion formula, because it does not compute terms which are guaranteed to be multiplied by zero.

Optimizations

edit

The algorithm above is not optimized for the implementation in a computer. It requires memory for   temporary control points  . Each temporary control point is written exactly once and read twice. By reversing the iteration over   (counting down instead of up), we can run the algorithm with memory for only   temporary control points, by letting   reuse the memory for  . Similarly, there is only one value of   used in each step, so we can reuse the memory as well.

Furthermore, it is more convenient to use a zero-based index   for the temporary control points. The relation to the previous index is  . Thus we obtain the improved algorithm:

Let   for  . Iterate for  :     Note that j must be counted down. After the iterations are complete, the result is  .

Example implementation

edit

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

edit
edit

Computer code

edit
  • 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

edit
  1. ^ a b C. de Boor [1971], "Subroutine package for calculating with B-splines", Techn.Rep. LA-4728-MS, Los Alamos Sci.Lab, Los Alamos NM; p. 109, 121.
  2. ^ Lee, E. T. Y. (December 1982). "A Simplified B-Spline Computation Routine". Computing. 29 (4). Springer-Verlag: 365–371. doi:10.1007/BF02246763. S2CID 2407104.
  3. ^ Lee, E. T. Y. (1986). "Comments on some B-spline algorithms". Computing. 36 (3). Springer-Verlag: 229–238. doi:10.1007/BF02240069. S2CID 7003455.
  4. ^ C. de Boor, p. 90

Works cited

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