Полтора дня разбирался в работе B-сплайна и как это применяется на практике. Хочу поделиться с вами результатами моих изысканий.
Немного теории
Допустим, у вас есть некий массив точек размером n. Что вам потребуется для работы B-сплайна:
- определить количество интерполируемых точек (от этого зависит скорость работы, гладкость получившейся кривой и размер формулы для вычисления, количество точек равно k+1)
- избавится от рекурсии формулы для вычисления
Общая формула для расчета коэффициентов:
N_{i,k}(x)=\frac{x-t_i}{t_{i+k-1}-t_i}N_{i,k-1}(x)+\frac{t_{i+k}-x}{t_{i+k}-t_{i+1}}N_{i+1,k-1}(x)
Формула рекурсивная и заканчивается на k равном единице:
N_{i,1}(x)=\begin{cases}1&\text{if }t_i\leq x\text{\textless} t_{i+1}\\0&\mathrm{otherwise}\quad\end{cases}
Пояснения к формуле:
- если делитель получается равным или близким по значению к нулю, значит вся дробь равна нулю
- параметр i — это индекс текущей точки
- x — значение от i до i + 1 с нужным вам шагом (например, i равно 7, и нам нужно создать 3 точки для интерполяции, значит берем x равным 7.0, 7.5, 8.0)
- t — неубывающий массив индексов
Вычисление элементов массива t:
t=\begin{cases}0,&\mathrm{if}\ i<k\\i-k+1,&\mathrm{if}\ k\leq{i}\leq{n}\\n-k+2,&\mathrm{if}\ i>n\end{cases}
Практика
Мне нужно было интерполировать по четырем точкам, поэтому я заранее на бумажке вычислил формулу без рекурсии:
N_{i,3}(x)=\begin{cases}0,&\mathrm{if}\ x\text{\textless} t_{i}\\\frac{x-t_i}{t_{i+2}-t_i}*\frac{x-ti}{t_{i+1}-ti},&\mathrm{if}\ t_i\leq x\text{\textless} t_{i+1}\\\frac{x-ti}{t_{i+2}-ti}*\frac{t_{i+2}-x}{t_{i+2}-t_{i+1}}+\frac{t_{i+3}-x}{t_{i+3}-t_{i+1}}*\frac{x-ti}{t_{i+2}-t_{i+1}},&\mathrm{if}\ t_{i+1}\leq x\text{\textless} t_{i+2}\\\frac{t_{i+3}-x}{t_{i+3}-t_{i+1}}*\frac{t_{i+3}-x}{t_{i+3}-t_{i+2}},&\mathrm{if}\ t_{i+2}\leq x\text{\textless} t_{i+3}\\0,&\mathrm{if}\ x\geq t_{i+3}\end{cases}Про бумажку я, конечно, соврал, мне было удобнее сделать функцию на питоне вот такого вида:
def div(a, b):
try:
return a / float(b)
except ZeroDivisionError:
return 0
def N(i, k, x, t):
if 1 == k:
if t[i] <= x < t[i + 1]:
return 1.0
return 0.0
a = div(x - t[i], t[i + k - 1] - t[i]) * N(i, k - 1, x, t)
b = div(t[i + k] - x, x[i + k] - x[i + 1]) * N(i + 1, k - 1, x, t)
return a + b
После этого добавил условие k==2 и все возможные возвращаемые значения в зависимости от x, затем условие k==3, значения которого зависят от k==2. Может вам будет удобнее все это сделать на бумажке.
Наконец можно приступать к интерполяции. Формула такая: \sum\limits_{i}^{i+k}=\vec{p_{i}}N_{i,k}(x)
Если индекс точки меньше нуля, берем нулевую точку, если индекс точки больше размера массива, то берем последнюю точку в массиве.
Небольшой пример
У нас есть некий массив точек, нам надо интерполировать с пятой по шестую точки с шагом 0.25, при этом $k$ равно 3 (значит для интерполяции нужно взять четыре точки).
p_0=p_5*N_{0,3}(0.0)+p_6*N_{1,3}(0.0)+p_7*N_{2,3}(0.0)+p_8*N_{3,3}(0.0)\newline p_1=p_5*N_{0,3}(0.25)+p_6*N_{1,3}(0.25)+p_7*N_{2,3}(0.25)+p_8*N_{3,3}(0.25)\newline p_2=p_5*N_{0,3}(0.5)+p_6*N_{1,3}(0.5)+p_7*N_{2,3}(0.5)+p_8*N_{3,3}(0.5)\newline p_3=p_5*N_{0,3}(1.0)+p_6*N_{1,3}(1.0)+p_7*N_{2,3}(1.0)+p_8*N_{3,3}(1.0)\newlineИтого четыре новых точки вместо двух.
Спасибо за внимание, надеюсь нигде не ошибся. На несложные вопросы по теме могу ответить в комментариях.