De Boor's algorithm: Difference between revisions

Revision as of 13:10, 31 May 2023

Latest revision as of 18:55, 9 August 2025

(5 intermediate revisions by 5 users not shown)

Line 1:

{{short description|Method of evaluating spline curves}}

In the [[mathematics|mathematical]] subfield of [[numerical analysis]] '''de Boor's algorithm'''<ref name="de_boor_paper">C. de Boor [1971], "Subroutine package for calculating with B-splines", Techn.Rep. LA-4728-MS, Los Alamos Sci.Lab, Los Alamos NM; p. 109, 121.</ref> is a polynomial-time and [[numerically stable]] [[algorithm]] for evaluating [[spline curve]]s in [[B-spline]] form. It is a generalization of [[de Casteljau's algorithm]] for [[Bézier curve]]s. 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.<ref>{{cite journal |last=Lee |first=E. T. Y. |date=December 1982 |title=A Simplified B-Spline Computation Routine |journal=Computing |volume=29 |issue=4 |pages=365–371 |publisher=Springer-Verlag|doi=10.1007/BF02246763|s2cid=2407104 }}</ref><ref>{{cite journal | author = Lee, E. T. Y. | journal = Computing | issue = 3 | pages = 229–238 | publisher = Springer-Verlag | doi=10.1007/BF02240069|title = Comments on some B-spline algorithms | volume = 36 | year = 1986| s2cid = 7003455 }}</ref>

In the [[mathematics|mathematical]] subfield of [[numerical analysis]], '''de Boor's algorithm'''<ref name="de_boor_paper">C. de Boor [1971], "Subroutine package for calculating with B-splines", Techn.Rep. LA-4728-MS, Los Alamos Sci.Lab, Los Alamos NM; p. 109, 121.</ref> is a [[polynomial-time]] and [[numerically stable]] [[algorithm]] for evaluating [[spline curve]]s in [[B-spline]] form. It is a generalization of [[de Casteljau's algorithm]] for [[Bézier curve]]s. 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.<ref>{{cite journal |last=Lee |first=E. T. Y. |date=December 1982 |title=A Simplified B-Spline Computation Routine |journal=Computing |volume=29 |issue=4 |pages=365–371 |publisher=Springer-Verlag | doi=10.1007/BF02246763|s2cid=2407104 }}</ref><ref>{{cite journal | author = Lee, E. T. Y. | journal = Computing | issue = 3 | pages = 229–238 | publisher = Springer-Verlag | doi=10.1007/BF02240069|title = Comments on some B-spline algorithms | volume = 36 | year = 1986| s2cid = 7003455 }}</ref>

== Introduction ==

A general introduction to B-splines is given in the [[B-spline|main article]]. Here we discuss de Boor's algorithm, an efficient and numerically stable scheme to evaluate a spline curve <math> \mathbf{S}(x) </math> at position <math> x </math>. The curve is built from a sum of B-spline functions <math> B_{i,p}(x) </math> multiplied with potentially vector-valued constants <math> \mathbf{c}_i </math>, called control points,

:<math display="block"> \mathbf{S}(x) = \sum_i \mathbf{c}_i B_{i,p}(x). </math>

:<math> \mathbf{S}(x) = \sum_i \mathbf{c}_i B_{i,p}(x). </math>

B-splines of order <math> p + 1 </math> are connected piece-wise polynomial functions of degree <math> p </math> defined over a grid of knots <math> {t_0, \dots, t_i, \dots, t_m} </math> (we always use zero-based indices in the following). De Boor's algorithm uses [[Big O notation|O]](p²) + [[Big O notation|O]](p) operations to evaluate the spline curve. Note: the [[B-spline|main article about B-splines]] and the classic publications<ref name="de_boor_paper"></ref> use a different notation: the B-spline is indexed as <math> B_{i,n}(x) </math> with <math>n = p + 1</math>.

== Local support ==

B-splines have local support, meaning that the polynomials are positive only in a finitecompact ___domain and zero elsewhere. The Cox-de Boor recursion formula<ref>C. de Boor, p. 90</ref> shows this:

:<math display="block">B_{i,0}(x) :=

:<math>B_{i,0}(x) :=

\begin{cases}

1 & \text{if } \quad t_i \leq x < t_{i+1} \\

Line 19 ⟶ 18:

\end{cases}

</math>

:<math display="block">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). </math>

:<math>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). </math>

Let the index <math> k </math> define the knot interval that contains the position, <math> x \in [t_{k},t_{k+1}) </math>. We can see in the recursion formula that only B-splines with <math> i = k-p, \dots, k </math> are non-zero for this knot interval. Thus, the sum is reduced to:

:<math display="block"> \mathbf{S}(x) = \sum_{i=k-p}^{k} \mathbf{c}_i B_{i,p}(x). </math>

It follows from <math> i \geq 0 </math> that <math> k \geq p </math>. Similarly, we see in the recursion that the highest queried knot ___location is at index <math> k + 1 + p </math>. This means that any knot interval <math> [t_k,t_{k+1}) </math> which is actually used must have at least <math> p </math> additional knots before and after. In a [[computer program]], this is typically achieved by repeating the first and last used knot ___location <math> p </math> times. For example, for <math> p = 3 </math> and real knot locations <math> (0, 1, 2) </math>, one would pad the knot vector to <math> (0,0,0,0,1,2,2,2,2) </math>.

:<math> \mathbf{S}(x) = \sum_{i=k-p}^{k} \mathbf{c}_i B_{i,p}(x). </math>

It follows from <math> i \geq 0 </math> that <math> k \geq p </math>. Similarly, we see in the recursion that the highest queried knot ___location is at index <math> k + 1 + p </math>. This means that any knot interval <math> [t_k,t_{k+1}) </math> which is actually used must have at least <math> p </math> additional knots before and after. In a computer program, this is typically achieved by repeating the first and last used knot ___location <math> p </math> times. For example, for <math> p = 3 </math> and real knot locations <math> (0, 1, 2) </math>, one would pad the knot vector to <math> (0,0,0,0,1,2,2,2,2) </math>.

== The algorithm ==

Line 33 ⟶ 30:

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

:<math display="block"> \mathbf{d}_{ji,r} := (1-\alpha_jalpha_{i,r}) \mathbf{d}_{ji-1,r-1} + \alpha_jalpha_{i,r} \mathbf{d}_{ji,r-1}; \quad ji=k-p+r, \dots, r \quadk </math> (<math> j </math> must be counted down)

:<math display="block"> \mathbf{d}_alpha_{i,r} = (1-\alpha_frac{i,rx-t_i}) \mathbf{d}_t_{i-+1,r+p-1} + \alpha_{i,r} \mathbf{d}_{i,r-1t_i}; \quad i=k-p+r,\dots,k .</math>

:<math> \alpha_{i,r} = \frac{x-t_i}{t_{i+1+p-r}-t_i}.</math>

Once the iterations are complete, we have <math>\mathbf{S}(x) = \mathbf{d}_{k,p} </math>, meaning that <math> \mathbf{d}_{k,p} </math> is the desired result.

Line 49 ⟶ 44:

Let <math> \mathbf{d}_{j} := \mathbf{c}_{j+k - p} </math> for <math> j = 0, \dots, p</math>. Iterate for <math> r = 1, \dots, p </math>:

<math display="block"> \mathbf{d}_{j} := (1-\alpha_j) \mathbf{d}_{j-1} + \alpha_j \mathbf{d}_{j}; \quad j=p, \dots, r \quad </math>

:<math display="block"> \alpha_{i,r}alpha_j := \frac{x-t_it_{j + k - p}}{t_{ij+1+pk-r}-t_it_{j+k-p}}. </math>

:<math> \mathbf{d}_{j} := (1-\alpha_j) \mathbf{d}_{j-1} + \alpha_j \mathbf{d}_{j}; \quad j=p, \dots, r \quad </math> (<math> j </math> must be counted down)

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

:<math> \alpha_j := \frac{x-t_{j + k - p}}{t_{j+1+k-r}-t_{j+k-p}}. </math>

After the iterations are complete, the result is <math>\mathbf{S}(x) = \mathbf{d}_{p} </math>.

== Example implementation ==

Line 60 ⟶ 52:

The following code in the [[Python (programming language)|Python programming language]] is a naive implementation of the optimized algorithm.

<syntaxhighlight lang="python" line>

def deBoor(k: int, x: int, t, c, p: int):

"""Evaluates S(x).

Line 80 ⟶ 72:

return d[p]

</syntaxhighlight>

== See also ==

* [[De Casteljau's algorithm]]

* [[Bézier curve]]

* [[Non-uniform rational B-spline]]

== External links ==

* [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]

== Computer code ==

Line 96 ⟶ 87:

* [https://www.gnu.org/software/gsl/ GNU Scientific Library]: C-library, contains a sub-library for splines ported from [[Netlib|PPPACK]]

* [https://www.scipy.org/ SciPy]: Python-library, contains a sub-library ''scipy.interpolate'' with spline functions based on [[Netlib|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

* [httphttps://einspline.sourceforge.net/ Einspline]: C-library for splines in 1, 2, and 3 dimensions with Fortran wrappers

== References ==

{{reflist}}

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

[[Category:Articles with example Python (programming language) code]]

[[Category:Numerical analysis]]

[[Category:Splines (mathematics)]]