De Boor's algorithm: Difference between revisions

Revision as of 10:22, 4 August 2017

Latest revision as of 18:55, 9 August 2025

(40 intermediate revisions by 20 users not shown)

Line 1:

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

In the [[mathematics|mathematical]] subfield of [[numerical analysis]] '''de Boor's algorithm''' is a fast 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}}</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}}</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 fast[[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 ==

TheA general settingintroduction to B-splines is asgiven in the [[B-spline|main followsarticle]]. WeHere wouldwe likediscuss tode constructBoor's aalgorithm, curvean whoseefficient shapeand isnumerically describedstable byscheme ato sequenceevaluate ofa ''p''spline pointscurve <math>\mathbf{d}_0, \mathbf{dS}_1,(x) \dots, \mathbf{d}_{p-1}</math>, whichat playsposition the<math> rolex of a ''control polygon''</math>. The curve can beis describedbuilt asfrom a functionsum of B-spline functions <math> \mathbfB_{si,p}(x) </math> ofmultiplied onewith parameterpotentially ''x''. To pass through the sequence of points, the curve mustvector-valued satisfyconstants <math> \mathbf{sc}(u_0)=\mathbf{d}_0_i </math>, \dotscalled control points,

:<math display="block"> \mathbf{sS}(x) = \sum_{i=0}^{p-1}sum_i \mathbf{dc}_i N_i^nB_{i,p}(x) ,. </math>

\mathbf{s}(u_{p-1})=\mathbf{d}_{p-1}</math> for a given knot sequence <math>u_0<u_1<\ldots<u_{p-1}</math>. But this is not quite the case: in general we accept some trade-off between the "approximation" <math>\mathbf{s}(u_k)\approx \mathbf{d}_k</math> of the control points and the smoothness of the curve <math>s(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 ==

One approach to solve this problem is by [[spline (mathematics)|spline]]s. A spline is a curve that is a piecewise ''n^th'' degree polynomial. This means that, on any interval ''<nowiki>[</nowiki>u_i, u_i+1)'', the curve must be equal to a polynomial of degree at most ''n''. It may be equal to different polynomials on different intervals. The polynomials must be ''synchronized'': when the polynomials from intervals ''<nowiki>[</nowiki>u_i-1, u_i)'' and ''<nowiki>[</nowiki>u_i, u_i+1)'' meet at the point ''u_i'', they must have the same value at this point and their derivatives must be equal up to order <math>n-1</math> (to ensure that the curve is as smooth as possible without restricting <math>s(x)</math> to a plain polynomial within <math>[u_{i-1},u_{i+1})</math>).

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

De Boor's algorithm is an algorithm which, given ''u₀, ..., u_p-1'' and <math>\mathbf{d}_0, \mathbf{d}_1, \dots, \mathbf{d}_{p-1}</math>, finds the value of spline curve <math>\mathbf{s}(x)</math> at a point ''x''. It uses [[Big O notation|O]](n²) + [[Big O notation|O]](n + p) operations where ''n'' is the degree and ''p'' the number of control points of ''s''.

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

\begin{cases}

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

0 & \text{otherwise}

\end{cases}

</math>

:<math display="block">N_i^nB_{i,p}(x) := \frac{x -\bar u_it_i}{\bar u_t_{i+np} -\bar u_it_i}N_i^ B_{ni,p-1}(x) + \frac{\bar u_t_{i+np+1} - x}{\bar u_t_{i+np+1} -\bar u_t_{i+1}}N_ B_{i+1}^{n,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:

== Outline of the algorithm==

:<math display="block"> \mathbf{sS}(x) = \sum_{i=\ellk-np}^{\ellk} \mathbf{dc}_i N_i^nB_{i,p}(x). </math>

Suppose we want to evaluate the spline curve for a parameter value <math> x \in [u_{0},u_{p-1}] </math>.

We can express the curve as

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=0}^{p-1} \mathbf{d}_i N_i^n(x) , </math>

where<ref>http://www.cs.mtu.edu/~shene/COURSES/cs3621/NOTES/spline/B-spline/bspline-basis.html</ref>

:<math>N_i^n(x)=\frac{x-\bar u_i}{\bar u_{i+n}-\bar u_i}N_i^{n-1}(x) + \frac{\bar u_{i+n+1}-x}{\bar u_{i+n+1}-\bar u_{i+1}}N_{i+1}^{n-1}(x) ,</math>

== The algorithm ==

:<math>N_i^0(x)=\left\{\begin{matrix} 1, & \mbox{if }x \in [\bar u_{i},\bar u_{i+1}) \\ 0, & \mbox{otherwise } \end{matrix}\right.</math>

With these definitions, we can now describe de Boor's algorithm. The algorithm does not compute the B-spline functions <math> B_{i,p}(x) </math> directly. Instead it evaluates <math> \mathbf{S}(x) </math> through an equivalent recursion formula.

with the extended knot sequence

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>\bar u_i := \left\{\begin{matrix} u_0& \mbox{for }i=0,\ldots,n\\

:<math display="block"> \mathbf{d}_i^_{[k]i,r} = (1-\alpha_{ki,ir}) \mathbf{d}_{i-1}^{[k,r-1]} + \alpha_{ki,ir} \mathbf{d}_i^_{[ki,r-1]}; \qquad k=1,\dots,n; \quad i=\ellk-np+kr,\dots,\ellk </math>

u_{i-n}&\mbox{for }i=n+1,\ldots,p+n-1\\

:<math display="block"> \alpha_{k,i,r} = \frac{x-\bar u_it_i}{\bar u_t_{i+n1+1p-kr}-\bar u_it_i}.</math>

u_{p-1}&\mbox{for }i=p+n,\ldots,p+2n.

\end{matrix}\right.</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.

While evaluating the algorithm quotients with vanishing numerator and denominator have to be set to zero.

De Boor's algorithm is more efficient than an explicit calculation of B-splines <math> B_{i,p}(x) </math> with the Cox-de Boor recursion formula, because it does not compute terms which are guaranteed to be multiplied by zero.

Let <math>\ell\in\{n,\ldots,p+n-2\}</math> be the index such that <math>x\in[\bar u_\ell,\bar u_{\ell+1})</math>. Due to the spline locality property,

:<math> \mathbf{s}(x) = \sum_{i=\ell-n}^{\ell} \mathbf{d}_i N_i^n(x) </math>

So the value <math>\mathbf{s}(x)</math> is determined by the control points <math> \mathbf{d}_{\ell-n},\mathbf{d}_{\ell-n+1},\dots,\mathbf{d}_{\ell} </math>; the other control points <math>\mathbf{d}_i</math> have no influence. De Boor's algorithm, described in the next section, is a procedure which efficiently calculates the expression for <math> \mathbf{s}(x) </math>.

== The algorithmOptimizations ==

The algorithm above is not optimized for the implementation in a computer. It requires memory for <math> (p + 1) + p + \dots + 1 = (p + 1)(p + 2)/2 </math> temporary control points <math> \mathbf{d}_{i,r} </math>. Each temporary control point is written exactly once and read twice. By reversing the iteration over <math> i </math> (counting down instead of up), we can run the algorithm with memory for only <math> p + 1 </math> temporary control points, by letting <math> \mathbf{d}_{i,r} </math> reuse the memory for <math> \mathbf{d}_{i,r-1} </math>. Similarly, there is only one value of <math> \alpha </math> used in each step, so we can reuse the memory as well.

We can compute the above <math>\mathbf{s}(x)</math> by defining some <math> x \in [\bar u_{\ell},\bar u_{\ell+1}) </math>, setting <math> \mathbf{d}_i^{[0]} = \mathbf{d}_i</math> for <math>i = \ell-n, \dots, \ell</math>, and with these, computing:

Furthermore, it is more convenient to use a zero-based index <math> j = 0, \dots, p </math> for the temporary control points. The relation to the previous index is <math> i = j + k - p </math>. Thus we obtain the improved algorithm:

:<math> \mathbf{d}_i^{[k]} = (1-\alpha_{k,i}) \mathbf{d}_{i-1}^{[k-1]} + \alpha_{k,i} \mathbf{d}_i^{[k-1]}; \qquad k=1,\dots,n; \quad i=\ell-n+k,\dots,\ell </math>

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

Where the ratio <math>\alpha</math> is described by:

<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_j := \frac{x-t_{j + k - p}}{t_{j+1+k-r}-t_{j+k-p}}. </math>

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

== Example implementation ==

:<math> \alpha_{k,i} = \frac{x-\bar u_i}{\bar u_{i+n+1-k}-\bar u_i}</math>

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

Doing so gives us <math>\mathbf{s}(x) = \mathbf{d}_{\ell}^{[n]} </math>

<syntaxhighlight lang="python" line>

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]

</syntaxhighlight>

== See also ==

* [[De Casteljau's algorithm]]

* [[Bézier curve]]

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

*[[NURBS]]

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

* [http://www.netlib.org/pppack/ PPPACK]: contains many spline algorithms in [[Fortran]]

[https://github.com/msteinbeck/tinyspline TinySpline: Open source C-library for splines which implements De Boor's algorithm]

* [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]: Open source C-library for splines whichwith implementsa DeC++ Boor'swrapper algorithmand bindings for C#, Java, Lua, [[PHP]], Python, and Ruby

* [https://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)]]