Main Content

B 型の作成と操作

この例では、Curve Fitting Toolbox™ で B 型のスプラインを作成および操作する方法を示します。

はじめに

Curve Fitting Toolbox では、しばしば B 型の区分的多項式 (pp) 関数を "スプライン" と呼びます。

"B 型" の (一変量) pp は、非減少の "節点" シーケンス t と B スプライン "係数" シーケンス a によって指定されます。

pp の節点シーケンスと係数シーケンスを指定すると、コマンド spmakfnval (関数の評価)、fnplt (関数のプロット)、fnder (関数の微分) および他の関係するコマンドで使用する、対応する B 型を返します。

結果のスプラインは "次数" 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}) は、与えられた "節点シーケンス" t に対する "次数" kj 番目の "B スプライン" です。つまり、節点 t_j, ..., t_{j+k} を持つ B スプラインです。たとえば、

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

Figure contains an axes object. The axes object with title A B-spline of Order 4 contains 4 objects of type line, text.

1 のローカル分割

次のスプラインの値

s(x) = sum B_{j,k}(x) a(:,j)

j

は節点区間 [t_i .. t_{i+1}] の任意の xk 係数 a(:.i-k+1), ..., a(:,i) の "凸" 結合となります。この区間では、k の B スプライン B_{i-k+1,k}, ..., B_{i,k} のみが非ゼロであり、そこで非負かつ合計が 1 になるためです (次の Figure を参照)。

これはしばしば、B スプラインで "1 がローカル (非負) 分割" されるとして説明されます。

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} }

が次の Figure で太線にて強調表示されていますが、これはその Figure の黄色で示された ka(:,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

ここで示した 2 次スプライン (つまり、k = 3) の場合、この曲線は "制御多角形" に対する正接 (破線で表示) となります。この破線は係数を接続しており、この接続の "制御点" と呼ばれます (白丸で表示)。

スカラー値のスプラインの制御多角形

"スカラー" 値のスプラインのグラフ

s = sum B_{j,k}*a(j)

j

について曲線 x |--> (x,s(x)) として検討します。ここで

x = sum B_{j,k}(x) t^*_j

j

が区間 [t_k .. t_{n+1}]x について成り立つため (ここで

t^*_i := (t_{i+1} + ... + t_{i+k-1})/(k-1) for i = 1:n

aveknt コマンドを用いて取得可能な節点平均)、スカラー値のスプラインの "制御多角形" は頂点 (t^*_i, a(i)), i=1:n をもつ破線となります。

次の例に、両端の節点が 4 次の "3 次" スプライン (k = 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 です。その他の部分は、関係する B スプラインまたは係数の "数" n、その多項式の項の "次数" k およびその係数 a の "次元" d です。特に、size(a)[d,n] と等しくなります。

もう 1 つ、"基本区間" [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'})

Figure contains an axes object. The axes object with title A Spline in B-form is Right(Left)-Continuous at Left(Right) Endpoint of its Basic Interval contains 2 objects of type line. One or more of the lines displays its values using only markers These objects represent Spline Values, Knots.

fnbrk は、B 型のいずれかの部分またはすべての部分を取得するために使用できます。たとえば、上に示した B 型スプラインの fnbrk で与えられる出力を次に示します。

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 型のスプラインを作成してから、fnvalfnpltfnder などのコマンドを使用し、作成したスプラインを利用するので、各部分を確認する必要はありません。

以降のセクションでは、B スプラインに関する詳細、特に節点の "多重度" が果たす重要な役割についての詳細情報を示します。

節点の多重度

ここでは、k = 2、3 および 4 の場合について、次数 k の B スプラインを示し、次にその 1 階微分および 2 階 (区分的) 微分を示して 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}) は、t_j, ..., t_{j+k} (のみ) にブレークがある次数 k の pp です。実際、その自明でない多項式の項はすべて次数 k-1 となります。

たとえば、上記の右端の B スプラインには 5 つの節点が含まれているため、次数 4、すなわち "3 次の" B スプラインとなります。それに応じて、この 2 階微分は区分的線形となります。

2. B_{j,k} は区間 (t_j .. t_{j+k}) では正になり、この区間外では 0 になります。また、その区間の端点が多重度 k の節点でない場合は、端点でも 0 になります (次の Figure の右端の例を参照)。

3.節点の "多重度" により、2 つの隣接する多項式がその節点をまたがって連結する際の平滑度が決まります。簡単に言えば、規則は次のとおりです。

knot multiplicity + number of smoothness conditions = order

この最後の点を説明するため、次の Figure に 4 つの 3 次 B スプラインを示し、次にそれらの最初の 2 つの微分を示します。節点のラインの長さで示されているように、それぞれのスプラインの節点の多重度は 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

たとえば、3 次の B スプラインの次数は 4 であるため、二重節点は 2 つの平滑性条件のみ、つまり関数の節点をまたがる連続性とその 1 階微分のみを意味します。

調整および節点挿入

どの B 型も、"節点挿入" により、"調整"、つまり同じ関数であってもより細かい節点シーケンスの B 型への変換が可能です。節点シーケンスが細かくなればなるほど、制御多角形は表現される関数に近くなります。

たとえば、次の Figure は、前に使用した 3 次スプラインの元の制御多角形 (黒) と調整後の制御多角形 (赤) を示しています。

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

2 番目の例として、まず制御点に標準的な菱形の頂点を使用し、このシーケンスを 2 回処理します。

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 になります。本当に必要なのは、区間 [0 .. 4] に対応するスプラインの部分のみで、この区間は次の Figure に太字でプロットされています。

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] に制限します。これは、pp 型への変換、[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 型は、節点がさまざまな一変量節点シーケンスを含む cell 配列となり、係数配列が適切な多次元配列となるため、少しだけ複雑度が増します。

たとえば、次のランダム スプライン曲面は最初の変数が 3 次 (この変数内に節点が 11 個、係数が 7 個) ですが、2 番目の変数は区分的定数のみです ((2+5+2)-8 = 1)。

fnplt( spmak({augknt(0:4,4),augknt(0:4,3)}, rand(7,8)) )

Figure contains an axes object. The axes object contains an object of type surface.