В этом примере показано, как создавать и работать с B-формой сплайна в 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 может использоваться для получения любой или всех частей В-формы. Например, вот выходные данные, предоставляемые 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
Сплайн в панели инструментов фитинга кривой также может быть многомерным, а именно тензорным произведением одномерных сплайнов. B-форма для такой функции является лишь немного более сложной, так как узлы теперь представляют собой матрицу ячеек, содержащую различные одномерные последовательности узлов, и матрицу коэффициентов, соответственно, многомерную.
Например, эта случайная сплайновая поверхность кубическая в первой переменной (в этой переменной есть 11 узла и коэффициенты 7), но только кусочно постоянная во второй переменной ((2 + 5 + 2) -8 = 1).
fnplt( spmak({augknt(0:4,4),augknt(0:4,3)}, rand(7,8)) )