今回はB-spline基底関数とB-spline曲線をMATLABで描写していきます。
- B-spline基底関数、B-spline曲線の具体的なコードが知りたい
- コード例を見てMATLABの組み込み関数について理解したい
B-spline曲線の概要については前回の記事で詳しく解説しているので、B-spline曲線の概要を理解していない方はまずはそちらからご覧ください。
MATLABコード公開
B-spline基底関数
まずはB-spline曲線を構成する基底関数であるB-spline基底関数のコードを公開します。
おさらいとして、B-spline基底関数は以下の式で表されます。
これをMATLABで新規関数として定義します。以下のコードをコピーしてmファイルとして保存して下さい。
保存する際のファイル名は以下の【】に示すように「bsplinebasisfunction.m」として下さい。
【bsplinebasisfunction.m】
function N = bsplinebasisfunction(i,p,xi,knot)
if p == 0
if knot(i) <= xi && xi < knot(i+1)
N = 1;
else
N = 0;
end
else
Np1 = bsplinebasisfunction(i,p-1,xi,knot);
Np2 = bsplinebasisfunction(i+1,p-1,xi,knot);
if knot(i+p) == knot(i) % define 0/0 = 0
fterm = 0;
else
fterm = ( (xi-knot(i)) / (knot(i+p)-knot(i)) ) * Np1;
end
if knot(i+p+1) == knot(i+1) % define 0/0 = 0
sterm = 0;
else
sterm = ( (knot(i+p+1)-xi) / (knot(i+p+1)-knot(i+1)) ) * Np2;
end
N = fterm + sterm;
end
end後で解説しますが、上記のコードのみだと正しい結果が出ないため、もうひとつ、補正用のコードとして「bsplinebasisfunction1.m」も併せて保存しておいて下さい。
【bsplinebasisifunction1.m】
function N = bsplinebasisfunction1(i,p,xi,knot)
if p == 0
if knot(i) <= xi && xi <= knot(i+1)
N = 1;
else
N = 0;
end
else
Np1 = bsplinebasisfunction1(i,p-1,xi,knot);
Np2 = bsplinebasisfunction1(i+1,p-1,xi,knot);
if knot(i+p) == knot(i) % define 0/0 = 0
fterm = 0;
else
fterm = ( (xi-knot(i)) / (knot(i+p)-knot(i)) ) * Np1;
end
if knot(i+p+1) == knot(i+1) % define 0/0 = 0
sterm = 0;
else
sterm = ( (knot(i+p+1)-xi) / (knot(i+p+1)-knot(i+1)) ) * Np2;
end
N = fterm + sterm;
end
end「bspinebasisfunction.m」と「bspinebasisfunction1.m」の違いは、3行目の if 文の後半が「xi < knot(i+1)」なのか「xi <= knot(i+1)」なのかの違いのみです。
メインプログラム
次にB-spline基底関数およびB-spline曲線のメインプログラムを公開します。
こちらもおさらいとして、B-spline曲線は以下の式で表されます。
メインプログラムのファイル名は何でも大丈夫です。
自分の分かりやすい名前をつけて、先ほど公開したB-spline基底関数のコードと同じファイル内に以下のコードを保存して下さい。
今回はB-spline基底関数 Ni,p(ξ) とB-spline曲線 C(ξ) のメインプログラムなので「main_N_C.m」とします。
【main_N_C.m】
clear
n = 8; % 制御点数 n
p = 2; % 曲線次数 p
knot = [0 0 0 1 2 3 4 4 5 5 5]; % ノットベクトル
xi = linspace(0,5,500); % size 500
P = [ 3 3 0; 6 1 0; 9 1 0;
9 5 0; 12 5 0; 16.5 9.5 0;
9 13 0; 6 8 0 ];
N = zeros( n, length(xi) );
Cx0 = zeros( size(xi) ); Cy0 = zeros( size(xi) );
Cx = zeros( size(xi) ); Cy = zeros( size(xi) );
for j = 1:length(xi)
Cx(j) = 0;
Cy(j) = 0;
for i = 1:n
if j == length(xi)
N(i,j) = bsplinebasisfunction1(i,p,xi(j),knot);
else
N(i,j) = bsplinebasisfunction(i,p,xi(j),knot);
end
Cx0(j) = N(i,j) * P(i,1);
Cy0(j) = N(i,j) * P(i,2);
Cx(j) = Cx(j) + Cx0(j);
Cy(j) = Cy(j) + Cy0(j);
end
end
% B-spline基底関数をプロットする
figure(1); cla
plot(xi,N,'LineWidth',1.5)
axis equal;
axis([0 5 0 1])
xlabel('ξ','FontSize',20); ylabel('N(ξ)','FontSize',20);
title('B-Spline Basis Function','FontSize',20)
set(gca,'FontSize',20);
hold off
% B-spline曲線をプロット
figure(2); cla
% 制御多角形をプロット
plot(P(:,1),P(:,2),'k');
hold on
% 制御点をプロット
plot(P(:,1),P(:,2),'ro','markerFaceColor','red','markersize',4);
% 曲線をプロット
plot(Cx,Cy,'b','LineWidth',2)
axis([0 18 0 15])
% 各軸の名前と図の凡例を追記
xlabel('x','FontSize',18)
ylabel('y','FontSize',18)
title('B-Spline Curve','FontSize',18)
set(gca,'FontSize',18);
hanrei = legend('Control Polygon','B-Spline Curve','Conrtol Points');
set(hanrei,'Location','NorthWest','FontSize',10)
hold off公開コードの解説
各変数の設定
メインプログラム冒頭の3-5行目で、B-spline基底関数を構成する以下の3つの変数を定義します。
- 制御点の数 n
- B-spline基底関数（B-spline曲線）の次数 p
- ノットベクトル Ξ
制御点の数 n は描写したい図の通り8個、曲線の次数 p は2次とします。
ノットベクトル Ξ については以下のルールを再確認しておきます。
- ノットベクトルの長さは n+p+1 とする。
- ノットベクトルの値は単純増加（値が減少しない）とする。
- 多重度を導入する場合を除き、増加間隔は等間隔とする。
- 領域の端部（始点と終点）はノットを p+1個重複（＝多重度 p+1）させる。
- 制御点を通過させたい場合はノットを p個重複（＝多重度 p）させる。
これらのルールを踏まえ、ノットベクトル Ξ = { 0, 0, 0, 1, 2, 3, 4, 4, 5, 5, 5 } としています。
MATLAB組み込み関数：linspace
【4行目】
xi = linspace(0,5,500);ここではMATLABの組み込み関数「linspace」（linear space：線形空間）を使って ξ の値を定義しています。
「linspace」を使うことで、指定した区間を等分割することができます。
今回の場合は「0から5までの区間を500分割する」という意味となり、ξ のサイズは1行500列となります。
似たような記述方法として、「xi = 0:0.01:500」という方法もあります。
これは「0から0.01刻みで500まで並べる」という意味になるため、 ξ のサイズは1行501列となります。
- ベクトル ξ のサイズをキリ良くしたいなら「xi = linspace(0,5,500)」
- ベクトル ξ の各値をキリ良くしたいなら「xi = 0:0.01:500」
といった感じで使い分けると良いでしょう。
また、変数 ξ の間隔は0.01か0.02刻みくらいだとキレイな曲線が描けます。
変数のサイズを定義する
【12-14行目】
N = zeros( n, length(xi) );
Cx0 = zeros( size(xi) ); Cy0 = zeros( size(xi) );
Cx = zeros( size(xi) ); Cy = zeros( size(xi) );12行目から14行目にかけて、この後計算する変数のサイズを予め定義しています。
この工程は必須ではありませんが、変数のサイズを予め定義しておくことで、解析時間を短縮することができます。
forループが進むたびに変数のサイズが変わる計算はCPUに負担がかかり、計算時間が長くなってしまいます。
最初に変数のサイズを指定（このとき、各値はすべて0）し、各値を計算して代入していくほうが計算速度が速くなります。
補正用コード：bsplinebasisifunction1.m
【16-20行目】
for j = 1:length(xi)
Cx(j) = 0;
Cy(j) = 0;
for i = 1:n
if j == length(xi)
N(i,j) = bsplinebasisfunction1(i,p,xi(j),knot);
else
N(i,j) = bsplinebasisfunction(i,p,xi(j),knot);
end
Cx0(j) = N(i,j) * B(i,1);
Cy0(j) = N(i,j) * B(i,2);
Cx(j) = Cx(j) + Cx0(j);
Cy(j) = Cy(j) + Cy0(j);
end
end先ほど公開した補正用のコードについて説明します。
実際に「bsplinebasisfunction.m」のみを用いた場合と、補正用コードを用いた場合のB-spline基底関数を比較します。
左図と右図の違い、分かりますか？
実はCox-de Boorの式をそのままコードとして書くと、左図のようなB-spline基底関数が返ってきます。
左図のB-spline基底関数のマズいところは、ξ = 5 のときの N(ξ) の値が 0 になってしまっているところです。
これに気付かずにB-spline曲線を描くと、下の左図のようにセグメントの最後の値が 0 になってしまいます。
補正関数を用いた正しい結果図が下の右図となります。
まとめ
今回はB-spline基底関数とB-spline曲線のMATLABコード例を公開し、コードについて解説しました。
B-spline基底関数はデザイン分野で最もよく使用されるNURBS関数やIsogeometric Analysisの基本となる非常に汎用性の高い関数です。
また、MATLABは結果の図示方法が非常に充実しています。
すべては紹介しきれませんが、これまで公開してきたコードの中で様々な使用例を載せていますので参考にしてください。
最後までお読みいただきありがとうございました。
コメント