User:Hdembinski/sandbox

In the mathematical subfield of numerical analysis de Boor's algorithm is a fast 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 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 B-spline curve $$ \mathbf{S}(x) $$ at position $$ x $$. The curve is build from a sum of potentially vector-valued constants $$ \mathbf{c}_i $$, called control points, and B-spline functions $$ B_{i,p}(x) $$,


 * $$ \mathbf{S}(x) = \sum_i \mathbf{c}_i B_{i,p}(x). $$

B-splines are piece-wise polynomial functions of order $$ 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.

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) := \left\{

\begin{matrix} 1 & \mathrm{if} \quad t_i \leq x < t_{i+1} \\ 0 & \mathrm{otherwise} \end{matrix} \right. $$


 * $$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 \in [k-p, 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,k-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 points 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 cell for $$ \mathbf{d}_{i,r-1} $$. Similarly, there is only one value of $$ \alpha $$ used in each step, so we can reuse the memory cell and drop the indices.

Furthermore, it is more convenient to use an index $$ j \in [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 $$ \mathbf{d}_{j} := \mathbf{c}_{j+k - p} $$ for $$ j = 0, \dots, p$$. Iterate for $$ r = 1, \dots, p $$:


 * $$ \alpha := \frac{x-t_{j + k - p}}{t_{j+1+k-r}-t_{j+k-p}}, $$


 * $$ \mathbf{d}_{j} := (1-\alpha) \mathbf{d}_{j-1} + \alpha \mathbf{d}_{j}; \quad j=p, \dots, r \quad $$ ($$ 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