Кривые B-spline: Вычисление Коэффициентов

Хотя алгоритм de Boor является стандартным способом нахождения точки на кривой B-spline, соответствующей данному u, во многих случаях (напр., интерполяция и аппроксимация кривых) нам действительно нужны коэффициенты. Покажем простой способ сделать это.

Дана фиксированая кривая B-spline степени p, построенная по n +1 контрольным точкам p ₀, p ₁,..., p _n, и m +1 узлам u ₀= u ₁=...= u_p =0, u_{p +1}, u_{p +2},..., u_{m-p -1}, u_m-p = u_{m-p +1}=...= u_m =1, нам нужно вычислить коэффициенты N _{0, p} (u), N _{1, p} (u),..., N_n,p (u) для любого данного u на [0,1]. Простейший способ - это использовать следующее рекурсивное выражение:

Тем не менее, это очень затратный по времени процесс. Чтобы найти N_i,p (u), нам нужно вычислить N_{i,p -1}(u) и N_{i +1, p -1}(u). Чтобы найти N_{i,p -1}(u), нам нужно вычислить N_{i,p -2}(u) и N_{i +1, p -2}(u). Чтобы найти N_{i +1, p -1}(u), нам нужно N_{i +1, p -2}(u) и N_{i +2, p -2}(u). Как видите N_{i +1, p -2}(u) участвует дважды, и, в результате, ее рекурсивные вычисления также повторятся. Когда рекурсия идет вглубь, появляется все больше повторяющихся вычислений. Это очень похоже на рекурсивную версию алгоритма de Casteljau, который обсуждался ранее. Следовательно, скорость вычисления очень мала.

Есть простой способ сделать то же самое. Допустим, u лежит на узловом интервале [ u_k, u_{k +1}). Важное свойство базисных функций B-spline говорит о том, что базисные функции самое большее p +1 степени не равны нулю на [ u_k, u_{k +1}), а именно: N_k-p,p (u), N_{k-p +1, p} (u), N_{k-p +2, p} (u),..., N_{k -1, p} (u), N_{k , p} (u). По определению, единственная ненулевая базисная функция 0 степени на [ u_k, u_{k +1}) - это N_{k ,0}(u). В итоге, коэффициенты можно найти из N_{k ,0}(u) в треугольной форме в виде веера:

Так как N_{k ,0}(u) = 1 на узловом интервале [ u_k, u_{k +1}), а другие базисные функции 0 степени равны нулю на [ u_k, u_{k +1}), мы можем начать с N_{k ,0}(u) и вычислить базисные функции 1 степени N_{k -1,1}(u) и N_{k ,1}(u). Из этих двух значений 1 степени можно найти значения базисных функций 2 степени N_{k -2,2}(u), N_{k -1,2}(u) и N_{k ,2}(u). Этот процесс повторяется, пока не будут найдены все p +1 ненулевых коэффициента.

В этом вычислении, "внутренние" значения, такие как N_{k -1,2}(u), имеют как соседа справа сверху (т.e. N_{k -1,1}(u)), так и соседа справа снизу (т.e. N_{k ,1}(u)); значения на границе, идущей по направлению вправо вниз, такие как N_{k -1,1}(u) имеют только соседа слева снизу (т.e., N_{k ,0}(u)); а значения на границе, идущей по направлению влево вниз, такие как N_{k ,2}(u) имеют только соседа слева сверху (т.e., N_{k ,1}(u)). Таким образом, значения, находящиеся на границе, идущей по направлению вправо вверх (соотв., вправо вниз), исползуют второй (соотв. первый) член рекурсивного уравнения в определении. Только внутренние значения зависят от обоих членов. Основываясь на этом наблюдении, имеем следующий алгоритм:

Вход: n, p, m, u, и m +1 фиксированных узлов { u ₀,..., u_m }
Выход: Коэффициенты N _{0, p} (u), N _{1, p} (u),..., N_n,p (u) в [in] N[0], N[1],..., N[ n ]
Алгоритм:

Задаем N[0.. n ] равными 0; // инициализация
If u = u ₀ then // проверяем особые случаи

N[0] = 1.0;
return

else u = u_m then

N[ n ] = 1.0
return

End if

// теперь u между u ₀ и u_m

Пусть u лежит на узловом интервале [ u_k, u_{k +1});
R[ k ]:= 1.0 // коэффициент 0 степени
for d:=1 to p do // изменяем степень d от 1 до p

Begin

N[ k-d ] = * N[(k-d)+1]; // только правый член (юго-западный угол)
for i:= k-d +1 to k -1 do // вычисляем значения внутри

N[ i ]:= * N[ i ] + * N[ i +1];

N[ k ] = * N[ k ]; // только левый член (северо-западный угол)

End

// в массиве N[0.. n ] нужные коэффициенты.

Вышеизложенный алгоритм не очень эффективен. Его цель - просто и понятно проиллюстрировать смысл вычислений. Массив N[ ] содержит все промежуточные и конечные результаты. Для степени d элемент N[ i ] хранит значение N_i,d (u), а в конце, N[ k-d ], N[ k-d +1],..., N[ k ] содержат все ненулевые коэффициенты. Вычисление начинается с d =1, потому что мы знаем, что единственная ненулевая базисная функция - это N_{k ,0}(u), если u на узловом интервале [ u_k, u_{k +1}). Внешний цикл изменяет степень d от 1 до p. Первое выражение после begin вычисляет N_{k-d , d} (u), используя только один член выражения (т.e., его юго-западный сосед в треугольной схеме, N_{k-d +1, d -1}(u)), внцтренний цикл for вычисляет "внутренние" члены, а последнее выражение во внешнем цикле вычисляет N_k,d (u), используя только один член выражения (т.e., его северо-западный сосед в треугольной схеме, N_{k,d -1}(u)).

Если вы хорошо разберетесь в этом аалгоритме, сможете ли вы сделать его более эффективным?

Дата добавления: 2015-10-29