В этом примере показано, как создать и работать с B-формой сплайна в Curve Fitting Toolbox™.
В Curve Fitting Toolbox кусочного полинома, или pp, функция в B-форме часто называется сплайном.
B-форма (одномерного) pp задается его (неразрешающей) последовательностью узлов t и своей последовательностью B-сплайн коэффициентов a.
Учитывая последовательность узлов и последовательность коэффициентов pp, команда spmak возвращает соответствующую B-форму для использования в таких командах, как fnval (вычислите функцию), fnplt (постройте график функции), fnder (дифференцирование функции) и другие связанные команды.
Получившийся сплайн имеет порядок k := length(t) - size(a,2). Это означает, что все его полиномиальные части имеют степень < k.
Сказать, что сплайн s имеет узлы t и коэффициенты a означает, что
n
s(x) := sum B_{j,k}(x) * a(:,j),
j=1
где B_(j,k) = B( . | t_j, ..., t_{j+k}) является j-й B-сплайн порядка k для заданной последовательности узлов t, т.е. B-сплайн с узлами t_j, ..., t_{j+k}. Для примера,
t = [.1 .4 .5 .8 .9]; a = 1; fnplt(spmak(t,a),2.5); tmp = repmat(t,3,1); ty = repmat(.1*[1;-1;NaN],1,5); hold on plot(tmp(:),ty(:)) text(.65,.5,'B( \cdot | .1, .4, .5, .8, .9)','FontSize',12) text(.05,1.,'s(x) = \Sigma_j B( x | t_j , \ldots, t_{j+k} ) a(:,j)', ... 'FontSize',16,'Color','r') axis([0 1 -.2 1.2]) title('A B-spline of Order 4') hold off
Значение сплайна
s(x) = sum B_{j,k}(x) a(:,j)
j
в любое x в узловом интервале [t_i .. t_{i+1}] является выпуклой комбинацией k коэффициенты a(:.i-k+1), ..., a(:,i) так как на этом интервале только k B-сплайны B_{i-k+1,k}, ..., B_{i,k} ненулевые, и они неотрицательны там и сумма до 1, как показано на следующем рисунке.
Это часто суммируется, говоря, что B-сплайны обеспечивают локальное (неотрицательное) разбиение единицы.
k = 3; n = 3; t = [1 1.7 3.2 4.2 4.8 6]; tt = (10:60)/10; vals = fnval(spmak(t,eye(k)),tt); plot(tt.',vals.'); hold on ind = find(tt>=t(3)&tt<=t(4)); plot(tt(ind).',vals(:,ind).','LineWidth',3) plot(t([3 4]),[1 1],'k','LineWidth',3) ty = repmat(.1*[1;-1;NaN],1,6); plot([0 0 -.2 0 0 -.2 0 0],[-.5 0 0 0 1 1 1 1.5],'k') text(-.5,0,'0','FontSize',12) text(-.5,1,'1','FontSize',12) tmp = repmat(t,3,1); plot(tmp(:),ty(:),'k'); yd = -.25; text(t(1),yd,'t_{i-2}','FontSize',12); text(t(3),yd,'t_i','FontSize',12); text(t(4),yd,'t_{i+1}','FontSize',12); text(1.8,.5,'B_{i-2,3}','FontSize',12); text(5,.45,'B_{i,3}','FontSize',12); axis([-.5 7 -.5 1.5]) title('B-splines Form a Local Partition of Unity') axis off hold off
Когда коэффициенты являются точками в плоскости и, соответственно, сплайном
s(x) = sum B_{j,k}(x) a(:,j)
j
трассирует кривую, это означает, что часть кривой
{ s(x) : t_i <= x <= t_{i+1} }
выделенный на рисунке ниже большей LineWidth, лежит в выпуклой оболочке, показанной желтым цветом на рисунке ниже, k точки a(:,i-k+1), ... a(:,i).
t = 1:9; c = [2 1.4;1 .5; 2 -.4; 5 1.4; 6 .5; 5 -.4].'; sp = spmak(t,c); fill(c(1,3:5),c(2,3:5),'y','EdgeColor','y'); hold on fnplt(sp,t([3 7]),1.5) fnplt(sp,t([5 6]),3) plot(c(1,:),c(2,:),':ok') text(2,-.55,'a(:,i-2)','FontSize',12) text(5,1.6,'a(:,i-1)','FontSize',12) text(6.1,.5,'a(:,i)','FontSize',12) title('The Convex-Hull Property') axis([.5 7 -.8 1.8]) axis off hold off
Для квадратичного сплайна (т.е. k = 3), как показано здесь, это даже означает, что кривая является касательной к управляющему многоугольнику (показан как пунктирные линии). Это ломаная линия, которая соединяет коэффициенты, которые в этой связи называются управляющими точками (показаны здесь как открытые круги).
Мы можем думать о графике скалярного сплайна
s = sum B_{j,k}*a(j)
j
как кривая x |--> (x,s(x)). С тех пор
x = sum B_{j,k}(x) t^*_j
j
для x в интервале [t_k .. t_{n+1}], где
t^*_i := (t_{i+1} + ... + t_{i+k-1})/(k-1) for i = 1:n
являются ли средние значения узлов, получаемые с помощью aveknt команда, управляющий многоугольник для скалярного сплайна является ломаной линией с вершинами (t^*_i, a(i)), i=1:n.
В примере ниже показан кубический сплайн (k = 4), с 4-кратными концевыми узлами, отсюда
t^*_1 = t_1 and t^*_n = t_{n+k}.
t = [0 .2 .35 .47 .61 .84 1]*(2*pi); s = t([1 3 4 5 7]); knots = augknt(s,4); sp = spapi(knots,t,sin(t)+1.8); fnplt(sp,2); hold on c = fnbrk(sp,'c'); ts = aveknt(knots,4); plot(ts,c,':ok'); tt = [s;s;NaN(size(s))]; ty = repmat(.25*[-1;1;NaN], size(s)); plot(tt(:),ty(:),'r') plot(ts(1,:),zeros(size(ts)),'*') text(knots(5),-.5,'t_5','FontSize',12) text(ts(2),-.45,'t^*_2','FontSize',12) text(knots(1)-.28,-.5,'t_1=t_4','FontSize',12) text(knots(end)-.65,-.45,'t_{n+1}=t^*_n=t_{n+4}','FontSize',12) title('A Cubic Spline and its Control Polygon') axis([-.72 7 -.5 3.5]) axis off hold off
savesp = sp;
Существенными частями B-формы являются последовательность узлов t и последовательность коэффициентов B-сплайна a. Другие части являются числом n из связанных B-сплайнов или коэффициентов порядок k его полиномиальных частей и размерности d его коэффициентов a. В частности, size(a) равен [d,n].
Существует еще одна часть, а именно основной интервал, [t(1) .. t(end)]. Он используется как интервал по умолчанию при построении графика функции. Кроме того, сплайн принимается непрерывным справа везде, кроме правой конечной точки основного интервала, где он принимается непрерывным слева. Это проиллюстрировано в примере ниже, для сплайна, созданного с помощью spmak.
b = 0:3; sp = spmak(augknt(b,3),[-1,0,1,0,-1]); x = linspace(-1,4,51); plot(x,fnval(sp,x),'x') hold on axis([-2 5,-1.5,1]) tx = repmat(b,3,1); ty = repmat(.5*[1;-1;NaN],1,length(b)); plot(tx(:),ty(:),'-r') legend({'Spline Values' 'Knots'}) hold off title({'A Spline in B-form is Right(Left)-Continuous ';... 'at Left(Right) Endpoint of its Basic Interval'})
fnbrk может использоваться для получения любой или всех частей B-формы. Для примера вот выход, предоставляемая fnbrk для B-формы сплайна, показанного выше.
fnbrk(sp)
The input describes a B- form
knots(1:n+k)
0 0 0 1 2 3 3 3
coefficients(d,n)
-1 0 1 0 -1
number n of coefficients
5
order k
3
dimension d of target
1
Однако обычно нет необходимости знать какие-либо из этих частей. Скорее вы используете такие команды, как spapi или spaps чтобы создать B-форму сплайна из некоторых данных, используйте такие команды, как fnval, fnplt, fnder, и т.д., чтобы использовать сплайн, построенный, без необходимости смотреть на его различные части.
В следующих разделах приводится более подробная информация о B-сплайнах, в частности о важной роли, которую играют узловые кратности.
Здесь, для k = 2, 3 и 4 являются B-сплайнами порядка kи ниже них их первая и вторая (кусочно) производные, чтобы проиллюстрировать некоторые факты о B-сплайнах. Опробуйте bpspligui инструмент, если вы хотите экспериментировать с собственными примерами.
cl = ['g','r','b','k','k']; v = 5.4; d1 = 2.5; d2 = 0; s1 = 1; s2 = .5; t1 = [0 .8 2]; t2 = [3 4.4 5 6]; t3 = [7 7.9 9.2 10 11]; tt = [t1 t2 t3]; ext = tt([1 end])+[-.5 .5]; plot(ext([1 2]),[v v],cl(5)) hold on plot(ext([1 2]),[d1 d1],cl(5)) plot(ext([1 2]),[d2 d2],cl(5)) ts = [tt;tt;NaN(size(tt))]; ty = repmat(.2*[-1;0;NaN],size(tt)); plot(ts(:),ty(:)+v,cl(5)) plot(ts(:),ty(:)+d1,cl(5)) plot(ts(:),ty(:)+d2,cl(5)) b1 = spmak(t1,1); p1 = [t1;0 1 0]; db1 = fnder(b1); p11 = fnplt(db1,'j'); p12 = fnplt(fnder(db1)); lw = 2; plot(p1(1,:),p1(2,:)+v,cl(2),'LineWidth',lw) plot(p11(1,:),s1*p11(2,:)+d1,cl(2),'LineWidth',lw) plot(p12(1,:),s2*p12(2,:)+d2,cl(2),'LineWidth',lw) b1 = spmak(t2,1); p1 = fnplt(b1); db1 = fnder(b1); p11 = [t2;fnval(db1,t2)]; p12 = fnplt(fnder(db1),'j'); plot(p1(1,:),p1(2,:)+v,cl(3),'LineWidth',lw) plot(p11(1,:),s1*p11(2,:)+d1,cl(3),'LineWidth',lw) plot(p12(1,:),s2*p12(2,:)+d2,cl(3),'LineWidth',lw) b1 = spmak(t3,1); p1 = fnplt(b1); db1 = fnder(b1); p11 = fnplt(db1); p12=[t3;fnval(fnder(db1),t3)]; plot(p1(1,:),p1(2,:)+v,cl(4),'LineWidth',lw) plot(p11(1,:),s1*p11(2,:)+d1,cl(4),'LineWidth',lw) plot(p12(1,:),s2*p12(2,:)+d2,cl(4),'LineWidth',lw) tey = v+1.5; text(t1(2)-.5,tey,'linear','FontSize',12,'Color',cl(2)) text(t2(2)-.8,tey,'quadratic','FontSize',12,'Color',cl(3)) text(t3(3)-.5,tey,'cubic','FontSize',12,'Color',cl(4)) text(-2,v,'B','FontSize',12) text(-2,d1,'DB','FontSize',12) text(-2,d2,'D^2B') axis([-1 12 -2 7.5]) title({'B-splines with Simple Knots and Their Derivatives'}) axis off hold off
1. B-сплайн B_{j,k} = B( . | t_j, ..., t_{j+k}) является pp порядка k с пропусками в t_j, ..., t_{j+k} (и нигде больше). На самом деле, его нетривиальные полиномиальные части являются точной степенью k-1.
Для примера самый правый B-сплайн выше включает 5 узлов, поэтому имеет порядок 4, т.е. кубический B-сплайн. Соответственно, его вторая производная является кусочно-линейной.
2. B_{j,k} положительно на интервале (t_j .. t_{j+k}) и нуль отключает этот интервал. Он также исчезает в конечных точках этого интервала, если только конечная точка не является узлом кратности k (см. самый правый пример на следующем рисунке).
3. Кратность узлов определяет гладкость, с которой два соседних полинома соединяются поперек этого узла. В сокращённом виде правило является:
knot multiplicity + number of smoothness conditions = order
Чтобы проиллюстрировать эту последнюю точку, на рисунке ниже показаны четыре кубических B-сплайна и, ниже них, их первые две производные. Каждый сплайн имеет определенный узел кратности 1, 2, 3, 4, что обозначено длинами узловых линий.
d2 = -1; t1 = [7 7.9 9.2 10 11]-7; t2 = [7 7.9 7.9 9 10]-2; t3 = [7 7.9 7.9 7.9 9]+2; t4 = [7 7.9 7.9 7.9 7.9]+5; [m,tt] = knt2mlt([t1 t2 t3 t4]); ext = tt([1 end])+[-.5 .5]; plot(ext,[v v],cl(5)) hold on plot(ext,[d1 d1],cl(5)) plot(ext,[d2 d2],cl(5)) ts = [tt;tt;NaN(size(tt))]; ty = .2*[-m-1;zeros(size(m));m]; plot(ts(:),ty(:)+v,cl(5)) plot(ts(:),ty(:)+d1,cl(5)) plot(ts(:),ty(:)+d2,cl(5)) b1 = spmak(t1,1); p1 = fnplt(b1); db1 = fnder(b1); p11 = fnplt(db1); p12 = [t1;fnval(fnder(db1),t1)]; plot(p1(1,:),p1(2,:)+v,cl(1),'LineWidth',lw) plot(p11(1,:),s1*p11(2,:)+d1,cl(1),'LineWidth',lw) plot(p12(1,:),s2*p12(2,:)+d2,cl(1),'LineWidth',lw) text(-2,v,'B'), text(-2,d1,'DB'), text(-2,d2,'D^2B') b1 = spmak(t2,1); p1 = fnplt(b1); db1 = fnder(b1); p11 = fnplt(db1); p12 = fnplt(fnder(db1),'j'); plot(p1(1,:),p1(2,:)+v,cl(2),'LineWidth',lw) plot(p11(1,:),s1*p11(2,:)+d1,cl(2),'LineWidth',lw) plot(p12(1,:),s2*s2*p12(2,:)+d2,cl(2),'LineWidth',lw) b1 = spmak(t3,1); p1 = fnplt(b1); db1 = fnder(b1); p11 = fnplt(db1,'j'); p12 = fnplt(fnder(db1),'j'); plot(p1(1,:),p1(2,:)+v,cl(3),'LineWidth',lw) plot(p11(1,:),s1*s2*p11(2,:)+d1,cl(3),'LineWidth',lw) plot(p12(1,:),s2*s2*p12(2,:)+d2,cl(3),'LineWidth',lw) b1 = spmak(t4,1); p1 = fnplt(b1); db1 = fnder(b1); p11 = fnplt(db1); p12 = fnplt(fnder(db1)); plot(p1(1,:),p1(2,:)+v,cl(4),'LineWidth',lw) plot(p11(1,:),s2*p11(2,:)+d1,cl(4),'LineWidth',lw) plot(p12(1,:),s2*s2*p12(2,:)+d2,cl(4),'LineWidth',lw) text(t2(2)-.5,tey,'2-fold','FontSize',12,'Color',cl(2)) text(t3(2)-.8,tey,'3-fold','FontSize',12,'Color',cl(3)) text(t4(3)-.8,tey,'4-fold','FontSize',12,'Color',cl(4)) axis([-1 14 -3 7.5]) title('Cubic B-splines With A Knot of Various Multiplicities') axis off hold off
Для примера, поскольку порядок кубического B-сплайна равен 4, двойной узел означает всего 2 условия гладкости, то есть просто непрерывность через этот узел функции и его первую производную.
Любая B-форма может быть уточнена, т.е. преобразована, посредством вставки узла, в B-форму для той же функции, но для более мелкой последовательности узлов. Чем ближе последовательность узлов, тем ближе многоугольник управления к представляемой функции.
Например, этот рисунок показывает исходные (черным) и уточненные (красным) многоугольники управления для кубического сплайна, используемого ранее.
sp = savesp; fnplt(sp,2.5); hold on c = fnbrk(sp,'c'); plot(aveknt(fnbrk(sp,'k'),4),c,':ok'); b = knt2brk(fnbrk(sp,'k')); spref = fnrfn(sp,(b(2:end)+b(1:end-1))/2); cr = fnbrk(spref,'c'); h2 = plot(aveknt(fnbrk(spref,'knots'),4),cr,':*r'); axis([-.72 7 -.5 3.5]) title('A Spline, its Control Polygon, and a Refined Control Polygon') axis off hold off
В качестве второго примера мы начинаем с вершин стандартного алмаза в качестве наших контрольных точек, но дважды пробегаем последовательность.
ozmz = [1 0 -1 0]; c = [ozmz ozmz 1; 0 ozmz ozmz]; circle = spmak(-4:8,c); fnplt(circle) hold on plot(c(1,:),c(2,:),':ok') axis(1.1*[-1 1 -1 1]) axis equal, axis off hold off
Однако, когда мы строим график получившегося сплайна, мы получаем кривую, которая начинается и заканчивается в источник, из-за того, что мы решили сделать последовательность узлов простой. Отсюда наш сплайн исчезает в конечных точках его основного интервала, [-4 .. 8]. Нам действительно нужна только та часть сплайна, которая соответствует интервалу [0 .. 4], более смело изображен на рисунке ниже.
fnplt(circle) hold on fnplt(circle,[0 4],4) plot(c(1,:),c(2,:),':ok') axis(1.1*[-1 1 -1 1]) title('A Circle as Part of a Spline Curve with a Diamond as Control Polygon') axis equal, axis off hold off
Чтобы получить только круг, мы ограничиваем наш сплайн интервалом [0 .. 4]. Мы делаем это, преобразовывая в ppform, ограничивая [0 .. 4], затем преобразование в B-форму.
circ = fn2fm(fnbrk(fn2fm(circle,'pp'),[0 4]),'B-'); fnplt(circ,2.5) hold on cc = fnbrk(circ,'c'); plot(cc(1,:),cc(2,:),':ok') axis(1.1*[-1 1 -1 1]) axis equal, axis off hold off
Уточнение получившейся последовательности узлов приводит к управляемому многоугольнику гораздо ближе к кругу.
ccc = fnbrk(fnrfn(circ,.5:4),'c'); hold on plot(ccc(1,:),ccc(2,:),'-r*') title('A Circle as a Spline Curve, its Control Polygon, and a Refinement') hold off
Сплайн в Curve Fitting Toolbox также может быть многомерным, а именно, тензорным продуктом одномерных сплайнов. B-форма для такой функции лишь немного сложнее, с узлами теперь массив ячеек, содержащий различные одномерные последовательности узлов, и массив коэффициентов, подходящим образом многомерный.
Для примера эта случайная поверхность сплайна кубическая по первой переменной (в этой переменной существует 11 узлов и 7 коэффициентов), но только кусочно-постоянная по второй переменной ((2 + 5 + 2) -8 = 1).
fnplt( spmak({augknt(0:4,4),augknt(0:4,3)}, rand(7,8)) )