De Boor's algorithm explained

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

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

S (x) at position x . The curve is built from a sum of B-spline functions B _{i,p} (x) multiplied with potentially vector-valued constants c _{i} , called control points,$$\backslash mathbf(x)\; =\; \backslash sum\_i\; \backslash mathbf\_i\; B\_(x).$$B-splines of order p + 1 are connected piece-wise polynomial functions of degree p defined over a grid of knots {t _{0,} ... , t _{i,} ... , t _{m}} (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 B _{i,n} (x) with n = p + 1 .

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^[4] shows this:$$B\_(x)\; :=\; \backslash frac\; B\_(x)\; +\; \backslash frac\; B\_(x).$$

Let the index

k define the knot interval that contains the position, x \in [t _{k} ,t _{k+1} ) . We can see in the recursion formula that only B-splines with i = k - p, ... , k are non-zero for this knot interval. Thus, the sum is reduced to:$$\backslash mathbf(x)\; =\; \backslash sum\_^\; \backslash mathbf\_i\; B\_(x).$$

It follows from

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

The algorithm

With these definitions, we can now describe de Boor's algorithm. The algorithm does not compute the B-spline functions

B _{i,p} (x) directly. Instead it evaluates S (x) through an equivalent recursion formula.

Let

d _{i,r} be new control points with d _{i,0} : = c _{i} for i = k - p, ... , k . For r = 1 , ... , p the following recursion is applied:$$\backslash mathbf\_\; =\; (1-\backslash alpha\_)\; \backslash mathbf\_\; +\; \backslash alpha\_\; \backslash mathbf\_;\; \backslash quad\; i=k-p+r,\backslash dots,k$$$$\backslash alpha\_\; =\; \backslash frac.$$

Once the iterations are complete, we have

S (x) = d _{k,p} , meaning that d _{k,p} is the desired result.

De Boor's algorithm is more efficient than an explicit calculation of B-splines

B _{i,p} (x) 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

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

Furthermore, it is more convenient to use a zero-based index

j = 0 , ... , p for the temporary control points. The relation to the previous index is i = j + k - p . Thus we obtain the improved algorithm:

Let

d _{j} : = c _{j+k} for j = 0 , ... , p . Iterate for r = 1 , ... , p :$$\backslash mathbf\_\; :=\; (1-\backslash alpha\_j)\; \backslash mathbf\_\; +\; \backslash alpha\_j\; \backslash mathbf\_;\; \backslash quad\; j=p,\; \backslash dots,\; r\; \backslash quad$$$$\backslash alpha\_j\; :=\; \backslash frac.$$Note that must be counted down. After the iterations are complete, the result is S (x) = d _{p} .

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

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

Notes and References

1. 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 . 365–371 . Springer-Verlag . 10.1007/BF02246763. 2407104 .
3. Lee, E. T. Y. . Computing . 3 . 229–238 . Springer-Verlag . 10.1007/BF02240069. Comments on some B-spline algorithms . 36 . 1986. 7003455 .
4. C. de Boor, p. 90