De Boor's algorithm

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 $$ \mathbf{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 $$ \mathbf{c}_i $$, called control points, $$ \mathbf{S}(x) = \sum_i \mathbf{c}_i B_{i,p}(x). $$ B-splines of order $$ p + 1 $$ are connected piece-wise polynomial functions of degree $$ p $$ defined over a grid of knots $$ {t_0, \dots, t_i, \dots, t_m} $$ (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 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 shows this: $$B_{i,0}(x) := \begin{cases} 1 & \text{if } \quad t_i \leq x < t_{i+1} \\ 0 & \text{otherwise} \end{cases} $$ $$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). $$

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, \dots, k $$ are non-zero for this knot interval. Thus, the sum is reduced to: $$ \mathbf{S}(x) = \sum_{i=k-p}^{k} \mathbf{c}_i B_{i,p}(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 $$ \mathbf{S}(x) $$ through an equivalent recursion formula.

Let $$ \mathbf{d}_{i,r} $$ be new control points with $$ \mathbf{d}_{i,0} := \mathbf{c}_{i} $$ for $$ i = k-p, \dots, k$$. For $$ r = 1, \dots, p $$ the following recursion is applied: $$ \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 $$ $$ \alpha_{i,r} = \frac{x-t_i}{t_{i+1+p-r}-t_i}.$$

Once the iterations are complete, we have $$\mathbf{S}(x) = \mathbf{d}_{k,p} $$, meaning that $$ \mathbf{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 + \dots + 1 = (p + 1)(p + 2)/2 $$ temporary control points $$ \mathbf{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 $$ \mathbf{d}_{i,r} $$ reuse the memory for $$ \mathbf{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, \dots, p $$ for the temporary control points. The relation to the previous index is $$ i = j + k - p $$. Thus we obtain the improved algorithm:

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

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

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