B-splineとNURBSの違い：重み関数を徹底解説（MATLABコード公開）

前回の記事ではB-spline曲線と、それを構成するB-spline基底関数について解説しました。

今回はB-spline曲線により汎用性を持たせたNURBS（ナーブス）曲線について解説します。

B-spline曲線と何が違うのか等、徹底解説していきますので、ぜひお読み下さい。

この記事はこんな人にオススメ

NURBS曲線について知りたい
B-spline曲線とNURBS曲線の違いが分からない
NURBS曲線に出てくる「重み」って何？

NURBSとは
B-splineからNURBSへ拡張する重み関数とは
1. 違い①：NURBSは、ノットベクトルの間隔が一様でなくてもよい
2. 違い②：NURBSでは、各制御点が”重み”を持つ
NURBS 基底関数
NURBS曲線
制御点の重みによる曲線形状の違い
MATLABコード公開
まとめ

NURBSとは

NURBSとは、Non-Uniform Rational B–Spline（非一様有理B-スプライン）の略称です。

NURBSは従来のベジェ（Bezier）やB-splineなどの曲線・曲面をより一般化したもので、円・円弧・楕円・楕円弧などの円錐曲線や球面・円柱面・円錐面などの曲面を近似なしに厳密に表現できるように改良されたものです。

さらに制御点と曲線・曲面との関係を拡張し、複数の制御点が重なることを許容するなど、より自由度の高い関数の表現が可能となっています。

以上のようなことから、今日のほとんどのCAD/CAEシステムでNURBSが使われており、形状モデリングを行う上でのスタンダードとなっています。

トライ

NURBSで厳密な円を表現した例を以下に示します。

NURBS曲線の次数 p = 2

ノットベクトル Ξ = { 0, 0, 0, 1, 1, 2, 2, 3, 3, 4, 4, 4 }

制御点の重み w = {1.0, √2/2, 1.0, √2/2, 1.0, √2/2, 1.0, √2/2, 1.0 }

Isogeometric Analysis NURBS — NURBSによる厳密円

B-splineからNURBSへ拡張する重み関数とは

これまでも取り上げてきたB-splineをNURBSに拡張する方法について説明します。

B-splineとNURBSでは具体的に以下の2点が異なります。

B-splineとNURBSの違い

違い①：NURBSは、ノットベクトルの間隔が一様でなくてもよい
違い②：NURBSでは、各制御点が”重み”を持つ

違い①：NURBSは、ノットベクトルの間隔が一様でなくてもよい

B-splineは始点と終点および多重度を導入する場合を除き、ノットベクトルの増加間隔は一定でなければなりません。

しかし、NURBSにはそのような制約はなく、ノットベクトルの値、つまりセグメントの位置を自由に設定することができます。

これがNURBS（Non-Uniform Rational B-Spline）のうち、「Non-Uniform」に該当する部分です。

ですが、これはそれほど重要ではないので「ふーん、そうなんだ」くらいで大丈夫です。

トライ

ノットベクトルの多重度の考え方やセグメントといった単語の意味については以前の記事で解説しているので、分からない人はこちらもあわせてお読み下さい。

違い②：NURBSでは、各制御点が”重み”を持つ

B-splineとNURBSの最も大きな違いであり、今回のメインテーマでもある違いの2つ目は

「NURBSでは、各制御点が”重み”を持つ」ということです。

制御点が重みを持つことで、Rational（有理な＝分数で表すことができる）B-Splineとなります。

重み関数 W (ξ) は、B-spline基底関数N_i,_p(ξ) に制御点の重みw_i を掛けることで次の式で表されます。

NURBS 基底関数

NURBS基底関数は、B-spline基底関数N_i,_p(ξ) に制御点の重みw_i を掛けたものを重み関数 W (ξ) で割ることで以下の式で表されます。

各制御点の重み w_i は，0から∞の範囲で各制御点ごとに自由に設定することができます。

トライ

ちなみに制御点の重み wi をすべて 1 とすると、NURBS基底関数とB-spline基底関数はまったく同じものとなります。

NURBS曲線

NURBS曲線は、NURBS基底関数 R(ξ) と制御点座標 P_i を用いて、以下の式で表されます。

制御点の重みによる曲線形状の違い

NURBS曲線は制御点座標 P_i と重み w_i の値により、形状を自由に変化させることができます。

例として、ある制御点の重みだけを変更した場合、NURBS曲線の形状がどのように変化するのかを見てみます。

制御点や曲線の次数などは以下のように設定しました。

制御点数 n = 4
NURBS曲線の次数 p = 3
ノットベクトル Ξ = { -1, -1, -1, -1, 1, 1, 1, 1 }
各制御点の重み w：制御点2の重みを1, 2, 5と変化させる（それぞれcase A～Cとする）
制御点2以外の制御点の重みはすべて1とする

では実際にNURBS基底関数とNURBS曲線を見てみましょう。

こうやって比較してみると違いが一目瞭然ですね。

制御点2の重みを増やしていくと、制御点2に対応する基底関数（左図でいうオレンジの関数）の値が大きくなり、NURBS曲線も重みの大きい制御点2の近くを通過する形状になります。

