Алгоритм де Бура (Glikjnmb ;y >rjg)

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

Введение

Алгоритм де Бура предназначен для вычисления сплайна $\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$ . В результате получается следующий усовершенствованный алгоритм:

Устанавливаем $\mathbf {d} _{j}:=\mathbf {c} _{j+k-p}$ для $j=0,\dots ,p$ .
Повторяем для всех $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}}}$ .
Возвращаем ответ: $\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]

См. также

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

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

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

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-мерные сплайны.

Примечания

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

Литература

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

[de_boor_paper-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. 29 (4). Springer-Verlag: 365—371. doi:10.1007/BF02246763. S2CID 2407104.

[3] Lee, E. T. Y. (1986). "Comments on some B-spline algorithms". Computing. 36 (3). Springer-Verlag: 229—238. doi:10.1007/BF02240069. S2CID 7003455.

[1]

[2]

[3]