Algoritmo de De Boor

En el subcampo matemático del análisis numérico, el algoritmo de De Boor^[1]​ es un algoritmo de tiempo polinomial y numéricamente estable para evaluar curvas spline en forma B-spline. Es una generalización del algoritmo Casteljau para las curvas de Bézier. El algoritmo fue ideado por Carl R. De Boor. Se han creado variantes simplificadas y potencialmente más rápidas del algoritmo de De Boor, pero sufren una estabilidad comparativamente menor.^[2]​^[3]​

El algoritmo de De Boor es un esquema eficiente y numéricamente estable para evaluar una curva spline ${\displaystyle \mathbf {S} (x)}$ en posición ${\displaystyle x}$. La curva se construye a partir de una suma de funciones B-spline ${\displaystyle B_{i,p}(x)}$ multiplicada con valores vectoriales potencialmente constantes ${\displaystyle \mathbf {c} _{i}}$, llamados puntos de control.

${\displaystyle \mathbf {S} (x)=\sum _{i}\mathbf {c} _{i}B_{i,p}(x).}$

Las B-splines de orden ${\displaystyle p+1}$ son funciones polinómicas unitarias de grado ${\displaystyle p}$ definidas sobre una cuadrícula de nudos ${\displaystyle {t_{0},\dots ,t_{i},\dots ,t_{m}}}$ (se utilizan índices basados en cero en adelante). El algoritmo de De Boor utiliza operaciones O(p²) + O(p) para evaluar la curva de spline. Nota: El artículo principal sobre B-splines y las publicaciones clásicas^[1]​ utilizan una notación diferente: la B-spline es indexada como ${\displaystyle B_{i,n}(x)}$${\displaystyle n=p+1}$.

Las B-splines tienen soporte local, lo que significa que los polinomios son positivos solamente en un ámbito finito y cero en otros lugares. La fórmula de recursión Cox-De Boor^[4]​ muestra esto::

${\displaystyle B_{i,0}(x):=\left\{{\begin{matrix}1&\mathrm {si} \quad t_{i}\leq x<t_{i+1}\\0&\mathrm {de\quad otra\quad forma} \end{matrix}}\right.}$

${\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).}$

Permitiendo que el índice ${\displaystyle k}$ defina el intervalo de nudos que contiene la posición, ${\displaystyle x\in [t_{k},t_{k+1}]}$. Podemos ver en la fórmula de recursión que solo B-splines con ${\displaystyle i=k-p,\dots ,k}$ no son ceros para este intervalo de nudos. Así, la suma se reduce a:

${\displaystyle \mathbf {S} (x)=\sum _{i=k-p}^{k}\mathbf {c} _{i}B_{i,p}(x).}$

De ${\displaystyle i\geq 0}$ se deduce que ${\displaystyle k\geq p}$. Del mismo modo, vemos en la recursividad que la ubicación del nudo más alto está en el índice ${\displaystyle k+1+p}$ . Esto significa que cualquier intervalo de nudos ${\displaystyle [t_{k},t_{k+1})}$ que se use realmente debe tener al menos ${\displaystyle p}$ nudos adicionales antes y después. En un programa de computadora, esto generalmente se logra repitiendo la primera y la última ubicación de nudo utilizada ${\displaystyle p}$ veces. Por ejemplo, para ${\displaystyle p=3}$ y ubicaciones de nudos reales ${\displaystyle (0,1,2)}$, uno podría rellenar el vector de nudos como ${\displaystyle (0,0,0,0,1,2,2,2,2)}$.

Con estas definiciones, ahora podemos describir el algoritmo de De Boor . El algoritmo no calcula las funciones de las B-splines ${\displaystyle B_{i,p}(x)}$ directamente. En cambio evalúa ${\displaystyle \mathbf {S} (x)}$ a través de una fórmula de recursión equivalente.

Deja a ${\displaystyle \mathbf {d} _{i,r}}$ ser un nuevo punto de control con ${\displaystyle \mathbf {d} _{i,0}:=\mathbf {c} _{i}}$ para ${\displaystyle i=k-p,\dots ,k}$. Para ${\displaystyle r=1,\dots ,p}$ la siguiente recursión es aplicada:

${\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}$

${\displaystyle \alpha _{i,r}={\frac {x-t_{i}}{t_{i+1+p-r}-t_{i}}}.}$

Una vez las iteraciones están completas, tenemos ${\displaystyle \mathbf {S} (x)=\mathbf {d} _{k,p}}$, esto significa que ${\displaystyle \mathbf {d} _{k,p}}$ es el resultado deseado.

El algoritmo De Boor es más eficaz que un cálculo explícito de B-splines ${\displaystyle B_{i,p}(x)}$ con la fórmula de recursión Cox-De Boor, porque no calcula términos que están garantizados para ser multiplicados por cero.

El algoritmo anterior no está optimizado para la implementación en un ordenador. Requiere memoria para ${\displaystyle (p+1)+p+\dots +1=(p+1)(p+2)/2}$ puntos de control provisional ${\displaystyle \mathbf {d} _{i,r}}$. Cada punto de control provisional es escrito exactamente una vez y leídos dos veces. Al invertir la iteración sobre ${\displaystyle i}$ (contando hacia abajo, en vez de hacia arriba), podemos ejecutar el algoritmo con memoria para solo ${\displaystyle p+1}$ puntos de control provisionales, permitiendo a ${\displaystyle \mathbf {d} _{i,r}}$ reutilizar la memoria para ${\displaystyle \mathbf {d} _{i,r-1}}$. De modo parecido, hay solo un valor de ${\displaystyle \alpha }$ utilizado en cada paso, de manera que podemos reutilizar la memoria también.

Además, es más conveniente utilizar un índice basado en cero ${\displaystyle j=0,\dots ,p}$ para los puntos de control provisionales. La relación al índice anterior es ${\displaystyle i=j+k-p}$. Por ello obtenemos el algoritmo mejorado:

Sea ${\displaystyle \mathbf {d} _{j}:=\mathbf {c} _{j+k-p}}$ para ${\displaystyle j=0,\dots ,p}$. Iterar para ${\displaystyle r=1,\dots ,p}$:

${\displaystyle \mathbf {d} _{j}:=(1-\alpha _{j})\mathbf {d} _{j-1}+\alpha _{j}\mathbf {d} _{j};\quad j=p,\dots ,r\quad }$(${\displaystyle j}$ debe ser contado hacia abajo)

${\displaystyle \alpha _{j}:={\frac {x-t_{j+k-p}}{t_{j+1+k-r}-t_{j+k-p}}}.}$

Después que las iteraciones se completan, el resultado es ${\displaystyle \mathbf {S} (x)=\mathbf {d} _{p}}$.

El código siguiente en el lenguaje de programación de Python es una implementación inocente ("naive") del algoritmo optimizado.

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]

• Algoritmo de De Casteljau
• Curva de Bézier
• NURBS

• De Boor's Algorithm

• PPPACK: Contiene muchos spline algoritmos en Fortran
• GNU Biblioteca científica: C-biblioteca, contiene un sub-biblioteca para splines ported de PPPACK
• SciPy: Pitón-biblioteca, contiene un sub-biblioteca scipy.interpolate Con spline las funciones basaron en FITPACK
• TinySpline: C-biblioteca para splines con un C++ wrapper y encuadernaciones para C#, Java, Lua, PHP, Pitón, y Ruby
• Einspline: C-biblioteca para splines en 1, 2, y 3 dimensiones con Fortran wrappers

• Carl de Boor (2003). A Practical Guide to Splines, Revised Edition. Springer-Verlag. ISBN 0-387-95366-3.