De Boor's Algorithm

In the mathematical subfield of numerical analysis, de Boor's algorithmC. de Boor 971 "Subroutine package for calculating with B-splines", Techn.Rep. LA-4728-MS, Los Alamos Sci.Lab, Los Alamos NM; p. 109, 121. 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(x)$ at position $x$. The curve is built from a sum of B-spline functions $B_(x)$ multiplied with potentially vector-valued constants $\mathbf_i$, called control points,
$\mathbf(x) = \sum_i \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 $$ (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_(x)$ with $n = p + 1$.

Local support

B-splines have local support, meaning that the polynomials are positive only in a compact domain and zero elsewhere. The Cox-de Boor recursion formulaC. de Boor, p. 90 shows this:
$B_(x) := \begin 1 & \text \quad t_i \leq x < t_ \\ 0 & \text \end$
$B_(x) := \frac B_(x) + \frac B_(x).$

Let the index $k$ define the knot interval that contains the position, $x \in$ \mathbf(x) = \sum_^ \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_)$ 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_(x)$ directly. Instead it evaluates $\mathbf(x)$ through an equivalent recursion formula.

Let $\mathbf_$ be new control points with $\mathbf_ := \mathbf_$ for $i = k-p, \dots, k$. For $r = 1, \dots, p$ the following recursion is applied:
$\mathbf_ = (1-\alpha_) \mathbf_ + \alpha_ \mathbf_; \quad i=k-p+r,\dots,k$
$\alpha_ = \frac.$

Once the iterations are complete, we have $\mathbf(x) = \mathbf_$, meaning that $\mathbf_$ is the desired result.

De Boor's algorithm is more efficient than an explicit calculation of B-splines $B_(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_$. 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_$ reuse the memory for $\mathbf_$. 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_ := \mathbf_$ for $j = 0, \dots, p$. Iterate for $r = 1, \dots, p$:
$\mathbf_ := (1-\alpha_j) \mathbf_ + \alpha_j \mathbf_; \quad j=p, \dots, r \quad$
$\alpha_j := \frac.$
Note that must be counted down. After the iterations are complete, the result is $\mathbf(x) = \mathbf_$.

Example implementation

The following code in the Python (programming language)">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 - 1+ alpha * d
return d

See also

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

External links

De Boor's 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

References

Works cited
* {{cite book , author = Carl de Boor , title = A Practical Guide to Splines, Revised Edition , publisher = Springer-Verlag , year = 2003, isbn=0-387-95366-3

Articles with example Python (programming language) code
Numerical analysis
Splines (mathematics)
Interpolation