De Boor's algorithm: Difference between revisions

Revision as of 17:10, 10 August 2017

Revision as of 10:04, 11 August 2017

Line 1:

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>

{{Multiple Issues|

{{Cleanup|reason=The article does not differentiate between control points and knots. This leads to index errors in the formulas (usage of undefined knots <math>u_i</math>).|date=August 2017}}

{{Cleanup|reason=The meaning of the degrees of freedom (control points) of the spline <math>d_i</math> is distorted in the introduction.|date=August 2017}}

}}

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>

== 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''B-spline pointscurve <math>\mathbf{d}_0, \mathbf{dS}_1,(x) \dots, \mathbf{d}_{p-1}</math>, whichat playsposition the<math> rolex of</math>. aThe ''controlcurve is polygon''.build Thefrom curvea cansum beof describedpotentially asvector-valued a functionconstants <math> \mathbf{sc}(x)_i </math>, ofcalled one parameter ''x''. To pass through the sequence ofcontrol points, theand curveB-spline must satisfyfunctions <math>\mathbf B_{si,p}(u_0x)=\mathbf{d}_0, \dots</math>,

\mathbf{s}(u_{p-1})=\mathbf{d}_{p-1}</math>. But this is not quite the case: in general we are satisfied that the curve "approximates" the control polygon. We assume that ''u₀, ..., u_p-1'' are given to us along with <math>\mathbf{d}_0, \mathbf{d}_1, \dots, \mathbf{d}_{p-1}</math>.

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

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>).

DeB-splines Boor'sare algorithmpiece-wise ispolynomial anfunctions algorithmof which,order given ''u<submath>0</sub>, ...,p u<sub>p-1</submath>'' anddefined over a grid of knots <math>\mathbf {d}_0t_0, \mathbf{d}_1dots, t_i, \dots, \mathbf{d}_{p-1t_m} </math>, finds(we thealways valueuse ofzero-based splineindices curvein <math>\mathbf{s}(xthe following)</math>. at a pointDe Boor''x''.s Italgorithm uses [[Big O notation|O]](np²) + [[Big O notation|O]](n + p) operations where ''n'' is the degree and ''p'' the number of control points of ''s''.

== OutlineLocal ofsupport the algorithm==

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

We can express the curve as

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<ref>de Boor, p. 131</ref> shows this:

:<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-u_i}{u_{i+n}-u_i}N_i^{n-1}(x) + \frac{u_{i+n+1}-x}{u_{i+n+1}-u_{i+1}}N_{i+1}^{n-1}(x) ,</math>

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

and

\begin{matrix}

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

0 & \mathrm{otherwise}

\end{matrix}

\right.

</math>

:<math>N_i^0B_{i,p}(x) := \left\frac{\beginx - t_i}{matrixt_{i+p} 1,- &t_i} \mboxB_{if i,p-1}(x) \in+ [u_\frac{t_{i+p+1} - x},u_{t_{i+p+1}) \\- 0, & \mboxt_{otherwise i+1}} \endB_{matrixi+1,p-1}\right(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 \in [k-p, k] </math> are non-zero for this knot interval. Thus, the sum is reduced to:

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

Due to the spline locality property,

:<math> \mathbf{s}(x) = \sum_{i=\ell-n}^{\ell} \mathbf{d}_i N_i^n(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>.

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

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.

We can compute the above <math>\mathbf{s}(x)</math> by defining some <math> x \in [u_{\ell},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:

:Let <math> \mathbf{d}_i^_{[k]i,r} =</math> (1-\alpha_{k,i})be new control points with <math> \mathbf{d}_{i-1}^{[k-1],0} + \alpha_{k,i}:= \mathbf{dc}_i^_{[k-1]i}; \qquad</math> for <math> i k=1 k-p, \dots,n; \quadk</math>. For <math> r i=\ell-n+k 1, \dots,\ell p </math> the following recursion is applied:

:<math> \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 </math>

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

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

DoingOnce sothe givesiterations usare complete, we have <math>\mathbf{sS}(x) = \mathbf{d}_{k,p} </math>, meaning that <math> \ellmathbf{d}^_{[n]k,p} </math> is the desired result.

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.

== Optimizations ==

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 points 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 cell 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 cell and drop the indices.

Furthermore, it is more convenient to use an index <math> j \in [0,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:

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

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

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

== Example implementation ==

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

<source lang="python">

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

"""

Evaluates S(x).

Args

----

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]

</source>

== See also ==

Line 54 ⟶ 92:

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

== References ==

<references/>

'''Works cited'''

*{{cite book | author = Carl de Boor | title = A Practical Guide to Splines | publisher = Springer-Verlag | year = 1978|ISBN=3-540-90356-9}}

[[Category:Numerical analysis]]