В этом примере показано, как создать и работать с B-формой сплайна в Curve Fitting Toolbox™.
В Curve Fitting Toolbox, кусочном полиноме или стр, функция в B-форме часто вызывается сплайн.
B-форма (одномерные) стр задана ее (не уменьшающейся) последовательностью узла t и его содействующей последовательностью B-сплайна a.
Учитывая последовательность узла и содействующую последовательность стр, команда 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-сплайн th порядка 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}) стр порядка k с пропусками в t_j, ..., t_{j+k} (и больше нигде). На самом деле его нетривиальные полиномиальные части являются всей точной степенью k-1.
Например, самый правый B-сплайн выше включает 5 узлов, следовательно имеет порядок 4, i.e., кубический 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 условия гладкости, i.e., только непрерывность через тот узел функции и ее первой производной.
Любая B-форма может быть усовершенствована, i.e., преобразованный, вставкой узла, в 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)) )