Введение

В 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)) )