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]
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 ) {\displaystyle \mathbf {S} (x)} at position x {\displaystyle x} . The curve is built from a sum of B-spline functions B i , p ( x ) {\displaystyle B_{i,p}(x)} multiplied with potentially vector-valued constants c i {\displaystyle \mathbf {c} _{i}} , called control points, S ( x ) = ∑ i c i B i , p ( x ) . {\displaystyle \mathbf {S} (x)=\sum _{i}\mathbf {c} _{i}B_{i,p}(x).} B-splines of order p + 1 {\displaystyle p+1} are connected piece-wise polynomial functions of degree p {\displaystyle p} defined over a grid of knots t 0 , … , t i , … , t m {\displaystyle {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[1] use a different notation: the B-spline is indexed as B i , n ( x ) {\displaystyle B_{i,n}(x)} with n = p + 1 {\displaystyle n=p+1} .
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 i , 0 ( x ) := { 1 if t i ≤ x < t i + 1 0 otherwise {\displaystyle 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 ) := x − t i t i + p − t i B i , p − 1 ( x ) + t i + p + 1 − x t i + p + 1 − t i + 1 B i + 1 , p − 1 ( x ) . {\displaystyle 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 {\displaystyle k} define the knot interval that contains the position, x ∈ [ t k , t k + 1 ) {\displaystyle x\in [t_{k},t_{k+1})} . We can see in the recursion formula that only B-splines with i = k − p , … , k {\displaystyle i=k-p,\dots ,k} are non-zero for this knot interval. Thus, the sum is reduced to: S ( x ) = ∑ i = k − p k c i B i , p ( x ) . {\displaystyle \mathbf {S} (x)=\sum _{i=k-p}^{k}\mathbf {c} _{i}B_{i,p}(x).}
It follows from i ≥ 0 {\displaystyle i\geq 0} that k ≥ p {\displaystyle k\geq p} . Similarly, we see in the recursion that the highest queried knot location is at index k + 1 + p {\displaystyle k+1+p} . This means that any knot interval [ t k , t k + 1 ) {\displaystyle [t_{k},t_{k+1})} which is actually used must have at least p {\displaystyle p} additional knots before and after. In a computer program, this is typically achieved by repeating the first and last used knot location p {\displaystyle p} times. For example, for p = 3 {\displaystyle p=3} and real knot locations ( 0 , 1 , 2 ) {\displaystyle (0,1,2)} , one would pad the knot vector to ( 0 , 0 , 0 , 0 , 1 , 2 , 2 , 2 , 2 ) {\displaystyle (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 ) {\displaystyle B_{i,p}(x)} directly. Instead it evaluates S ( x ) {\displaystyle \mathbf {S} (x)} through an equivalent recursion formula.
Let d i , r {\displaystyle \mathbf {d} _{i,r}} be new control points with d i , 0 := c i {\displaystyle \mathbf {d} _{i,0}:=\mathbf {c} _{i}} for i = k − p , … , k {\displaystyle i=k-p,\dots ,k} . For r = 1 , … , p {\displaystyle r=1,\dots ,p} the following recursion is applied: d i , r = ( 1 − α i , r ) d i − 1 , r − 1 + α i , r d i , r − 1 ; i = k − p + r , … , k {\displaystyle \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} α i , r = x − t i t i + 1 + p − r − t i . {\displaystyle \alpha _{i,r}={\frac {x-t_{i}}{t_{i+1+p-r}-t_{i}}}.}
Once the iterations are complete, we have S ( x ) = d k , p {\displaystyle \mathbf {S} (x)=\mathbf {d} _{k,p}} , meaning that d k , p {\displaystyle \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 ) {\displaystyle 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 + ⋯ + 1 = ( p + 1 ) ( p + 2 ) / 2 {\displaystyle (p+1)+p+\dots +1=(p+1)(p+2)/2} temporary control points d i , r {\displaystyle \mathbf {d} _{i,r}} . Each temporary control point is written exactly once and read twice. By reversing the iteration over i {\displaystyle i} (counting down instead of up), we can run the algorithm with memory for only p + 1 {\displaystyle p+1} temporary control points, by letting d i , r {\displaystyle \mathbf {d} _{i,r}} reuse the memory for d i , r − 1 {\displaystyle \mathbf {d} _{i,r-1}} . Similarly, there is only one value of α {\displaystyle \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 {\displaystyle j=0,\dots ,p} for the temporary control points. The relation to the previous index is i = j + k − p {\displaystyle i=j+k-p} . Thus we obtain the improved algorithm:
Let d j := c j + k − p {\displaystyle \mathbf {d} _{j}:=\mathbf {c} _{j+k-p}} for j = 0 , … , p {\displaystyle j=0,\dots ,p} . Iterate for r = 1 , … , p {\displaystyle r=1,\dots ,p} : d j := ( 1 − α j ) d j − 1 + α j d j ; j = p , … , r {\displaystyle \mathbf {d} _{j}:=(1-\alpha _{j})\mathbf {d} _{j-1}+\alpha _{j}\mathbf {d} _{j};\quad j=p,\dots ,r\quad } α j := x − t j + k − p t j + 1 + k − r − t j + k − p . {\displaystyle \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 S ( x ) = d p {\displaystyle \mathbf {S} (x)=\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]
Works cited