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

Введение

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

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

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

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

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

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

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

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

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

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

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

Алгоритм

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

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

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

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

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

Оптимизации

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

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

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

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

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

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

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

Примечания

1. 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.
2. Lee, E. T. Y. (December 1982). "A Simplified B-Spline Computation Routine". Computing . Springer-Verlag. 29 (4): 365—371. doi : . S2CID .
3. Lee, E. T. Y. (1986). "Comments on some B-spline algorithms". Computing . Springer-Verlag. 36 (3): 229—238. doi : . S2CID .

Литература

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