逆に重み wi の値を小さくすると、それに対応するNURBS基底関数の値は小さくなり、NURBS曲線はその制御点から遠ざかります。

この性質を利用すれば、制御点座標を変更することなく局所的に曲線の形状を変更することができます。

トライ

ちなみにある制御点の重みを 0 にすると、その制御点が存在しないものとしてNURBS曲線が描かれます。

実際に制御点2の重みを0にした場合を見てみます。

下の図のように、制御点2に対応するオレンジの基底関数の値が全範囲で0、NURBS曲線は制御点1, 3, 4しかない場合の曲線形状となっています。

MATLABコード公開

最後に、先ほどお見せしたfigure(1)～figure(6)を作図するためのMATLABコードを無料公開します。

以下のコードをコピーし、MALTABもしくはOctaveのエディタ画面にペーストして保存して下さい。

保存するときの mファイルの名前は何でもOKです。

トライ

オススメの名前の付け方としては、以下のようなメインプログラムは「main_NURBS.m」といった感じで頭に「main」とつけておくと後で分かりやすいかと思います。

【main_NURBS.m】

clear

P = [ 0  1  0; 8  6  0; 11  0  0; 18  3  0 ]; % 制御点座標
n = 4; % 制御点数
p = 3; % NURBS曲線の次数
xi = -1:0.01:1; % 局所座標刻み
knot = [ -1 -1 -1 -1  1  1  1  1 ]; % ノットベクトル
wa = [ 1  1  1  1 ]; % 重み case A
wb = [ 1  2  1  1 ]; % 重み case B
wc = [ 1  5  1  1 ]; % 重み case C
% N:B-spline基底関数
N = zeros( n, length(xi) );
for i = 1:n
    for j = 1:length(xi)
        if j == length(xi)
            N(i,j) = bsplinebasisfunction1(i,p,xi(j),knot);
        else
            N(i,j) = bsplinebasisfunction(i,p,xi(j),knot);
        end
    end
end

% W:NURBSの重み関数  R:NURBS基底関数
Wa0 = zeros( size(xi) );  Wa = zeros( size(xi) );  Ra = zeros( n, length(xi) );
Wb0 = zeros( size(xi) );  Wb = zeros( size(xi) );  Rb = zeros( n, length(xi) );
Wc0 = zeros( size(xi) );  Wc = zeros( size(xi) );  Rc = zeros( n, length(xi) );
for j = 1:length(xi)
    for i = 1:n
        Wa0(j) = N(i,j) * wa(i);   Wa(j) = Wa(j) + Wa0(j);
        Wb0(j) = N(i,j) * wb(i);   Wb(j) = Wb(j) + Wb0(j);
        Wc0(j) = N(i,j) * wc(i);   Wc(j) = Wc(j) + Wc0(j);
    end
    for i = 1:n
        if Wa(j) == 0
            Ra(i,j) = 0;
        else
            Ra(i,j) = ( N(i,j) * wa(i) ) / Wa(j);
        end
        if Wb(j) == 0
            Rb(i,j) = 0;
        else
            Rb(i,j) = ( N(i,j) * wb(i) ) / Wb(j);
        end
        if Wc(j) == 0
            Rc(i,j) = 0;
        else
            Rc(i,j) = ( N(i,j) * wc(i) ) / Wc(j);
        end
    end
end
% NURBS基底関数を作図
figure(1)
plot(xi,Ra,'LineWidth',1.5);
axis equal;
axis([-1 1 0 1])
xlabel('ξ','FontSize',18)
ylabel('R(ξ)','FontSize',18)
title('NURBS Basis Function(case A)','FontSize',18)
set(gca,'FontSize',18);

figure(2)
plot(xi,Rb,'LineWidth',1.5);
axis equal
axis([-1 1 0 1])
xlabel('ξ','FontSize',18)
ylabel('R(ξ)','FontSize',18)
title('NURBS Basis Function(case B)','FontSize',18)
set(gca,'FontSize',18);

figure(3)
plot(xi,Rc,'LineWidth',1.5);
axis equal
axis([-1 1 0 1])
xlabel('ξ','FontSize',18)
ylabel('R(ξ)','FontSize',18)
title('NURBS Basis Function(case C)','FontSize',18)
set(gca,'FontSize',18);

% NURBS曲線を作図
NCx0a = zeros( size(xi) );  NCxa = zeros( size(xi) ); % case A
NCy0a = zeros( size(xi) );  NCya = zeros( size(xi) );
NCx0b = zeros( size(xi) );  NCxb = zeros( size(xi) ); % case B
NCy0b = zeros( size(xi) );  NCyb = zeros( size(xi) );
NCx0c = zeros( size(xi) );  NCxc = zeros( size(xi) ); % case C
NCy0c = zeros( size(xi) );  NCyc = zeros( size(xi) );
for j = 1:length(xi)
    for i = 1:n
        NCx0a(j) = Ra(i,j) * P(i,1); % case A
        NCy0a(j) = Ra(i,j) * P(i,2);
        NCxa(j) = NCxa(j) + NCx0a(j);
        NCya(j) = NCya(j) + NCy0a(j);
        
        NCx0b(j) = Rb(i,j) * P(i,1); % case B
        NCy0b(j) = Rb(i,j) * P(i,2);
        NCxb(j) = NCxb(j) + NCx0b(j);
        NCyb(j) = NCyb(j) + NCy0b(j);
        
        NCx0c(j) = Rc(i,j) * P(i,1); % case C
        NCy0c(j) = Rc(i,j) * P(i,2);
        NCxc(j) = NCxc(j) + NCx0c(j);
        NCyc(j) = NCyc(j) + NCy0c(j);
    end
