平面のスプライン

スクリプトを開く

この例では、Curve Fitting Toolbox™ の spmak、spcrv、cscvn および rscvn の各コマンドを使用して、平面のスプライン曲線を作成する方法を示します。また、正接のプロットおよび曲線で囲まれた面積の計算も行います。

簡単なスプライン曲線

Curve Fitting Toolbox は "ベクトル値の" スプラインを処理できます。d ベクトル値の一変量スプラインは d 空間の曲線を提供します。このモードでは、平面曲線が得られるので d = 2 が最も一般的です。

次に、2 次元の係数をもつスプラインを作成し、プロットする例を示します。

knots = [1,1:9,9];
curve = spmak( knots, repmat([ 0 0; 1 0; 1 1; 0 1 ], 2,1).' );

t = linspace(2,8,121);
values = fnval(curve,t);
plot(values(1,:),values(2,:),'LineWidth',2);
axis([-.4 1.4 -.2 1.2]), axis equal
title('A Spline Curve');
hold on

注意

この例では、fnplt を使用して曲線をプロットせず、代わりに fnval で得られた曲線に対して点をプロットしたことがわかります。次にそのコードを再度示します。

  t = linspace(2,8,121);
  values = fnval(curve,t);
  plot(values(1,:),values(2,:),'LineWidth',2)

この固有スプライン曲線に fnplt を直接使用すると、次の Figure の赤色の曲線が得られます。

fnplt(curve,'r',.5);
title('The Full Spline Curve, in Red')

どうしてこうなるのでしょうか。

このスプラインの次数は 4 ですが、次の節点シーケンス内の両端の節点を見ると、

knots

knots =

     1     1     2     3     4     5     6     7     8     9     9

多重度は 2 です。したがって、この節点シーケンスの次数 4 のすべての B スプラインは、基本区間の端点で 0 となります。このため、この曲線は (0,0) で開始および終了することになります。

修復

この場合、実際に関心があるのはパラメーター区間 [3 .. 7] に対応する曲線セグメントのみであるため、fnbrk を使用してその部分を抽出すると、fnplt でそれを容易にプロットできます (黄色)。

mycv = fnbrk(curve,[3 7]);
fnplt(mycv,'y',2.5);
title('The Spline Curve of Interest, in Yellow')

この曲線で囲まれた部分の面積

曲線 (のみ) を記述するスプライン、つまり mycv があるので、この閉じた曲線で囲まれた面積を、次のように簡単に計算できます。

area = diff(fnval(fnint( ...
       fncmb(fncmb(mycv,[0 1]),'*',fnder(fncmb(mycv,[1 0]))) ...
                        ),fnbrk(mycv,'interval')))

area =

   -0.8333

少し工夫すると、これを次の積分の値として認識できます。

  int y(t) d(x(t)) = int y(t) Dx(t) dt

このとき、スプライン mycv の基本区間で、(x(t),y(t)) := fnval(mycv,t) をパラメーター値 t に対応する曲線上の点とします。ここで、fncmb(mycv,[1,0])、fncmb(mycv,[0,1]) はスプライン曲線の 2 つの成分、すなわちスカラー値のスプライン x と y を表します。

また、この曲線は半径 1/2 の円に近くなっています。そのため、面積をおよそ次のように見積もることができます。

disp(pi/4)

    0.7854

計算された面積が "負" になるのはなぜでしょうか。その理由としては、t を増やしつつ曲線上を移動すると、曲線で囲まれた面積が左側にくるからです。これを確認するため、いくつかの接ベクトルを描画します。

接ベクトルの追加

もう一度曲線を描画し、何点かにおいてその曲線に対する接ベクトルを描きます。

hold off
fnplt(mycv,'y',2.5); hold on
t = 3:.4:6.2;
cv = fnval(curve, t);
cdv = fnval(fnder(curve), t);
quiver(cv(1,:),cv(2,:), cdv(1,:),cdv(2,:));
title('A Spline Curve With Some Tangents')
axis([-.4 1.4 -.2 1.2]), axis equal

曲線と直線の交差部

このスプライン曲線と直線 y = x が交差する点を決定する場合は、次のコードでその点を求め、曲線の内部にその直線のセグメントをプロットします。

cuts = fnval(mycv, ...
    mean(fnzeros(fncmb(fncmb(mycv,[0,1]),'-',fncmb(mycv,[1,0])))));
plot(cuts(1,:), cuts(2,:),'y','LineWidth',2.5)
hold off
title('A Spline Curve With Some Tangents and a Cut Across')

SPCRV: 制御多角形および対応するスプライン曲線

スプライン曲線は、ある形状を大ざっぱに仮定する滑らかな曲線があればよい図を生成する場合に、広く使用されています。このため、Curve Fitting Toolbox には、ツールボックスの他の部分とは独立して使用できる特別なコマンド spcrv が含まれています。

平面内の点シーケンスおよびオプションとして次数 k を指定すると、spcrv は中間点での節点挿入を繰り返し、指定したシーケンスにより制御多角形が指定される、次数 k のスプライン曲線を生成します。

次の Figure に、そのような制御多角形とそれに対応する次数 3 のスプライン曲線を示します。

points = [0 0; 1 0; 1 1; 0 2; -1 1; -1 0; 0 -1; 0 -2].';
values = spcrv(points,3);

plot(points(1,:),points(2,:),'k');
axis([-2 2.25 -2.1 2.2]);
hold on
plot(values(1,:),values(2,:),'r','LineWidth',1.5);
legend({'Control Polygon' 'Quadratic Spline Curve'},  'location','SE');

この曲線は制御多角形の各セグメントの中点で接し、制御多角形の描く形状に沿うことに注意してください。

次数の増加

次数 k を大きくすると制御多角形から曲線が離れ、曲線が平滑化かつ短縮化されます。ここでは、次数 4 の対応するスプライン曲線を追加しました。

value4 = spcrv(points,4);
plot(value4(1,:),value4(2,:),'b','LineWidth',2);
legend({'Control Polygon' 'Quadratic Spline Curve' ...
        'Cubic Spline Curve'}, 'location','SE');

CSCVN

一方、内挿曲線を求めるには、'自然な' パラメトリック 3 次スプライン曲線を提供する cscvn コマンドを使用します。

fnplt(cscvn(points), 'g',1.5);
legend({'Control Polygon' 'Quadratic Spline Curve' ...
        'Cubic Spline Curve' 'Interpolating Spline Curve'}, ...
        'location','SE');

2 番目の制御点 (1,0) 付近の点 (.95,-.05) を追加すると、そこでより急に曲がる内挿スプライン曲線を作成できます。

np = size(points, 2);
fnplt( cscvn([ points(:,1) [.95; -.05] points(:,2:np) ]), 'm',1.5);
plot(.95,-.05,'*');
legend({'Control Polygon' 'Quadratic Spline Curve' ...
        'Cubic Spline Curve' 'Interpolating Spline Curve' ...
        'Faster Turning Near (1,0)'}, ...
        'location','SE');
hold off

RSCVN

平面の任意の点シーケンスを通過し、オプションでその点における任意の法線方向に直交する、円弧で構成された接線連続曲線を求めることもできます。コマンド rscvn がその曲線を提供します。

たとえば、次により、円が生成されます。

c = rscvn([-1 1 -1;0 0 0],[1 1;0 0]);

次のようにプロットします。

fnplt(c);
axis([-1.05 1.05 -1.05 1.05]), axis equal, axis off

c は、2 つの部分のみで構成される 2 次有理スプラインです。次のコマンドを見るとはっきりします。

[form, order, breaks] = fnbrk(c,'f','o','b')

form =

    'rBform'


order =

     3


breaks =

     0     2     4

このツールを使用すると、数点のデータ点だけで印象的なパターンを簡単に生成できます。たとえば、次に、ベルファストのアルスター博物館にある三脚ともえ紋の青銅メダルのデザインを変えたものを示します。これは、はるか昔に円弧の断片によりデザインされたと考えられています。

pp =[zeros(1,7); 5.4, 3, 6.9, 2.75, 2.5, .5, 5];
alpha = 2*pi/3; ca = cos(alpha); sa = sin(alpha); c = [ca sa;-sa ca];
d = [0 0 .05 -.05;1 -1 .98 .98]; d = [d c*d];
yin = rscvn([pp(:,[7,1:3]),c*pp(:,3:4),pp(:,3)], d(:,[1 2 1 4 7 5 1]));
fnplt(yin), hold on, fnplt(fncmb(yin,c)), fnplt(fncmb(yin,c'))
yang = rscvn([pp(:,6),-pp(:,6),pp(:,5),c*pp(:,4)],[d(:,[2 1 1]),c(:,2)]);
fnplt(yang), fnplt(fncmb(yang,c)), fnplt(fncmb(yang,c'))
axis([-7.2 7.2 -7.2 7.2]), axis equal, axis off, hold off