De Boor's algorithm

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 $\backslash 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 $\backslash mathbf\{c\}\_i$, called control points,

$$\backslash mathbf\{S\}(x)\; =\; \backslash sum\_i\; \backslash 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,\; \backslash dots,\; t\_i,\; \backslash dots,\; 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$.

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\_\{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)\; :=\; \backslash frac\{x\; -\; t\_i\}\{t\_\{i+p\}\; -\; t\_i\}\; B\_\{i,p-1\}(x)\; +\; \backslash 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\; \backslash in\; [t\_\{k\},t\_\{k+1\})$. We can see in the recursion formula that only B-splines with $i\; =\; k-p,\; \backslash dots,\; k$ are non-zero for this knot interval. Thus, the sum is reduced to:

$$\backslash mathbf\{S\}(x)\; =\; \backslash sum\_\{i=k-p\}^\{k\}\; \backslash mathbf\{c\}\_i\; B\_\{i,p\}(x).$$

It follows from $i\; \backslash geq\; 0$ that $k\; \backslash 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)$.

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 $\backslash mathbf\{S\}(x)$ through an equivalent recursion formula.

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

$$\backslash mathbf\{d\}\_\{i,r\}\; =\; (1-\backslash alpha\_\{i,r\})\; \backslash mathbf\{d\}\_\{i-1,r-1\}\; +\; \backslash alpha\_\{i,r\}\; \backslash mathbf\{d\}\_\{i,r-1\};\; \backslash quad\; i=k-p+r,\backslash dots,k$$

$$\backslash alpha\_\{i,r\}\; =\; \backslash frac\{x-t\_i\}\{t\_\{i+1+p-r\}-t\_i\}.$$

Once the iterations are complete, we have $\backslash mathbf\{S\}(x)\; =\; \backslash mathbf\{d\}\_\{k,p\}$, meaning that $\backslash 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.

The algorithm above is not optimized for the implementation in a computer. It requires memory for $(p\; +\; 1)\; +\; p\; +\; \backslash dots\; +\; 1\; =\; (p\; +\; 1)(p\; +\; 2)/2$ temporary control points $\backslash 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 $\backslash mathbf\{d\}\_\{i,r\}$ reuse the memory for $\backslash mathbf\{d\}\_\{i,r-1\}$. Similarly, there is only one value of $\backslash 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,\; \backslash 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 $\backslash mathbf\{d\}\_\{j\}\; :=\; \backslash mathbf\{c\}\_\{j+k\; -\; p\}$ for $j\; =\; 0,\; \backslash dots,\; p$. Iterate for $r\; =\; 1,\; \backslash dots,\; p$:

$$\backslash mathbf\{d\}\_\{j\}\; :=\; (1-\backslash alpha\_j)\; \backslash mathbf\{d\}\_\{j-1\}\; +\; \backslash alpha\_j\; \backslash mathbf\{d\}\_\{j\};\; \backslash quad\; j=p,\; \backslash dots,\; r\; \backslash quad$$

$$\backslash alpha\_j\; :=\; \backslash frac\{x-t\_\{j\; +\; k\; -\; p\}\}\{t\_\{j+1+k-r\}-t\_\{j+k-p\}\}.$$

Note that {{mvar|j}} must be counted down. After the iterations are complete, the result is $\backslash mathbf\{S\}(x)\; =\; \backslash mathbf\{d\}\_\{p\}$.

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]

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

• [http://www.cs.mtu.edu/~shene/COURSES/cs3621/NOTES/spline/de-Boor.html De Boor's Algorithm]
• [http://www.idav.ucdavis.edu/education/CAGDNotes/Deboor-Cox-Calculation/Deboor-Cox-Calculation.html The DeBoor-Cox Calculation]

• [http://www.netlib.org/pppack/ PPPACK]: contains many spline algorithms in Fortran
• [https://www.gnu.org/software/gsl/ GNU Scientific Library]: C-library, contains a sub-library for splines ported from PPPACK
• [https://www.scipy.org/ SciPy]: Python-library, contains a sub-library scipy.interpolate with spline functions based on FITPACK
• [https://github.com/msteinbeck/tinyspline TinySpline]: C-library for splines with a C++ wrapper and bindings for C#, Java, Lua, PHP, Python, and Ruby
• [http://einspline.sourceforge.net/ Einspline]: C-library for splines in 1, 2, and 3 dimensions with Fortran wrappers

{{reflist}}

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}}

Category:Articles with example Python (programming language) code

Category:Numerical analysis

Category:Splines (mathematics)

Category:Interpolation