end

figure(4) % case A
plot(NCxa,NCya,'b','LineWidth',1.5); % NURBS曲線
hold on
plot(P(:,1),P(:,2),'k'); % 制御多角形
plot(P(:,1),P(:,2),'ro','LineWidth',1.5,'MarkerSize',4,'MarkerFacecolor','red'); % 制御点
axis equal % 図のX軸とY軸の比率を等しくする
axis([ -1 19 -1 10 ]) % X軸とY軸の範囲を指定
xlabel('x','FontSize',18)
ylabel('y','FontSize',18)
title('NURBS Curve(case A)','FontSize',18)
set(gca,'FontSize',18);
hanrei = legend('NURBS Curve','Control Polygon','Control Points');
set(hanrei,'Location','NorthEast','FontSize',10)
hold off

figure(5) % case B
plot(NCxb,NCyb,'b','LineWidth',1.5);
hold on
plot(P(:,1),P(:,2),'k');
plot(P(:,1),P(:,2),'ro','LineWidth',1.5,'MarkerSize',4,'MarkerFacecolor','red');
axis equal
axis([ -1 19 -1 10 ])
xlabel('x','FontSize',18)
ylabel('y','FontSize',18)
title('NURBS Curve(case B)','FontSize',18)
set(gca,'FontSize',18);
hanrei = legend('NURBS Curve','Control Polygon','Control Points');
set(hanrei,'Location','NorthEast','FontSize',10)
hold off

figure(6) % case C
plot(NCxc,NCyc,'b','LineWidth',1.5);
hold on
plot(P(:,1),P(:,2),'k');
plot(P(:,1),P(:,2),'ro','LineWidth',1.5,'MarkerSize',4,'MarkerFacecolor','red');
axis equal
axis([ -1 19 -1 10 ])
xlabel('x','FontSize',18)
ylabel('y','FontSize',18)
title('NURBS Curve(case C)','FontSize',18)
set(gca,'FontSize',18);
hanrei = legend('NURBS Curve','Control Polygon','Control Points');
set(hanrei,'Location','NorthEast','FontSize',10)
hold off

figure(7) % case A-C 重ね合わせ
plot(NCxa,NCya,'b','LineWidth',1.5); hold on
plot(NCxb,NCyb,'b','LineWidth',1.5);
plot(NCxc,NCyc,'b','LineWidth',1.5);
plot(P(:,1),P(:,2),'k');
plot(P(:,1),P(:,2),'ro','LineWidth',1.5,'MarkerSize',4,'MarkerFacecolor','red');
axis equal
axis([-1 19 -1 7 ])
xlabel('x','FontSize',18)
ylabel('y','FontSize',18)
title('NURBS Curve(case A-C)','FontSize',18)
set(gca,'FontSize',18);
hold off

このメインプログラムのでは、以下の2つを新規関数として用いています。

bsplinebasisfunction.m
bsplinebasisfunction1.m

いずれもB-spline基底関数を計算するための新規関数です。

以前の記事で数式とコードを公開していますが、この記事が初見の方は以下の2つのコードをコピーし、エディタ画面にペーストして保存して下さい。

ただし、保存するときの mファイルの名前は【】内と同じにして下さい。

B-spline基底関数を求める新規関数1【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

B-spline基底関数を求める新規関数2【bsplinebasisfunction1.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

まとめ

今回はB-splineからNURBSの最も大きな違いである「重み関数」について解説しました。

今回の内容をまとめるとこのようになります。

今回の内容まとめ

NURBSはノットベクトルの値が非一様（Non-Uniform）でよい
NURBSは各制御点が重みを持つ、有理（Rational）なB-splineである
NURBSは各制御点の重みを変更することで、制御点座標を変更することなく局所的に曲線の形状を変更することができる

CADやRhinocerosを使っていると当たり前に使っているNURBS曲線ですが、どういった理論で構成されているかを理解することは重要です。

今回の記事を読んで、

「制御点の重みを変えれば、わざわざ制御点を動かさなくても局所的に曲線形状を変えることができる」

ということを知っておけば、より自由な形状をスタディできるようになることでしょう。

今後も有限要素法（FEM）や Isogeometric Analysis（IGA）の理論解説とMATLABコードの無料公開を続けていきますので、今回の記事が「役に立った」と思っていただけた方は、ぜひ他の記事もご覧いただけると幸いです。

最後までお読みいただきありがとうございました。