Алгоритм де Бура

Алгоритм де Бура — это численный метод вычисления значения B-сплайна в заданной точке; является обобщением алгоритма де Кастельжо для кривых Безье. Оригинальная версия алгоритма, разработанная Карлом де Буром в 1971 году, имеет полиномиальную вычислительную сложность, но отличается хорошей численной устойчивостью. Последующие модификации, появившиеся в попытке упростить и ускорить алгоритм, демонстрируют сравнительно худшую устойчивость.

Введение

Основная статья: B-сплайн

Алгоритм де Бура предназначен для вычисления сплайна $\mathbf {S} (x)$ при фиксированном значении аргумента $x$. Сам сплайн представляется в виде линейной комбинации контрольных точек $\mathbf {c} _{i}$ (векторов, не зависящих от $x$), коэффициенты которой суть базисные сплайн-функции $B_{i,p}(x)$:

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

Базисные сплайны по определению являются кусочно-полиномиальными функциями: они составлены из сегментов — многочленов степени $p$, — которые стыкуются между собой в узлах ${t_{0},\dots ,t_{i},\dots ,t_{m}}$. Алгоритм де Бура вычисляет $\mathbf {S} (x)$ за $O(p^{2})$ элементарных операций.

Примечание. В основной статье и в классических публикациях на эту тему приняты немного другие обозначения: B-сплайны индексируются как $B_{i,n}(x)$ с $n=p+1$. В данной статье второй индекс совпадает со степенью полиномов, а все индексы начинаются с нуля.

Ограниченность носителя

Каждый базисный сплайн имеет компактный (ограниченный) носитель, что означает, что он принимает ненулевое значения только на конечном отрезке. Это видно из рекурсивной формулы Кокса—де Бура:

$B_{i,0}(x):={\begin{cases}1&{\text{при }}t_{i}\leq x<t_{i+1},\\0&{\text{в противном случае.}}\end{cases}}$

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

Пусть индекс $k$ обозначает сегмент, в который попадает интересующее нас значение аргумента: $x\in [t_{k},t_{k+1})$. Тогда, согласно рекурсивной формуле, ненулевое значение могут иметь только B-сплайны с индексами $i=k-p,\dots ,k$. Таким образом, их линейную комбинацию можно упростить, исключив нулевые слагаемые:

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

Заметим, что индекс сумирования принимает осмысленные значения ($i\geq 0$) только при $k\geq p$. С другой стороны, рекурсивная формула для $k$-го сегмента включает узлы с индексами до $k+1+p$. Следовательно, чтобы находить значение сплайна для любого $x$ из отрезка $[t_{0},t_{m})$, необходимо иметь не менее $p$ дополнительных узлов до и после него. На практике этого обычно достигают, дублируя первый и последний узел нужное число раз. Например, если $p=3$ и даны узлы в точках $(0,1,2)$, то расширенный массив узлов будет иметь вид:$(0,0,0,0,1,2,2,2,2)$.

Алгоритм

Теперь мы готовы непосредственно описать алгоритм де Бура. Он не вычисляет B-сплайны $B_{i,p}(x)$ напрямую, а находит значение $\mathbf {S} (x)$ по эквивалентной рекурсивной формуле.

Введём новый набор контрольных точек $\mathbf {d} _{i,r}$; положим $\mathbf {d} _{i,0}:=\mathbf {c} _{i}$ для $i=k-p,\dots ,k$, а остальные $\mathbf {d} _{i,r}$найдём, последовательно применяя для $r=1,\dots ,p$ слудующую формулу:

$\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);$

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

После всех итераций, мы получим искомое значение в виде $\mathbf {S} (x)=\mathbf {d} _{k,p}$.

Алгоритм де Бура эффективнее, чем вычисление B-сплайнов $B_{i,p}(x)$ по формуле Кокса—де Бура, поскольку он заранее отбрасывает те слагаемые, которые гарантированно будут равны нулю.

Оптимизации

Приведённые выше алгоритм не оптимизирован для практической реализации на компьютере. Он предполагает хранение$(p+1)+p+\dots +1=(p+1)(p+2)/2$ промежеточных контрольных точек $\mathbf {d} _{i,r}$, каждая из которых считывается из памяти только два раза. Если пробегать по $i$ в обратном порядке, от большего значения к меньшему, можно обойтись $p+1$ изменяемой переменной для промежуточных точек, помещая $\mathbf {d} _{i,r}$ в ту переменную, которую занимала $\mathbf {d} _{i,r-1}$. Кроме того, поскольку на каждом шаге требуется единственное значение $\alpha$, его тоже можно хранить в одной и той же переменной.

Далее, вместо индекса $i$ удобнее использовать индекс $j$, который начинается с нуля: $j=i+p-k$. В результате получается следующий усовершенствованный алгоритм:

1. Устанавливаем $\mathbf {d} _{j}:=\mathbf {c} _{j+k-p}$ для $j=0,\dots ,p$.
2. Повторяем для всех $r=1,\dots ,p$:
   $\mathbf {d} _{j}:=(1-\alpha )\mathbf {d} _{j-1}+\alpha \mathbf {d} _{j}\quad (j=p,p-1,\dots ,r)$, где $\alpha ={\frac {x-t_{j+k-p}}{t_{j+1+k-r}-t_{j+k-p}}}$.
3. Возвращаем ответ: $\mathbf {S} (x)=\mathbf {d} _{p}$.

Пример реализации

Ниже приводится код усовершенствованного алгоритма на языке Python.

def deBoor(k: int, x: int, t, c, p: int):
    """Вычисление сплайна S(x) по алгоритму де Бура.

    Входные параметры
    -----------------
    k: Индекс сегмента, содержащего x.
    x: Точка, в которой требуется вычислить сплайн.
    t: Массив узлов (расширенный дополнительными узлами, см. статью).
    c: Массив контрольных точек.
    p: Степень сплайна.
    """
    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]

См. также

• Алгоритм де Кастельжо
• Кривая Безье
• NURBS

Внешние ссылки

• De Boor's Algorithm
• The DeBoor-Cox Calculation

Промышленные реализации

Алгоритм де Бура реализован во многих популярных библиотеках на разных языках программирования:

• PPPACK: содержит алгоритмы для работы со сплайнами на языке Fortran.
• GNU Scientific Library: библиотека на языке C, включает в себя портированную версию предыдущей библиотеки.
• SciPy: библиотека для Python; функции для работы со спалайнами содержатся в модуле scipy.interpolate, они основываются на FITPACK.
• TinySpline: специализированная библиотека на C, имеет также интерфейсы для C++, C#, Java, Lua, PHP, Python и Ruby.
• Einspline: ещё реализация на C с обёрткой для Fortran; поддерживает 1-, 2- и 3-мерные сплайны.

Примечания

Литература

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